Finite-temperature topological invariant for interacting systems

TL;DR

This paper extends the Ensemble Geometric Phase to finite-temperature interacting 1D systems, establishing a topological invariant that remains quantized below a critical temperature and demonstrating its robustness through numerical simulations of a Bose-Hubbard model.

Contribution

It introduces a finite-temperature topological invariant for interacting 1D systems that generalizes previous zero-temperature invariants, applicable even with degeneracies and fractional fillings.

Findings

The finite-temperature invariant matches the zero-temperature invariant below a critical temperature.

Numerical simulations show the invariant remains quantized despite non-quantized particle transport at higher temperatures.

The approach applies to systems with degenerate ground states and fractional fillings.

Abstract

We generalize the Ensemble Geometric Phase (EGP), recently introduced to classify the topology of density matrices, to finite-temperature states of interacting systems in one spatial dimension (1D). This includes cases where the gapped ground state has a fractional filling and is degenerate. At zero temperature the corresponding topological invariant agrees with the well-known invariant of Niu, Thouless and Wu. We show that its value at finite temperatures is identical to that of the ground state below some critical temperature $T_c$ larger than the many-body gap. We illustrate our result with numerical simulations of the 1D extended super-lattice Bose-Hubbard model at quarter filling. Here a cyclic change of parameters in the ground state leads to a topological charge pump with fractional winding $\nu=1/2$. The particle transport is no longer quantized when the temperature becomes comparable to the many-body gap, yet the winding of the generalized EGP is.

Figures

Equations

$$\hat{U}(\theta )=e^{i\theta \hat{X}},\text{with}\hat{X}=\frac{1}{L}\sum _{j=1}^L\sum _{s=1}^n\left(j+r_s\right){\hat{n}}_{js}$$

$$\Delta n$$

$${\nu }_{\text{tot}}=\frac{1}{2\pi }{\int }_0^{\mathcal{T}}dt\frac{\partial }{\partial t}\text{Im}\log \det \mathsf{W}(t),\mathsf{W}=\mathcal{P}{e}^{i{\int }_0^{2\pi }d\theta \mathsf{A}(\theta )}.$$

$${\nu }_{\text{tot}}=\frac{1}{2\pi }{\int }_0^{\mathcal{T}}dt{\int }_0^{2\pi d}d\theta \text{Im}⟨{\partial }_t{\Psi }_0^{\mu }|{\partial }_{\theta }{\Psi }_0^{\mu }⟩.$$

$$\nu_{\textrm{tot}}=\oint\frac{d\lambda}{2\pi}\frac{\partial\phi^{(d)}_{\mu}}{\partial\lambda},\quad\textrm{with}\enspace\phi^{(d)}_{\mu}=\textrm{Im}\log\bigl{\langle}\Psi_{0}^{\mu}\bigr{|}({\hat{U}})^{d}\bigl{|}\Psi_{0}^{\mu}\bigr{\rangle},$$

$$\phi_{\textrm{EGP}}^{(d)}=\textrm{Im}\log\textrm{Tr}\bigl{\{}\rho\,(\hat{U})^{d}\bigr{\}}$$

$$H=-t_1\sum _{j,\text{even}}{\hat{a}}_j^{†}{\hat{a}}_{j+1}-t_2\sum _{j,\text{odd}}{\hat{a}}_j^{†}{\hat{a}}_{j+1}+h.a.$$

$$+\frac{\Delta }{2}\sum _j(-1{)}^j{\hat{a}}_j^{†}{\hat{a}}_j+\frac{U}{2}\sum _j{\hat{n}}_j({\hat{n}}_j-1)+\sum _{j,d}{V}_d{\hat{n}}_i{\hat{n}}_{i+d}.$$

$$\Delta {\varphi }_{\text{EGP}}=2\pi \mathrm{deg}z=\mathrm{Im}\oint \frac{1}{z\left(\lambda \right)}\frac{dz\left(\lambda \right)}{d\lambda }d\lambda$$

$$\bigl{|}z(T_{c},\lambda)\bigr{|}=0,$$

$$T_0\sim {\Delta }_{\text{gap}}/\ln L.$$

$$\tilde{H}=E_0\sum _{\mu =1}^d|{\Psi }_0^{\mu }⟩⟨{\Psi }_0^{\mu }|+E_1\sum _{j\ne (0,\mu )}|{\Psi }_j⟩⟨{\Psi }_j|,$$

$$T_c=\frac{{\Delta }_{\text{gap}}}{\ln \left(1+|\text{Tr}\{(\hat{U}{)}^d\}|/d\right)}.$$

$$\hat{U}=e^{i2\pi \hat{X}},\text{with}\hat{X}=\frac{1}{L}\sum _{j=1}^L\sum _{s=1}^n\left(j+r_s\right){\hat{n}}_{js},$$

$$\hat{U}=\prod _{j=1}^L{\hat{U}}_j=\prod _{j=1}^L{e}^{i2\pi {\hat{X}}_j},\text{with}{\hat{X}}_j=\frac{1}{L}\sum _{s=1}^n\left(j+r_s\right){\hat{n}}_{js}.$$

$${\rho }_{\mathrm{local}}[{\bar{m}}_L,{\bar{\sigma }}_L,{\bar{\sigma }}_R,m_L,{\sigma }_L,{\sigma }_R]=\sum _{m_R}{C}^{\ast }[{\bar{m}}_L,{\bar{\sigma }}_L,{\bar{\sigma }}_R,m_R]C[m_L,{\sigma }_L,{\sigma }_R,m_R]$$

$$\begin{array}{rl}(\rho U_0{)}_{\mathrm{extend}}[{\bar{m}}_L^{\prime },{\bar{\sigma }}_L^{\prime \prime },{\bar{\sigma }}_L^{\prime },m_L^{\prime },{\sigma }_L^{\prime \prime },{\sigma }_L^{\prime }]=\sum _{{\bar{m}}_L,m_L,{\bar{\sigma }}_L,{\bar{\sigma }}_R,{\sigma }_R,{\sigma }_L}& L^{\ast }[{\bar{m}}_L^{\prime },{\bar{\sigma }}_L^{\prime \prime },{\bar{\sigma }}_L^{\prime },{\bar{m}}_L]L[m_L^{\prime },{\sigma }_L^{\prime \prime },{\sigma }_L^{\prime },m_L]\times \\ & \times {\rho }_{\mathrm{local}}[{\bar{m}}_L,{\bar{\sigma }}_L,{\bar{\sigma }}_R,m_L,{\sigma }_L,{\sigma }_R]U_0[{\bar{\sigma }}_L,{\bar{\sigma }}_R,{\sigma }_L,{\sigma }_R].\end{array}$$

$$\text{Tr}\{(\hat{U}{)}^d\}=\text{Tr}\left\{\exp \left(\frac{2\pi i}{L}d\sum _{j=1}^L\sum _s(j+r_s){\hat{n}}_{js}\right)\right\}.$$

$$\text{Tr}\{(\hat{U}{)}^d\}=\frac{1}{2\pi }{\int }_0^{2\pi }d\theta e^{iN\theta }\prod _{j,s}\sum _n{e}^{-in(\theta -{\varphi }_s)}\exp \left\{\frac{2\pi i}{N}nj\right\}$$

$$\text{Tr}\{(\hat{U}{)}^d\}=\frac{1}{2\pi }{\int }_0^{2\pi }d\theta e^{iN\theta }\prod _s\left[\prod _{j=1}^L\left(1+e^{-i(\theta -{\varphi }_s)}{e}^{2\pi ij/N}\right)\right]$$

$$=\frac{1}{2\pi }{\int }_0^{2\pi }d\theta e^{iN\theta }\prod _s{\left[\prod _{j=1}^N\left(1+e^{-i(\theta -{\varphi }_s)}{e}^{2\pi ij/N}\right)\right]}^d$$

$$\prod _{j=1}^N\left(1+e^{-i\left(\theta -{\varphi }_s\right)}{e}^{2\pi ij/N}\right)=\prod _{j=1}^N\left(1-z{e}^{2\pi ij/N}\right)=$$

$$\prod _{j=1}^N\left(e^{-2\pi ij/N}-z\right)=1-z^N.$$

$$\displaystyle\textrm{Tr}\{(\hat{U})^{d}\}=\frac{1}{2\pi}\int_{0}^{2\pi}\!\!\!d\theta\,e^{iN\theta}\prod_{s}\Bigl{[}1-(-1)^{N}e^{-i\theta N}e^{i\phi_{s}N}\Bigr{]}^{d}$$

$$|\text{Tr}\{(\hat{U}{)}^d\}|=d\sum _s{e}^{2\pi i{r}_s}=\mathcal{O}(d).$$

$$z=\frac{⟨R⟩}{\mathcal{Z}},$$

$$⟨R⟩=\mathtt{tr}\left[\prod _{k=1}^N\exp \left(\beta {\overrightarrow{\alpha }}_k\overrightarrow{\sigma }\right)\right],$$

$$\mathcal{Z}=\mathtt{tr}\exp \left(-\beta H\right)=\prod _{k=1}^N\left(1+e^{\beta |{\overrightarrow{\alpha }}_k|}\right)\prod _{k=1}^N\left(1+e^{-\beta |{\overrightarrow{\alpha }}_k|}\right).$$

$$⟨R⟩=\mathtt{tr}\left[\prod _{k=1}^N\exp \left(\beta \left|{\overrightarrow{\alpha }}_k\right|{\sigma }_z\right)\cdot u_{k-1}{u}_k^{†}\right],$$

$$u_{k-1}{u}_k^{†}\approx 1+\frac{2\pi }{N}{u}_{k-1}\frac{\partial u_k^{†}}{\partial k}+...\underset{N>>1}{\to }1.$$

$$\left|⟨R⟩\right|\approx \exp \left(\sum _{k=1}^N\beta \left|{\overrightarrow{\alpha }}_k\right|\right)+\exp \left(-\sum _{k=1}^N\beta \left|{\overrightarrow{\alpha }}_k\right|\right)+\mathcal{O}\left(\frac{1}{N}\right).$$

Full text

Introduction. –

Starting with the discovery of the quantum Hall effect [1, 2, 3, 4, 5] topology has become an important paradigm for the understanding and classification of phases of matter [6, 7, 8]. Topology is characterized by integer-valued invariants which describe global properties of the system and are responsible for the robustness of characteristic features like quantized bulk transport, or edge states and edge currents. Invariants such as the winding of the geometric Zak phase or the Chern number of single-particle Bloch functions are defined in terms of the wave-function of ground states and are thus restricted to cases where the system is in a pure state.

In recent years several attempts have been made to generalize the concept of topology to finite-temperatures and to non-equilibrium steady states of non-interacting fermions [9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. For example, it was shown that the generalization of geometric phases to density matrices based on the Uhlmann construction [20] leads to consistent topological invariants in 1D [11, 12]. Its application to higher dimensions [13] is however faced with difficulties [21]. Other approaches predict an unphysical extensive number of topological phases for arbitrarily small temperatures [16].

Recently it was shown that the winding of the many-body polarization introduced by Resta [22] is an alternative topological invariant for Gaussian mixed states of fermions in 1D, termed ensemble geometric phase (EGP) [15, 17, 19], which can also be applied in 2D [23]. It was shown that the EGP of finite-temperature states in non-interacting systems is reduced to the ground-state Zak phase in the thermodynamic limit $L\to\infty$ and thus these states have the same topological classification as the corresponding ground states (following the Altland-Zirnbauer classification [24, 25, 26]). Despite being a genuine many-body quantity the EGP can be measured directly [17]. Furthermore a non-trivial winding of the EGP of a finite-temperature or non-equilibrium steady state upon cyclic parameter variations has direct physical consequences. E.g. it can lead to quantized transport in a weakly coupled auxiliary system initially prepared in a low temperature state [27]. For non-interacting bosons the EGP winding is always zero [28].

In the present paper we extend this concept to the case of interactions between particles including the possibility of interaction-induced fractionalization and degeneracy. We show that a generalization of the EGP to systems with a gapped ground state of fractional filling allows to define a topological invariant for finite-temperature states of one-dimensional systems of interacting particles. The winding of the EGP reduces to the well-known Niu-Thouless-Wu (NTW) invariant [29, 30] at $T=0$ and has the same value for all temperatures below a certain critical value. The EGP also provides a theoretical tool to detect topological order present in the ground state in cases where the gap is small or numerical calculations are restricted to non-zero temperatures. We illustrate our results with numerical simulations of the extended super-lattice Bose-Hubbard model (Ext-SLBHM) [31] at quarter filling, where the ground state is a doubly degenerate Mott insulator (MI).

Ensemble Geometric Phase: integer case. –

We consider 1D lattice models with a many-body hamiltonian $H(\lambda)$ that has a periodic dependence on some variable $\lambda$. We assume periodic boundary conditions in space with $L$ unit cells and allow for interactions among the particles. Since single-particle Bloch states $|u_{k}\rangle$ are no longer a good eigen-basis, topological invariants must be defined in terms of the many-body ground state $|\Phi_{0}\rangle$. As suggested by Niu, Thouless and Wu [30], the Zak phase of Bloch wavefunctions can be generalized by replacing the single-particle crystal momentum $k$ by a twist angle $\theta$ of boundary conditions $\Phi(x_{1},\dots x_{j}+L,\dots x_{N})=e^{i\theta}\Phi(x_{1},\dots x_{N})$. Since $\theta$ and $\theta+2\pi$ define the same boundary conditions, the parameter space $(\theta,\lambda)$ is a torus ${\sf T}^{2}$. Twisted boundary conditions can be removed and replaced by periodic ones via a canonical transformation to a twisted Hamiltonian, $\bar{H}(\theta)=\hat{U}(\theta)H\hat{U}(\theta)^{-1}$ and $|\Psi_{0}(\theta)\rangle=\hat{U}(\theta)|\Phi_{0}(\theta)\rangle$. Here

[TABLE]

is the momentum shift operator with $\hat{n}_{js}$ denoting the particle number operator of the $s$th site $(s\in\{1,2,\dots,n\})$ in the $j$th unit cell, and the lattice constant is $a=1$. $0\leq r_{s}\leq 1$ characterizes the position within the unit cell. In terms of the $|\Psi\rangle$’s the many-body equivalent of the Zak phase then reads \phi_{\textrm{MB}}=i\int_{0}^{2\pi}\!\!\!d\theta\,\,\bigl{\langle}\Psi_{0}(\theta)\bigr{|}\partial_{\theta}\Psi_{0}(\theta)\bigr{\rangle}.

If $|\Psi_{0}(\theta)\rangle$ is a gapped and non-degenerate ground state, a slow (adiabatic) change of the parameter $\lambda=\lambda(t)$ in a closed loop, such that the many-body gap does not close, induces a Thouless pump, described by a current density \langle\hat{j}\rangle=\partial_{\theta}E(\theta)/\hbar+i\bigl{(}\langle\partial_{t}\Psi_{0}|\partial_{\theta}\Psi_{0}\rangle-\langle\partial_{\theta}\Psi_{0}|\partial_{t}\Psi_{0}\rangle\bigr{)}. Then following Niu, Thouless and Wu, averaging over twisted boundary conditions one finds a strictly integer-quantized particle transport over one time period ${\cal T}$. The transported charge $\Delta n$ is then directly related to the NTW invariant $\nu$

[TABLE]

which is the generalization of the celebrated invariant of Thouless, Kohmoto, Nightingale and den Nijs (TKNN) for free fermions [32, 2] to the case of interactions and disorder. In Ref.[17] the TKNN invariant was generalized to Gaussian mixed states of fermions using the King-Smith – Vanderbildt relation [33] between changes of $\phi_{\textrm{MB}}$ to those of the many-body polarization P=\frac{1}{2\pi}\textrm{Im}\log\bigl{\langle}\Psi_{0}\bigr{|}\hat{U}\bigl{|}\Psi_{0}\bigr{\rangle} introduced by Resta [22], where $\hat{U}\equiv\hat{U}(2\pi)$: $\partial_{\lambda}\phi_{\textrm{MB}}=2\pi\partial_{\lambda}P$. This then allowed to replace the ground-state average by the trace over a density matrix \phi_{\textrm{EGP}}=2\pi P=\textrm{Im}\log\textrm{Tr}\bigl{\{}\rho\,\hat{U}\bigr{\}} defining the so-called ensemble geometric phase in [17].

Ensemble Geometric Phase: fractional case.–

In the absence of interactions, gapped ground states (of fermions) occur only at integer fillings per unit cell. This changes with interactions. Here gapped ground states can exist which have fractional fillings and the Lieb-Schulz-Mattis theorem and its generalizations [34, 35] tell us that they attain a ”topological order” accompanied by fractionalization and degeneracy [36]. In such a case the above arguments do not hold as a parameter loop of $\lambda$ does not return the initial state to itself (up to a phase), but in general to an orthogonal state in the ground-state manifold. For a $d$-fold degenerate subspace the NTW invariant (2) must instead be replaced by the gauge-invariant determinant of a Wilson loop [37]

[TABLE]

Here {\sf A}_{\mu\nu}(\theta)=i\bigl{\langle}\Psi_{0}^{\mu}\bigr{|}\partial_{\theta}\Psi_{0}^{\nu}\bigr{\rangle} is a $d\times d$ matrix, and ${\cal P}$ denotes path ordering.

We now argue that also $\nu_{\textrm{tot}}$ can be related to an expectation value of a unitary operator, which then allows for a generalization to mixed states. To see this, we note that following Niu, Thouless and Wu [30], $\nu_{\textrm{tot}}$ can be expressed as integral of the Berry curvature corresponding to any one of the degenerate ground states $|\Psi_{0}^{\mu}\rangle$ over an enlarged torus, extending either the time integration to ${\cal T}\cdot d$ or the $\theta$ integration to $2\pi\cdot d$, ($\mu=0,\dots,d-1$)

[TABLE]

As a consequence particle transport is integer quantized only after $d$ cycles of a Thouless pump [32].

As shown by Aligia and Ortiz [38, 39, 40], eq.(4) gives also the winding number of a many-body Berry phase

[TABLE]

which does not depend on the particular ground state. To see this, we note that the lattice Hamiltonian is invariant under spatial translation by one unit cell, described by the unitary lattice shift operator $\hat{T}$. For a $d$-fold ground-state degeneracy one can construct a basis set of ground states $\{|\Phi_{0}^{\mu}\rangle\}$ with $|\Phi_{0}^{\mu}\rangle=\hat{T}^{\mu}|\Phi_{0}^{0}\rangle$ and $\mu=0,\dots,d-1$. Since $\hat{T}^{-1}\hat{U}^{d}\hat{T}=\hat{U}^{d}e^{2\pi idN/L}$ with $N$ being the total number of fermions, and $dN/L$ is an integer for a fractional filling $1/d$, the Berry phases $\phi^{(d)}_{\mu}$ of all states in this basis have the same value. It should be noted that the finiteness of the absolute value $|\langle\Psi|(\hat{U})^{d}|\Psi\rangle|$ in the thermodynamic limit is an indicator of a localized, i.e. insulating ground state with filling $1/d$ per unit cell [41, 38].

Eq.(5) then allows to define a generalized ensemble geometric phase for mixed states

[TABLE]

which we now use to construct a topological winding number for mixed states.

Extended super-lattice Bose-Hubbard model. –

Let us discuss a specific example with interaction-induced fractional topological charges and associated fractional winding number [42, 43] , the one-dimensional extended super-lattice Bose Hubbard model (Ext-SLBHM). As shown in Fig.1(b), bosons at lattice site $j$, described by annihilation and creation operators $\hat{a}_{j}$ and $\hat{a}^{\dagger}_{j}$, and with onsite interaction strength $U$ move along a 1D lattice with alternating hopping $t_{1}$ and $t_{2}$ and a staggered potential $\Delta$. In addition there is a nearest-neighbor (NN) and next-nearest neighbor (NNN) interaction $V_{1}$ and $V_{2}$ respectively. The Hamiltonian reads in second quantization

[TABLE]

Here $\hat{n}_{j}=\hat{a}^{\dagger}_{j}\hat{a}_{j}$. The system has a unit cell of two sites. With periodic boundary conditions and in the absence of interactions one finds two single-particle energy bands with a finite gap, which closes only for $\Delta=t_{1}-t_{2}=0$. For sufficiently strong interactions there are further gap openings and Mott insulating (MI) states with fractional fillings emerge. The ground-state phase diagram, obtained from DMRG simulations [44] is shown in Fig.1(a), where MI phases with fractional fillings are indicated.

In the following we are interested in the phase with average filling of $\rho=1/4$ per lattice site, or $1/2$ per unit cell. Here the ground state is doubly degenerate for periodic boundary conditions if the number $L$ of unit cells is even. If one starts in one of the two many-body ground states, say $|\Psi_{1}\rangle$, and changes $t_{1}-t_{2}$ and $\Delta$ in a closed loop in parameter space such that the many-body gap remains open at all points, the ground-state manifold $\{|\Psi_{1}\rangle,|\Psi_{2}\rangle\}$ returns to itself up to a $U(2)$ rotation. Then $2$ loops in parameter space need to be performed for the initial state to return to itself modulo a phase. Similarly the average current $\langle\hat{j}\rangle$ needs to be integrated over $2$ periods of length ${\cal T}$ to lead to an integer-quantized number $\Delta n$ of pumped particles, which is verified by our numerical simulations in Fig.2a. This does no longer hold true, however, for finite temperatures. As expected and shown in the same figure the number of transported particles deviates substantially from unity as soon as the temperature approaches the many-body gap, since higher energy states are occupied.

Remarkably and in sharp contrast, the winding of the generalized EGP $\phi^{(2)}_{\textrm{EGP}}$ remains strictly unity even at temperatures on the order of the many-body gap as can be seen from Fig.2b. (It should be noted that increasing the temperature becomes numerically more demanding.)

It was shown in Ref.[17] for non-interacting fermions, that the winding of the EGP $\phi_{\textrm{EGP}}$ remains the same for all temperatures $T<\infty$. In the following we will give some arguments that also for interacting systems the EGP winding of a thermal state is identical to that of the ground state below some critical temperature which is larger than the many-body gap.

Finite-temperature winding. –

We first show that there exists a critical temperature $T_{c}$ different from zero, below which the winding of the EGP coincides with that of the many-body ground state $\Delta\phi_{\textrm{EGP}}|_{T<T_{c}}=\Delta\phi_{\textrm{EGP}}|_{T=0}$. Thus, different e.g. from the prediction in [16], a temperature-induced topological transition can only occur at a finite, nonzero temperature.

To this end we note that the change of $\phi_{\textrm{EGP}}$

[TABLE]

with $z(\lambda)\neq 0$ for $\lambda\in\delta{\cal C}$, is just the winding number or degree of a smooth map $z(\lambda)=$Tr$\{\rho(\lambda)\,(\hat{U})^{d}\}$ from a closed loop $\delta{\cal C}$ in parameter space $(\lambda)$ into the complex plane $\mathbb{C}$ of the polarization amplitude $z(\lambda)$. It measures the algebraic change of phase of $z$ as the variable $\lambda$ goes around the loop once. The degree of $z$ is by definition the number of solutions of $z=0$ inside $\delta{\cal C}$ taking into account their algebraic multiplicity (for more details about degree theory see [51]). Suppose for all temperatures $T$ with $T_{1}<T<T_{2}$, $z\neq 0$ on $\delta{\cal C}$, then $\Delta\phi_{\textrm{EGP}}$ is independent of temperature in the interval $(T_{1},T_{2})$. These arguments are in parallel with Hopf’s homotopy theorem [51]. Now, the condition $z(T,\lambda)\neq 0$ everywhere on $\delta{\cal C}$ can always be fulfilled in a finite range of temperatures starting at $T=0$ until for a critical temperature $T_{c}$

[TABLE]

for some $\lambda$ on $\delta{\cal C}$. This defines a temperature-induced topological phase transition.

Critical temperature. –

As shown in [17], $T_{c}=\infty$ for non-interacting fermions. To investigate the possibility of a finite-temperature topological phase transition in an interacting system, we have numerically calculated $|z(T)|$ for the extended SLBHM. The results are shown in Fig.3 in dependence of $T$ for different system sizes. One recognizes that $|z|$ remains approximately unity below temperatures that are a sizable fraction of the many-body gap $\Delta_{\textrm{gap}}$ followed by a fall-off, which becomes more pronounced with increasing system size. This behavior is very similar to the non-interacting case, for which one finds a double-exponential scaling (see supplementary material) |z(T)|\sim\exp\bigl{(}-2L\exp(-\beta\Delta_{\textrm{gap}})\bigr{)}. Note, that although $z(T)\to 0$ for $L\to\infty$, the temperature $T_{0}$, where $|z|$ starts to deviate from unity only scales logarithmically with system size $L$.

[TABLE]

From our numerical data it remains unclear if $|z(T)|$ has a strict zero, indicating a topological phase transition, or not and it would be interesting to investigate other interacting superlattice models such as [49, 50].

In order to show that interactions may lead to a finite value of $T_{c}$ let us consider the flattened Hamiltonian

[TABLE]

with $\Delta_{\textrm{gap}}=E_{1}-E_{0}$. Here all excited states have the maximum weight compatible with temperature and gap. Since the total Chern number of all excited bands must be opposite to that of the ground state, flattening a many-body Hamiltonian to such a form is expected to lead to the ”worst” case. Then for the polarization amplitude $z(T)=\textrm{Tr}\{\exp(-\beta\tilde{H})(\hat{U})^{d}\}/Z$ holds $\textrm{Tr}\{e^{-\beta\tilde{H}}(\hat{U})^{d}\}=z(0)d\left(e^{-\beta E_{0}}-e^{-\beta E_{1}}\right)+\textrm{Tr}\{(\hat{U})^{d}\}e^{-\beta E_{1}}$, where $z(0)=\frac{1}{d}\sum_{\mu}\langle\Psi_{0}^{\mu}|(\hat{T})^{d}|\Psi_{0}^{\mu}\rangle$ is the zero-temperature polarization amplitude. As shown in [41, 38] $|z(0)|\to 1$ in the thermodynamic limit $L\to\infty$ if the ground state is insulating. Since $\textrm{Tr}\{(\hat{U})^{d}\}$ does not depend on the Hamiltonian, its phase is fixed and does not change upon parameter variations. Thus the winding of $\phi_{\textrm{EGP}}^{(d)}$ for the flattened Hamiltonian $\tilde{H}$ at temperature $T$ remains equal to that in the ground state as long as $e^{-\beta E_{0}}-e^{-\beta E_{1}}>|\textrm{Tr}\{(\hat{U})^{d}\}|e^{-\beta E_{1}}/d$. With this we find for the critical temperature $T_{c}$, defined through eq.(9), (with $k_{B}=1$)

[TABLE]

As shown in the supplementary material $|\textrm{Tr}\{(\hat{U})^{d}\}|$ is intensive and bounded by a value of order ${\cal O}(d)$.

Summary. –

We have introduced a many-body topological invariant for finite-temperature states of interacting many-body systems with fractional filling and ground-state degeneracy. The invariant is based on a generalization of the EGP introduced in [17]. In the limit $T\to 0$ it coincides with the well-known NTW invariant [30]. We showed that there exists a critical temperature $T_{c}$, defined by a vanishing polarization amplitude, which generically is larger than the many-body gap. Below $T_{c}$ the EGP-based topological winding number is identical to the ground-state invariant. The generalized EGP also allows to detect topological properties from finite-temperature measurements. It can be measured either directly (see [17]) or obtained from the full counting statistics, measurable e.g. in ultra-cold atom experiments with a gas microscope. A non-trivial topological winding of the EGP can furthermore induce a quantized transport in a weakly coupled auxiliary system, prepared at a low temperature [23]. We illustrated our results for the Ext-SLBHM model at quarter filling, which has a two-fold degenerate ground state and an associated fractional topological charge of $1/2$. The arguments given here can be extended to interacting two-dimensional systems of the Chern class with translational invariance, and thus to systems like fractional Chern insulators with intrinsic topological order. These systems can be mapped to independent one-dimensional systems by transforming to momentum space in one of the spatial directions.

acknowledgement

Financial support from the DFG through SFB TR 185, project number 277625399 is gratefully acknowledged. The authors thank Sebastian Diehl and Matthias Moos for stimulating discussions and Richard Jen for help in the DMRG calculation of the ground-state phase diagram. We also thank J. Sirker for support with the LCRG code. M.F. thanks the KITP for hospitality, where part of this work was done and supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.

I Appendix

II Numerical methods

For computations beyond exact diagonalization routines we used the light cone renormalization group (LCRG) [47], which allows to calculate observables directly in the thermodynamic limit. To this end the LCRG restricts the time-evolution of an operator to the Hilbert space of an effective light cone using the Lieb-Robinson bound for a Hamiltonian with short range interactions [52], see Fig. 4. Besides applying the Lieb-Robinson bounds, the Hilbert space is truncated similarly to time-adaptive density matrix renormalization group (t-)DMRG algorithms, in its traditional formulation by S.R. White [53, 54]. The formalism can however also be understood in the language of matrix product states (MPS) [44]. To time-evolve the system stepwise the LCRG algorithm uses transfer matrices to enlarge the system size and a second order Trotter-Suzuki decomposition [55, 56]. Connecting to conventional MPS language, a tensor network can be easily transformed to the light cone shape as adjacent plaquettes would just trace out, see Fig 4. A description of the full algorithm can be found in Ref. [47] as well as a working code in the supplementary material of Ref. [47]. The algorithm has been generalized to open systems in Refs. [48, 57]. In particular temperature effects are included by extending the Hilbert space by auxiliary degrees of freedom and subsequent purification. This is explained in detail in Refs. [46, 58].

In this supplementary we will not go into further details of the LCRG time-evolution, for which we refer to [47, 48, 57], but will focus on the measurement of the momentum shift operator $\hat{U}$ from the main text. We will follow the notation and graphical representation of [47, 48, 57].

The LCRG has three major components. The two transfer matrices for enlarging the system to the left L and right R and the actual state in a light cone shape named C. All components have 4 indices, which are explained in the following. From a numerical viewpoint, the components are therefore refered to as tensor objects of the LCRG. The time-evolution is done by Trotter-Suzuki plaquettes applied to the initial tensor objects [47, 48], see fig. 5, (red and blue plaquettes).

The standard tensor then has two local bonds, denoted by $\sigma_{L}$, $\sigma_{R}$ in the index and two superblocks, denoted by $m_{L}$, $m_{R}$ in the index. Blue dotted lines in Fig.5 imply the summation over pairwise indices. In the numerics of the main text one bond contains two sites, as well as two additional sites for the auxilliary system needed for the purification ansatz. Thus a bond has dimension $d=2^{4}$, since we are only considering hardcore bosons in the numerics. We keep around 300 states for the superblocks $m_{L}$ and $m_{R}$, where we have checked convergence and made sure not to exceed a discarded weigth of $10^{-10}$ for the first period of the Thouless pump and $10^{-8}$ for the second. The momentum shift operator, needed to calculate the EGP of the main text, reads

[TABLE]

where $\hat{n}_{js}$ denotes the particle number operator of the $s$th site $(s\in\{1,2,\dots,n\})$ in the $j$th unit cell. $0\leq r_{s}\leq 1$ characterizes the position within the unit cell. Here we set the lattice constant to 1. When using the purification ansatz to account for finite temperature, we need to measure $\hat{n}_{js}\otimes\mathbb{I}$, where $\mathbb{I}$ is the identity matrix applied to the auxiliary system effectively tracing out the ancillary sites. Since the operator $\hat{n}_{js}$ is diagonal in real space, we can rewrite $\hat{U}$ as a product

[TABLE]

In the LCRG we can compute $\hat{U}$ e.g. by multiplying left transfer matrices to extend the local density matrix.

[TABLE]

[TABLE]

Applying eq. (16) twice, changing $U_{0}$ to $U_{1},U_{2}$ etc. and then tracing out the most left superblock, will leave us with the tensor object presented in fig. 5. By applying eq. (16) more than twice we can receive any wanted block size.

III Flattened Hamiltonian

Here we proof that $|\textrm{Tr}\{(\hat{U})^{d}\}|\sim{\cal O}(d)$, as argued in the main text. For this we consider a system of $N$ particles in a lattice with $L$ unit cells and a fractional filling of $1/d$ per unit cell, i.e. $N=L/d$.

[TABLE]

The constraint of constant particle number can most conveniently be written in terms of an integral of an angle $\theta$:

[TABLE]

where $\phi_{s}=2\pi r_{s}/N$. In the case of fermions the sums over occupation numbers $n$ extends over [math] and $1$ and we find

[TABLE]

By introducing $z=-e^{-i\left(\theta-\phi_{s}\right)}$ the product incide the bracket can be simplified to

[TABLE]

[TABLE]

We have used the fact that $e^{-2\pi ij/N}$ are complex roots of $1-z^{N}=0$. Thus we find

[TABLE]

The only non-vanishing contribution to the $\theta$ integration are terms in the product that contain $e^{-i\theta N}$ to the first power. Thus we find

[TABLE]

We note that for special cases the sum over the sites within a unit cell can vanish, but generically $|\textrm{Tr}\{(\hat{U})^{d}\}|$ is on the oder of $d$. A similar calculation can be made for bosons where the sums over occupation numbers runs from [math] to $N$.

IV Scaling of $|z(T,L)|$ for non-interacting fermions

Here we discuss the scaling of $|z(T,L)|$ for the case of a two-band model of non-interacting fermions at half filling. ($N=L,$ where $N$ is the total number of fermions and $L$ is the number of unit cells in the system) As shown in Ref.[17] of the main paper, the many-body polarization can be represented in the following form

[TABLE]

where

[TABLE]

and $\mathcal{Z}$ is the parition function

[TABLE]

The Hamiltonian $H$ can be represented as a direct sum of the the momentum-space Hamiltonian matrix $h_{k}=\overrightarrow{\alpha}_{k}\cdot\overrightarrow{\sigma}$, where $\overrightarrow{\sigma}$ is the vector of Pauli matrices and $\overrightarrow{\alpha}_{k}$ are real vectors.

Our ultimate goal is to investigate the amplitude of $z$ in the thermodynamic limit i.e. when $N=\frac{M}{2}>>1$. Following Ref.[17] of the main paper, expression (19) can be represented as

[TABLE]

where the unitary matrix $u_{k}$ transforms the Bloch Hamiltonian $h_{k}$ into a diagonal matrix i.e.$\overrightarrow{\alpha}_{k}\overrightarrow{\sigma}=\left|\overrightarrow{\alpha}_{k}\right|u_{k}^{\dagger}\sigma_{z}u_{k}$. The trace in eq.(21) can be calculated easily in the thermodynamic limit i.e. $N\rightarrow\infty$. Indeed,

[TABLE]

Using this asymptotic form the module of $\left\langle R\right\rangle$ takes the following simple form

[TABLE]

Then, using Eq. (20), one has

[TABLE]

We proceed by giving estimates for the lower and upper bounds of $\left|z\right|$.

Lower bound:

[TABLE]

where in the first line we have used the arithmetic-geometric inequality.

*Upper bound: *To derive an upper bound on $\left|z\right|$, we note first that the function $\ln\left(2e^{-x}+1+e^{-2x}\right)$ is a convex function and by Jensen’s inequality

[TABLE]

where $\Delta=\frac{1}{N}{\sum\limits_{k}}\left|\overrightarrow{\alpha}_{k}\right|$ is the spectral gap of the Hamiltonian. We thus have the following upper bound

[TABLE]

Therefore,

[TABLE]

We verify immediately that if

[TABLE]

the absolute value of $z$ never vanishes, i.e. $\left|z\right|>0$ even when $N\rightarrow\infty$. However, in general condition (28) is violated and there might be an exact zero of $|z|$.

Taking $\left|\overrightarrow{\alpha}_{k}\right|$outside the summation in the above inequality corresponds to assuming a flat single particle spectrum. In the low temperature limit, $\beta\Delta>>1$, this is well justified and it is easily seen that upper and lower bounds of $\left|z\right|$ are the same giving the general scaling

[TABLE]

Hence, despite of the fact that $\left|z\right|$ is exponentially small in system size, the relevant temperature, where $|z(T)|$ deviates from unity scales only logarithmically in system size i.e. $T\sim\Delta/\ln\left(2N\right)$.