Finite-temperature topological invariant for interacting systems
Razmik Unanyan Maximilian Kiefer-Emmanouilidis Michael Fleischhauer







Department of Physics and Research Center OPTIMAS, University of Kaiserslautern, 67663 Kaiserslautern, Germany

Department of Physics and Research Center OPTIMAS, University of Kaiserslautern, 67663 Kaiserslautern, Germany
Department of Physics, University of Manitoba, Winnipeg, Manitoba, Canada R3T 2N2

Department of Physics and Research Center OPTIMAS, University of Kaiserslautern, 67663 Kaiserslautern, Germany
Kavli Institute for Theoretical Physics, University of California, Santa Barbara, California 93106, USA






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 ν=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. 






Finite-temperature topological invariant for interacting systems
    Michael Fleischhauer
    December 30, 2023
================================================================








  
Introduction. –


Starting with the discovery of the quantum Hall effect <cit.>
topology has become an important paradigm for the understanding and classification of phases of matter <cit.>. 
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 <cit.>. For example, it was shown that the generalization of geometric
phases to density matrices based on the Uhlmann construction <cit.> leads to consistent topological
invariants in 1D <cit.>. Its application to higher dimensions <cit.> is however faced with
difficulties <cit.>. Other approaches predict an unphysical extensive number of topological phases for arbitrarily small temperatures
<cit.>.

Recently it was shown that the winding of the many-body polarization introduced by Resta <cit.> is an alternative
 topological invariant for Gaussian mixed states of fermions in 1D, termed ensemble geometric phase (EGP) <cit.>, which can also be
applied in 2D <cit.>.
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→∞ and thus these states have the same topological 
classification  as the corresponding ground states (following the Altland-Zirnbauer classification <cit.>). Despite being a genuine many-body quantity the EGP can be measured directly <cit.>. 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 <cit.>. For non-interacting bosons the EGP winding is always zero <cit.>.

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 <cit.> 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)  <cit.> 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(λ) that has a periodic dependence on some variable λ. 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⟩ are no longer a good eigen-basis, topological invariants must be defined
in terms of the many-body ground state |Φ_0⟩.  As suggested by Niu, Thouless and Wu <cit.>,  the Zak phase of Bloch wavefunctions
can be generalized by replacing the single-particle crystal momentum k 
by a twist angle θ 
of boundary conditions Φ(x_1,… x_j + L ,… x_N) = e^i θΦ(x_1,…  x_N). Since θ and θ+2π define the same boundary conditions, the parameter space (θ,λ) is a torus T^2. Twisted boundary conditions can be removed and replaced by periodic ones via a canonical transformation to a twisted Hamiltonian, H̅(θ) = Û(θ) H Û(θ)^-1 and |Ψ_0(θ)⟩ = Û(θ) |Φ_0(θ)⟩. 
Here


    Û(θ) = e^ i θX̂,  with  X̂ = 1/L∑_j=1^L∑_s=1^n (j + r_s) n̂_js


is the momentum shift operator with n̂_js denoting the particle number operator of the sth site (s∈{1,2,…,n}) in the jth unit cell, and the
lattice constant is a=1. 0≤ r_s ≤ 1 characterizes the position within the unit cell.
In terms of the |Ψ⟩'s the many-body equivalent of the Zak phase then reads
ϕ_MB = i ∫_0^2π dθ  ⟨Ψ_0(θ)|∂_θΨ_0(θ)⟩.


If |Ψ_0(θ)⟩ is a gapped and non-degenerate ground state, a slow (adiabatic) change of the parameter λ=λ(t) in a closed loop, such that the many-body gap does not close, induces a Thouless pump, described by a current density
⟨ĵ⟩  = ∂_θ E(θ) / ħ + i (⟨∂_t Ψ_0|∂_θΨ_0⟩ -⟨∂_θΨ_0|∂_t Ψ_0⟩).
Then following Niu, Thouless and Wu, averaging over twisted boundary conditions one finds a strictly integer-quantized particle transport 
over one time period T. The transported charge Δ n is then directly
related to the NTW invariant ν


    Δ n     =   1/4π∫_0^ T dt ∫_0^2π dθ    i [⟨∂_t Ψ_0|∂_θΨ_0⟩ -
    ⟨∂_θΨ_0|∂_t Ψ_0⟩]
       =   1/2π∫_0^ T dt  ∂ϕ_MB/∂ t=1/2π∮ dλ ∂ϕ_MB/∂λ=
    ν,


which is the generalization of the celebrated invariant of Thouless, Kohmoto,  Nightingale and den Nijs (TKNN) for free fermions <cit.> 
to the case of interactions and disorder.
In Ref.<cit.> the TKNN invariant was generalized to Gaussian mixed states of fermions using 
the King-Smith – Vanderbildt relation <cit.> between
changes of ϕ_MB to those of the many-body polarization 
P = 1/2πImlog⟨Ψ_0|Û|Ψ_0⟩
introduced by Resta <cit.>, where  Û≡Û(2π):
∂_λϕ_MB=2 π∂_λ P. 
This then allowed to replace the ground-state average  by the trace
over a density matrix
 ϕ_EGP = 2π P = ImlogTr{ρ Û} defining the so-called ensemble geometric phase in <cit.>.





  
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 <cit.>
tell us that they attain a "topological order" accompanied by fractionalization and degeneracy
<cit.>. 
In such a case the above arguments do not hold as a parameter loop
of λ 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 (<ref>) must instead be replaced by the gauge-invariant determinant of a
Wilson loop <cit.>


    ν_tot = 1/2π∫_0^ T dt ∂/∂ tImlog W(t),   W = P e^i∫_0^2π dθ  A(θ).


Here A_μν(θ) = i⟨Ψ_0^μ|∂_θΨ_0^ν⟩ is a d× d matrix,
and P denotes path ordering.  

We now argue that also ν_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 <cit.>, ν_tot can be expressed
as integral of the Berry curvature corresponding to any one of the degenerate ground states |Ψ_0^μ⟩ over an enlarged torus, extending either the
time integration to T· d or the θ integration to 2π· d, (μ=0,…,d-1)


    ν_tot = 1/2π∫_0^ T dt   ∫^2π d_0 dθ  Im ⟨∂_t Ψ_0^μ|∂_θΨ_0^μ⟩.


As a consequence particle transport 
is integer quantized only after d cycles of a Thouless pump <cit.>. 

As shown by Aligia and Ortiz <cit.>, eq.(<ref>) gives also the winding number of 
a many-body Berry phase


    ν_tot =∮dλ/2π∂ϕ^(d)_μ/∂λ,  withϕ^(d)_μ =  Imlog⟨Ψ_0^μ| (Û)^d |Ψ_0^μ⟩,


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 T̂. For a d-fold ground-state degeneracy one can construct a basis set of ground states {|Φ_0^μ⟩} with |Φ_0^μ⟩ = T̂^μ|Φ_0^0⟩ and μ =0,…, d-1. Since T̂^-1Û^dT̂ = Û^d e^2π i d N/L with N being
the total number of fermions, and d N/L is an integer for a fractional filling 1/d, the Berry phases ϕ^(d)_μ  of all states in this basis have the same value. 

It should be noted that the finiteness of the absolute value |⟨Ψ| (Û)^d|Ψ⟩| in the thermodynamic limit 
is an indicator of a localized, i.e. insulating ground state  with filling 1/d per unit cell <cit.>.

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


    ϕ_EGP^(d) =  ImlogTr{ρ  (Û)^d}


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 <cit.>
, the one-dimensional extended super-lattice Bose Hubbard model (Ext-SLBHM).  As shown in
Fig.<ref>(b), bosons at lattice site j, described by annihilation and creation operators _j and _j, and with onsite interaction strength U move along a 1D lattice with alternating hopping  t_1 and t_2 and a staggered potential Δ. 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


    H = -t_1 ∑_j,even_j _j+1 - t_2 ∑_j,odd_j _j+1 + h.a. 
           + Δ/2∑_j (-1)^j _j_j+ U/2∑_j 
     n̂_j(n̂_j-1) + ∑_j,d V_d  n̂_i n̂_i+d.


Here n̂_j = _j _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 Δ=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 <cit.> is
shown in Fig.<ref>(a), where MI phases with fractional fillings are indicated.





In the following we are interested in the phase with average
filling of ρ=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 |Ψ_1⟩, and changes  t_1-t_2 and  Δ in a closed loop in parameter space such that the many-body gap remains open at all points, the
ground-state manifold {|Ψ_1⟩, |Ψ_2⟩} 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 ⟨ĵ⟩ needs to be integrated over 2 periods of length T
to lead to an integer-quantized number Δ n of pumped particles, which is verified by our numerical simulations
in Fig.<ref>a. 
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 ϕ^(2)_EGP remains strictly unity even at temperatures on the order of the many-body gap as can be seen from Fig.<ref>b.  (It should be noted that increasing the temperature becomes numerically more demanding.)





It was shown in Ref.<cit.> for non-interacting fermions, that the winding of the EGP ϕ_EGP remains the same for all temperatures T<∞. 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 Δϕ_EGP|_T<T_c = Δϕ_EGP|_T=0.
Thus, different e.g. from the prediction in <cit.>, a temperature-induced topological transition can
only occur at a finite, nonzero temperature.

To this end we note that the change of ϕ_EGP


    Δϕ_EGP=2π z=Im∮1/z(  λ)  dz( λ)  /dλ
    dλ


with z(λ) 0 for λ∈δ C,
is just the winding number or degree of a smooth map z(λ)=Tr{ρ(λ) 
(Û)^d} from a closed loop δ C in parameter space (λ) 
into the complex plane ℂ of the polarization amplitude z(λ). It
measures the algebraic change of phase of z as the variable λ goes around
the loop once. The degree of z is by definition the number of
solutions of z=0 inside δ C taking into account their algebraic
multiplicity (for more details about degree theory see <cit.>). Suppose
for all temperatures T with T_1< T < T_2, z≠0 on δ C, then 
Δϕ_EGP is independent of temperature in the interval (T_1,T_2). 

These arguments are in parallel with Hopf's homotopy theorem
<cit.>. Now, the condition z(T,λ)≠0 everywhere on δ C 
can always be fulfilled in a finite range of temperatures starting at T=0 until for a critical temperature T_c 


    | z(T_c,λ)| =0,


for some λ on δ C. This defines a temperature-induced topological phase transition.







  
Critical temperature. –



As shown in <cit.>, T_c=∞ 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.<ref> 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  Δ_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)|∼exp(-2Lexp(-βΔ_gap)). Note, that although z(T)→ 0 for  L→∞, 
the temperature T_0, where | z| starts to deviate from unity only scales logarithmically 
with system size L.


    T_0  ∼ Δ_gap/ln L.
 

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 <cit.>.

In order to show that interactions may lead to a finite value of T_c let us consider the flattened Hamiltonian 


    H̃ = E_0 ∑_μ=1^d |Ψ_0^μ⟩⟨Ψ_0^μ| + E_1 ∑_j (0,μ)|Ψ_j⟩⟨Ψ_j|,


with  Δ_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)=Tr{exp(-βH̃) (Û)^d}/Z holds
Tr{e^-βH̃ (Û)^d} =  z(0) d (e^-β E_0 -e^-β E_1) 
+ Tr{ (Û)^d} e^-β E_1,
where z(0) =  1/d∑_μ⟨Ψ_0^μ| (T̂)^d|Ψ_0^μ⟩ is the zero-temperature polarization amplitude.
As shown in <cit.> | z(0)|→ 1 in the thermodynamic limit L→∞ if the ground state is insulating.
Since Tr{(Û)^d} does not depend on the Hamiltonian, its phase is fixed and does not change upon parameter variations. 
Thus the winding of ϕ_EGP^(d) for the flattened Hamiltonian H̃ at temperature T remains equal to that in the ground state as long as
e^-β E_0-e^-β E_1 > |Tr{(Û)^d}| e^-β E_1/d. With this we find for the critical temperature T_c, defined
through eq.(<ref>), (with k_B=1)


    T_c = Δ_gap/ln(1+|Tr{(Û)^d}| /d).


As shown in the supplementary material |Tr{(Û)^d}| is intensive and bounded by a value of order 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 <cit.>. 
In the limit T→ 0 it coincides with the well-known NTW invariant <cit.>. 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 <cit.>) 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 <cit.>.
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.




§ APPENDIX






§ NUMERICAL METHODS





For computations beyond exact diagonalization routines we used the light cone renormalization group (LCRG) <cit.>, 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 <cit.>, see Fig. <ref>. 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 <cit.>. The formalism can however also be understood  in the language of matrix product states (MPS) <cit.>. To  time-evolve the system stepwise the LCRG algorithm uses transfer matrices to enlarge the system size and a second order Trotter-Suzuki decomposition <cit.>. 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 <ref>. A description of the full algorithm can be found in Ref. <cit.> as well as a working code in the supplementary material of Ref. <cit.>. 
The algorithm has been generalized to open systems in Refs. <cit.>. 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. <cit.>.
 

In this supplementary we will not go into further  details of the LCRG time-evolution, for which we refer to <cit.>, but  will  focus on the measurement of the momentum shift operator Û from the main text. We will follow the notation and graphical representation of <cit.>.	







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 <cit.>, see fig. <ref>, (red and blue plaquettes).

The standard tensor then has two local bonds, denoted by σ_L, σ_R in the index and two superblocks, denoted by m_L, m_R in the index. Blue dotted lines in Fig.<ref> 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


    Û = e^ i 2πX̂,  with  X̂ = 1/L∑_j=1^L∑_s=1^n (j + r_s) n̂_js,


where n̂_js denotes the particle number operator of the sth site (s∈{1,2,…,n}) in the jth unit cell. 0≤ r_s ≤ 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 n̂_js⊗𝕀, where 𝕀 is the identity matrix applied to the auxiliary system effectively tracing out the ancillary sites. 
Since the operator n̂_js is diagonal in real space, we can rewrite Û as a product 


    Û = ∏ _j=1 ^ L Û _j = ∏ _j=1 ^ L e^ i 2πX̂_j,  with  X̂_j = 1/L∑_s=1^n (j + r_s) n̂_js.


In the LCRG we can compute Û e.g. by multiplying left transfer matrices to extend the local density matrix.


    ρ_local[m̅_L,σ̅_L,σ̅_R,m_L,σ_L,σ_R]  = ∑_m_R C^*[m̅_L,σ̅_L,σ̅_R,m_R]  C[m_L,σ_L,σ_R,m_R]



    (ρ U_0)_extend[m̅_L', σ̅_L”, σ̅_L', m_L', σ_L”, σ_L' ] = ∑_m̅_L,m_L,σ̅_L,σ̅_R, σ_R,σ_L    L^*[m̅_L', σ̅_L”, σ̅_L',m̅_L]  L[ m_L',  σ_L”,  σ_L', m_L] ×
       ×ρ_local[m̅_L,σ̅_L,σ̅_R,m_L,σ_L,σ_R]  U_0 [σ̅_L,σ̅_R, σ_L,σ_R].


Applying eq. (<ref>) 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. <ref>. By applying eq. (<ref>) more than twice we can receive any wanted block size. 









§ FLATTENED HAMILTONIAN



Here we proof that |Tr{(Û)^d}|∼ 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. 


    Tr{(Û)^d} = Tr{exp(2π i/L d ∑_j=1^L∑_s (j + r_s) n̂_js)}.


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


    Tr{(Û)^d} =  1/2π∫_0^2π dθ  e^iNθ∏_j,s∑_n 
    	e^-i n(θ -ϕ_s)exp{2π i/N n  j}


where ϕ_s = 2π r_s/N. In the case of fermions the sums over occupation numbers n extends over 0 and 1 and we find


    Tr{(Û)^d} =1/2π∫_0^2π dθ  e^iNθ∏_s[∏_j=1^L(1+e^-i(θ-ϕ_s) e^2π i j/N)]
            = 1/2π∫_0^2π dθ  e^iNθ∏_s[∏_j=1^N(1+e^-i(θ-ϕ_s) e^2π i j/N)]^d


By introducing z=-e^-i(  θ-ϕ_s) the product incide the bracket can be
simplified to

    ∏_j=1^N
    (  1+e^-i(  θ-ϕ_s)  e^2π ij/N)  =
    
    
    ∏_j=1^N
    (  1-ze^2π ij/N)  =


    ∏_j=1^N
    (  e^-2π ij/N-z)  =1-z^N.

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


    Tr{(Û)^d} =1/2π∫_0^2π dθ  e^iNθ∏_s[1-(-1)^N e^-iθ N e^iϕ_s N]^d


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


    |Tr{(Û)^d}| = d ∑_s e^2π i r_s =  O(d).


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




§ 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


    z=⟨ R⟩/𝒵,


where


    ⟨ R⟩ =𝚝𝚛[
    ∏_k=1^Nexp(  βα_kσ)
    ]  ,


and 𝒵 is the parition function


    𝒵=𝚝𝚛exp(  -β H)  =
    ∏_k=1^N(  1+e^β|α_k|)
    ∏_k=1^N(  1+e^-β|α_k|)  .


The Hamiltonian H can be represented as a direct sum of the the
momentum-space Hamiltonian matrix h_k=α_k·σ,
where σ is the vector of Pauli matrices and α_k
are real vectors.

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


    ⟨ R⟩ =𝚝𝚛[
    ∏_k=1^Nexp(  β|α_k|σ
    _z)  · u_k-1u_k^†]  ,


where the unitary matrix u_k transforms the Bloch Hamiltonian
h_k into a diagonal matrix i.e.α_kσ=|α_k| u_k^†σ_zu_k.
The trace in eq.(<ref>) can be calculated easily in the thermodynamic limit
i.e. N→∞. Indeed,


    u_k-1u_k^†≈1+2π/Nu_k-1∂
    	u_k^†/∂ k+...N>>1→
    1.


Using this asymptotic form the module of ⟨ R⟩ takes
the following simple form


    |⟨ R⟩|≈exp(
    	∑_k=1^Nβ|α_k|)  +exp(  -∑_k=1^Nβ|α_k|)  +𝒪
    	(  1/N)  .


Then, using Eq. (<ref>), one has


    | z|≈exp(
    	∑_kβ|α_k|)  +exp(  -
    	∑_kβ|α_k|)  /
    	∏_k(  2+e^β|α_k|+e^-β|α_k|)
    .


We proceed by giving estimates for the lower and upper bounds of
| z|.




  
Lower bound: 


    | z| =1+exp(  -2β∑_k=1|α_k|)  /
    	∏_k(  2e^-β|α_k|+1+e^-2β|α_k|)
    ≥1/∏_k(  2e^-β|α_k|+1+e^-2β|α_k|)
    ≥
          1/(  1+2
    		∑_k=1
    		e^-β|α_k|+
    		∑_k=1
    		e^-2β|α_k|/N)
    	^NN_F>>1→1/exp(  2
    	∑_k
    	e^-β|α_k|+
    	∑_k
    	e^-2β|α_k|)
    ,


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


Upper bound: To derive an upper bound on | z|,
we note first that the function ln(  2e^-x+1+e^-2x) is a
convex function and by Jensen's inequality


    ln(
    ∏_k(  2e^-β|α_k|+1+e^-2β|α_k|)
    )  =
          =
    ∑_kln(  2e^-β|α_k|+1+e^-2β|α_k|)
    ≥ln(  1+e^-βΔ)  ^2N,


where Δ=1/N∑_k|α_k|
is the spectral gap of the Hamiltonian. We thus have the following upper bound


    | z|≤1+exp(  -2NβΔ)
    /(  1+e^-βΔ)  ^2N.


Therefore, 


    exp(  -
    ∑_k
    e^-2β|α_k|-2∑_k
    e^-β|α_k|)
    ≤| z|≤1+exp(  -2NβΔ)
    /(  1+e^-βΔ)  ^2N.


We verify immediately that if


    ∑_k
    e^-β|α_k|<∞
    ,


the absolute value of z never vanishes, i.e. | z| >0 even when N→∞. However, in general 
condition (<ref>) is violated and there might be an exact zero of | z|. 



Taking |α_k|outside the
summation in the above inequality corresponds to assuming a flat single particle spectrum.
In the low temperature limit, βΔ>>1, this is well justified and it is easily seen that upper and lower
bounds of | z| are the same giving the general scaling


    | z|≈exp(  -2Ne^-βΔ)
    .


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



99


Klitzing-PRL-1980 K. V. Klitzing, G. Dorda, and M. Pepper, New Method for High-Accuracy Determination of the Fine-Structure Constant Based on Quantized Hall Resistance, Phys. Rev. Lett. 45, 494 (1980).

TKNN-PRL-1982 D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Quantized Hall Conductance in a Two-Dimensional Periodic Potential, Phys. Rev. Lett. 49, 405 (1982).

Tsui-PRL-1982 D. C. Tsui, II. L. Störmer, and A. C. Gossard, Two-Dimensional Magnetotransport in the Extreme Quantum Limit, Phys. Rev. Lett. 48, 1559 (1982)

Laughlin-PRL-1983 R. B. Laughlin, Anomalous Quantum Hall Effect: An Incompressible Quantum Fluid with Fractionally Charged Excitations,
Phys.Rev.Lett. 50, 1395 (1983).

Arovas-PRL-1984 D. Arovas, J. R. Schrieffer, and F. Wilczek, Fractional Statistics and the Quantum Hall Effect, Phys. Rev. Lett. 53, 722 (1984).

Xiao-RMP-2010 Di Xiao, Ming-Che Chang, and Qian Nu, 
Berry phase effects on electronic properties, Rev. Mod. Phys. 82, 1959 (2010).

Hazan-Kane-RMP-2010 M. Z. Hazan and C.L. Kane, Colloquium: Topological insulators, Rev.Mod. Phys. 82, 3045 (2010).

Wen-RMP-2017 Xiao-Gang Wen Colloquium: Zoo of quantum-topological phases of matter, Rev. Mod. Phys. 89, 041004 (2017).

Avron-NJP-2011 J. E. Avron, M. Fraas, G. M. Graf, and O. Kenneth, Quantum response of dephasing open systems, New J. Phys. 13 053042 (2011).

Bardyn-NJP-2013 C.-E. Bardyn, M. A. Baranov, C. V. Kraus, E. Rico, A. Imamoglu, P. Zoller, and S. Diehl, Topology by dissipation,
New J. Phys. 15, 085001 (2013).

Viyuela-PRL-2014 O. Viyuela, A. Rivas, and M. A. Martin-Delgado, Uhlmann
Phase as a Topological Measure for One-Dimensional
Fermion Systems, Phys. Rev. Lett. 112, 130401 (2014).

Huang-PRL-2014 Z. Huang and D. P. Arovas, Topological Indices for Open
and Thermal Systems via Uhlmann's Phase, Phys. Rev.
Lett. 113, 076407 (2014).

Viyuela-PRL-2014b O. Viyuela, A. Rivas, and M. A. Martin-Delgado, Two- Dimensional Density-Matrix Topological Fermionic Phases: Topological Uhlmann Numbers, Phys. Rev. Lett.
113, 076408 (2014).

Nieuwenburg-PRB-2014 E. P. L. van Nieuwenburg and S. D. Huber, Classification of
mixed-state topology in one dimension, Phys. Rev. B 90, 075141 (2014).

Linzner-PRB-2016 D. Linzner, L. Wawer, F. Grusdt, M. Fleischhauer, Reservoir-induced Thouless pumping and symmetry protected topological order in open quantum chains, Phys. Rev. B (R) 94, 201105 (2016).

Grusdt-PRB-2017 F. Grusdt Topological order of mixed states in quantum many-body systems Phys. Rev. B 95, 075106 (2017).

Bardyn-PRX-2018 C. E. Bardyn, L. Wawer, A. Altland, M. Fleischhauer, S. Diehl, Probing the topology of density matrices, Phys. Rev. X 8, 011035 (2018).

Altland2020 A. Altland, M. Fleischhauer, and S. Diehl, Symmetry classes of open fermionic quantum matter, arxiv:2007.10448

Bardyn-2017 C. E. Bardyn, A recipe for topological observables of density matrices, arxiv:1711.09735.


Uhlmann-Rep-Math-Phys-1986
A.Uhlmann, Parallel Transport and "Quantum Holonomy" along Density Operators, Rep. Math. Phys. 24, 229 (1986).

Budich-Diehl-PRB-2015 J.C. Budich and S. Diehl, Topology of density matrices, Phys. Rev. B 91, 165140 (2015).

Resta-PRL-1998 R. Resta Quantum Mechanical Position Operator in Extended Systems, Phys. Rev. Lett. 80, 1800 (1998).

Wawer L. Wawer, R. Li, and M. Fleischhauer, Quantized transport induced by topology transfer from interacting to 
non-interacting fermion chains, arxiv:2008.xyz



Altland-PRB-1997 A. Altland, and M. Zirnbauer Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures, Phys. Rev. B 55, 1142 (1997).

Schnyder-PRB-2008 Andreas P. Schnyder, Shinsei Ryu, Akira Furusaki, and Andreas W. W. Ludwig
Classification of topological insulators and superconductors in three spatial dimensions, Phys. Rev. B 78, 195125 (2008).

Ryu-NJPhys-2010 
Shinsei Ryu, Andreas P Schnyder, Akira Furusaki and Andreas W W Ludwig, Topological insulators and superconductors: ten-fold way and dimensional hierarchy, New J. of Phys. (2010).

Wawer-2 L. Wawer, R. Li, and M. Fleischhauer, Quantized transport induced by topology transfer between coupled one-dimensional lattice systems, preprint, arxiv:2009.04149

Mink-PRB-2019 C. Mink, M. Fleischhauer, R.G. Unanyan, Absence of topology in Gaussian mixed states of bosons
Phys. Rev. B 100, 014305 (2019).


Niu-JPhysA-1984 Q. Niu, and D.J. Thouless, Quantised adiabatic charge transport in the presence of substrate disorder and many-body interactions, J. Phys. A, 17, 2453 (1984).

Niu-PRB-1985 Q. Niu, D. J. Thoulsess and Y.-S. Wu, Quantized Hall conductance as a topological invariant, Phys. Rev. B 31, 3372 (1985).


ExtSLBHM F. J. Burnell, Meera M. Parish, N. R. Cooper, and S. L. Sondhi, Devil's staircases and supersolids in a one-dimensional dipolar bose gas, Phys. Rev. B, 80,174519 (2009).

Thouless-PRB-1983 D. J. Thouless, Quantization of particle transport, Phys. Rev. B 27, 6083 (1983).

King-Smith-PRB-1993 R. D. King-Smith and David Vanderbilt Theory of polarization of crystalline solids, Phys. Rev. B 47, 1651 (1993)

Lieb-AnnPhys-1961 E. Lieb, T. Schultz, and D. Mattis, Two soluble models of an antiferromagnetic chain,
Ann. Phys. (N.Y.) 16, 407 (1961).

Oshikawa-PRL-2000 Masaki Oshikawa
Commensurability, Excitation Gap, and Topology in Quantum Many-Particle Systems on a Periodic Lattice
Phys. Rev. Lett. 84, 1535 (2000).

Tao-PRB-1984 R. Tao and Yong-Shi Wu, Gauge invariance and fractional quantum Hall effect, Phys. Rev. B 30, 1097 (1984).

Wilczek-PRL-1984 F. Wilczek and A. Zee, Appearance of Gauge Structure in Simple Dynamical Systems, Phys. Rev. Lett. 52, 2111
(1984).


Aligia-PRL-1999 A.A. Aligia and G. Ortiz, Quantum Mechanical Position Operator and Localization in Extended Systems, Phys. Rev. Lett. 82, 2560 (1999) .

Aligia-EPL-1999 A. A. Aligia, Berry phases in superconducting transitions Europhys. Lett. 45, 411 (1999).

Ortiz-pss-2000 G. Ortiz, and A. A. Aligia, How Localized is an Extended Quantum System?, phys. stat. sol. (b) 220, 737 (2000).

Resta-Sorella-PRL-1999 R. Resta and S. Sorella, Electron Localization in the Insulating State, Phys. Rev. Lett. 82, 370 (1999).

Zheng-PRB-2016 Tian-Sheng Zeng, W. Zhu, and D. N. Sheng, Fractional charge pumping of interacting bosons in one-dimensional superlattice,
Phys. Rev. B 94, 235139  (2016).

Li-PRB-2017 R. Li and M. Fleischhauer, Finite-size corrections to quantized particle transport in topological charge pump, Phys. Rev. B 96, 085444 (2017).

DMRG see e.g. U. Schollwöck The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011).


supplement See appendix which includes Refs. <cit.> for more details on the
numerical methods, the flattened Hamiltonian and the scaling of | z(T,L)|.  

LCRG1 Adrian E. Feiguin and Steven R. White,
 Finite-temperature density matrix renormalization using an enlarged Hilbert space, Phys. Rev. B 72, 220401(R) (2005).
 
 LCRG2
Tilman Enss and Jesko Sirker Light cone renormalization and quantum quenches in one-dimensional Hubbard models, New J. Phys. 14, 023008 (2012). 

LCRG3 M. Kiefer-Emmanouilidis and J. Sirker Current reversals and metastable states in the infinite Bose-Hubbard chain with local particle loss, Phys. Rev. A 96, 063625 (2017).

Torio-PRB-2006 M. E. Torio, A. A. Aligia, G. I. Japaridze, and B. Normand, Quantum phase diagram of the generalized ionic Hubbard model for AB_n chains Phys. Rev. B 73, 115109 (2006).

Stenzel-PRA-2019 L. Stenzel, A. L. C. Hayward, C. Hubig, U. Schollwöck, and F. Heidrich-Meisner, Quantum phases and topological properties of interacting fermions in one-dimensional superlattices, Phys. Rev. A 99, 053614 (2019).

ReganDonal O'Regan, Yeol Je Cho, and Yu-Qing Chen, Topological
	Degree Theory and Applications, Taylor & Francis Group, LLC, 2006
 

lieb1972
E. H. Lieb  and D. W. Robinson, The finite group velocity of quantum spin systems, Comm. Math. Phys., 28 (3), 251 (1972).


white1992
S. R. White  Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett., 69, 2863 (1992).
white2004
S. R. White and A. E. Feiguin Real-Time Evolution Using the Density Matrix Renormalization Group, Phys. Rev. Lett., 93, 076401 (2004). 


Trotter
H. F. Trotter,
Proc. Amer. Math. Soc. 10,
545 (1959).

Suzuki1
M.Suzuki,
Commun. Math. Phys. 51, 183 (1976).

LCRG4 M. Kiefer-Emmanouilidis Renormalization Group Algorithms for Open
	One-Dimensional Quantum Systems in the
	Thermodynamic Limit, unpublished, Diploma Thesis, Technische Universität Kaiserslautern (2017).


LCRG5 F. Verstraete and J. J. García-Ripoll  and J. I. Cirac
Matrix Product Density Operators: Simulation of Finite-Temperature and Dissipative Systems, PhysRevLett.93.207204 (2004)





