Probing the topology of density matrices

TL;DR

This paper introduces the ensemble geometric phase (EGP), a new observable that preserves topological invariants in mixed quantum states, enabling topological detection even at finite temperature.

Contribution

It proposes the EGP as a measurable topological invariant for mixed states and demonstrates its robustness in one-dimensional Gaussian quantum systems.

Findings

EGP remains quantized in mixed Gaussian states

EGP detects spectral singularities in density matrices

Proposed interferometric measurement setup for ultracold atoms

Abstract

The mixedness of a quantum state is usually seen as an adversary to topological quantization of observables. For example, exact quantization of the charge transported in a so-called Thouless adiabatic pump is lifted at any finite temperature in symmetry-protected topological insulators. Here, we show that certain directly observable many-body correlators preserve the integrity of topological invariants for mixed Gaussian quantum states in one dimension. Our approach relies on the expectation value of the many-body momentum-translation operator, and leads to a physical observable --- the "ensemble geometric phase" (EGP) --- which represents a bona fide geometric phase for mixed quantum states, in the thermodynamic limit. In cyclic protocols, the EGP provides a topologically quantized observable which detects encircled spectral singularities ("purity-gap" closing points) of density matrices. While we identify the many-body nature of the EGP as a key ingredient, we propose a conceptually simple, interferometric setup to directly measure the latter in experiments with mesoscopic ensembles of ultracold atoms.

Figures

Equations

$${\phi }_{\text{E}}\equiv \mathrm{Im}\ln \mathrm{Tr}\left(\rho e^{i\delta k\hat{X}}\right),$$

$${\phi }_{\text{E}}={\phi }_{\text{Z}}+\Delta (N),\Delta (N)\stackrel{N\to \infty }{\to }0,$$

$$\Delta {\phi }_{\text{E}}=\oint d\varphi {\partial }_{\varphi }{\phi }_{\text{E}}=C,$$

$$P=\frac{1}{2\pi }\mathrm{Im}\ln ⟨{\psi }_0|\hat{T}|{\psi }_0⟩,\hat{T}\equiv e^{i\delta k\hat{X}},$$

$$P$$

$$\simeq \frac{1}{2\pi }\oint dk\mathrm{Tr}(A_k)\equiv \frac{{\phi }_{\text{Z}}}{2\pi },$$

$$\Delta P=\frac{i}{2\pi }\iint dkd\varphi \mathrm{Tr}\left(⟨{\partial }_{\varphi }{\psi }_0|{\partial }_k{\psi }_0⟩-⟨{\partial }_k{\psi }_0|{\partial }_{\varphi }{\psi }_0⟩\right),$$

$${\phi }_{\text{E}}=\mathrm{Im}\ln ⟨\hat{T}⟩,$$

$$\rho =\frac{1}{\mathcal{Z}}\exp \left(-\sum _{i,j}{\hat{a}}_i^{†}{G}_{ij}{\hat{a}}_j\right),$$

$$⟨{\hat{a}}_i^{†}{\hat{a}}_j⟩=[f(G){]}_{ij}.$$

$$⟨[{\hat{a}}_j,{\hat{a}}_i^{†}]⟩={\left[\tanh (\frac{G}{2})\right]}_{ij},$$

$$⟨\hat{O}({\hat{a}}^{†},\hat{a})⟩=\mathcal{N}\int d(\bar{\psi },\psi )e^{\bar{\psi }{f}^{-1}(G)\psi }\hat{O}(\bar{\psi },\psi ),$$

$$\hat{T}({\hat{a}}^{†},\hat{a})$$

$$=\prod _i\left(1-{\hat{a}}_i^{†}{\hat{a}}_i+{\hat{a}}_i^{†}{e}^{i\delta k{x}_i}{\hat{a}}_i\right)$$

$$=\prod _i{e}^{-{\hat{a}}_i^{†}(1-t_i){\hat{a}}_i},$$

$$⟨\hat{T}({\hat{a}}^{†},\hat{a})⟩\equiv ⟨\hat{T}⟩$$

$$=\det [-f(G)]\det [-f^{-1}(G)+1-T]$$

$$=\det [1-f(G)+f(G)T],$$

$$G=\sum _k{G}_k|k⟩⟨k|,$$

$$G_k=U_k{B}_k{U}_k^{†},B_k\equiv {\mathrm{diag}}_s({\beta }_{k,s}),$$

$$f(G)=\sum _k[U_{k}f(B_k)U_k^{†}]|k⟩⟨k|,$$

$$T=\sum _k|k+1⟩⟨k|,$$

$${\phi }_{\text{E}}$$

$$=\mathrm{Im}\ln \det [1+(1-f(G){)}^{-1}f(G)T],$$

$${\phi }_{\text{E}}$$

$$=\mathrm{Im}\ln \left(\exp \mathrm{Tr}\ln [1+(1-f(B){)}^{-1}f(B)U^{†}TU]\right)$$

$$=\mathrm{Im}\ln \left[\exp \mathrm{Tr}\ln (1+M_T)\right]$$

$$=\mathrm{Im}\ln \det (1+M_T),$$

$$M_T$$

$$=(-1{)}^{N+1}\prod _k{e}^{-B_k}{U}_{k+1,k},$$

$$(A_k{)}_{s{s}^{\prime }}=i⟨u_{k,s}|{\partial }_k{u}_{k,s^{\prime }}⟩.$$

$${\phi }_{\text{E}}\simeq \mathrm{Im}\ln \prod _k{e}^{-{\beta }_{k,0}}{e}^{i\delta k{A}_{k,0}},A_{k,0}\equiv i⟨u_{k,0}|{\partial }_k{u}_{k,0}⟩.$$

$${\phi }_{\text{E}}$$

$$\stackrel{N\to \infty }{\to }\oint dk{A}_{k,0}={\phi }_{\text{Z}},$$

$$U_{k+1,k}\simeq e^{i\delta k{A}_k}\equiv e^{i\delta k{\mathcal{A}}_k\cdot \sigma },$$

$$U_{k+1,k}$$

Full text

I Introduction

Topology has emerged as an important paradigm in the classification of ground states in many-particle quantum systems. Metaphorically speaking, topology enters the stage when the ground-state wavefunction of a complex quantum system contains “twists” (or “knots”), as a function of defining system parameters. Relevant parameters sets may include, e.g., the quasiparticle momenta labeling single-particle states in a translationally invariant system, the collective phase governing the macroscopic ground-state wavefunction of a superconductor, or the parameters controlling an external drive or pump.

Where topology is present, it is characterized by integer-valued invariants with a high degree of robustness with regard to perturbations, including parametric deformations of a Hamiltonian, translational symmetry breaking, or even the addition of particle interactions to symmetry-protected topological states in noninteracting systems. These invariants are generally formulated in terms of the zero-temperature ground state of topological quantum systems, which implicitly assumes that the system can be described by a single pure-state wavefunction. Realistic systems, however, are generally characterized by mixed states corresponding to thermal, or even more exotic distributions. Therefore, an obvious set of questions presents itself: to what extent can the concept of topology be generalized to finite-temperature states and, more broadly, to arbitrary mixed states (described by a density matrix instead of a state vector)? And, even more ambitiously, are possible formal generalizations connected to topologically quantized observables?

The fact that density matrices lend themselves to topological classification is reflected in the definition of various geometric phases and corresponding topological invariants. An important example is given by the Uhlmann phase Uhlmann (1986), a formal generalization of the geometric Berry phase Berry (1984); Simon (1983); Wilczek and Zee (1984). The definition of this phase is based on a gauge structure in the space of positive-definite density matrices Uhlmann (1986); Viyuela et al. (2014a, b); Huang and Arovas (2014); Budich and Diehl (2015). Recently, a different approach has been proposed Grusdt (2017) to generalize the concept of topological order in the sense of Ref. Chen et al. (2010) to mixed states, based on the equivalence of topologically identical states under local unitary transformations. While these approaches are formally elegant, they do not directly relate to observables that are readily representable in terms of system correlators 111Ref. Viyuela et al. (2016) proposes to use an ancillary system in topological insulators that are simulated by single qubits to implement the formal state purification required to access the Uhlmann phase..

In this work, we explore a distinct notion of parallel transport for mixed states within the comparatively simple class of symmetry-protected topological (SPT) insulators. Our approach is conceptually different from previous works in that the starting point of our construction is a many-body correlator which is directly observable, instead of the entire density matrix. One of its defining features is its reduction to common classifying topologies in the zero-temperature limit where thermal mixed states become projectors onto quantum ground states. Our construction can therefore be regarded as an extension of such notions to finite temperatures or, more broadly, to ensembles of mixed states. Most importantly, it shows that topological quantization can survive mixedness at the cost of dealing with genuine many-body correlators. Specifically, focusing on translation-invariant one-dimensional (1D) lattice systems of fermions in Gaussian states, we identify a physical many-body observable $\varphi_{\text{E}}$ which approaches, in the thermodynamic limit, a well-defined geometric phase for mixed states. Although realistic distributions will often be generated by thermalization, the concept applies to generic density matrices (subject to a few conditions spelled out below). For this reason, we dub $\varphi_{\text{E}}$ the “ensemble geometric phase” (EGP).

Before describing our key results, we start by defining the aforementioned notion of topological “twists” in the context of SPT ground states. There, distinct topological sectors are identified as distinct homotopy classes characterizing how single-particle states $\{|\psi_{\alpha}\rangle\}$ forming the many-particle ground state of a SPT insulator vary with parameter(s) $\alpha\equiv\{\alpha_{i}\}$. Depending on the context, the relevant parameters may be “internal”, such as the crystal-momentum components $\alpha_{i}\equiv k_{i}$ of a translationally invariant system, or “external”, such as the parameters $\alpha_{i}\equiv\phi_{i}$ of an imposed drive or pump protocol. The homotopy classes characterizing the map $\alpha\mapsto|\psi_{\alpha}\rangle$ can be described in terms of a $U(n)$ Berry connection or gauge field $(A_{i})_{ss^{\prime}}\equiv i\langle\psi_{\alpha,s}|\partial_{\alpha_{i}}\psi_{\alpha,s}\rangle$, where $n$ is the number of bands (indexed by $s$) composing the ground state. In this picture, topological invariants can be understood as Chern classes of this gauge field (or quantities that depend on the latter Ryu et al. (2010)). More importantly, invariants are often related to (thus quantized) measurable observables, such as Hall transport coefficients. Where a direct connection to an observable exists, high levels of stability, e.g., with regard to disorder or particle interactions, are to be expected.

All these concepts are beautifully exemplified in the Rice-Mele model — a paradigmatic model for noninteracting topological insulators in 1D Rice and Mele (1982) (with two bands and periodic boundary conditions; see Ref. Asboth et al. (2016) for a recent review). In this model, a two-dimensional toroidal parameter space $\alpha\equiv(k,\phi)$ is defined by the system 1D Bloch momentum $k$ and an external pump parameter $\phi$. The relevant topological invariant then describes how many times the parameter-dependent ground state of the system completely covers an effective Bloch sphere during a full cycle of the parameters. When the external parameter $\phi$ is periodically and adiabatically varied in time, this invariant describes the quantized charge pumped through the system. Formally, the invariant is related to the Berry connection via a Chern class. It can be understood as the winding in $\phi$ space of the geometric (Berry) phase defined by the loop integral of the connection over $k$ (i.e., over the Brillouin zone), which in the 1D context is commonly called a Zak phase Zak (1989). The Zak phase and the corresponding topological invariant are related to physical observables, i.e., to the zero-temperature (ground-state) polarization, and to the associated current flow.

Key results and outline

In this work, we consider the thermal equilibrium and nonequilibrium analogs of the afore-discussed noninteracting 1D systems, and extract topological information from their mixed states. The backbone of our construction is the expectation value of a many-particle momentum-translation operator Linzner et al. (2016), which was considered for pure states in a seminal work by Resta on the polarization of periodic systems Resta (1998). Specifically, for reasons that we motivate in Sec. II, we consider the ensemble geometric phase (EGP)

[TABLE]

where $\rho$ is the relevant density matrix, $\delta k\equiv 2\pi/L$ is the smallest possible momentum in a periodic system of size $L$, and $\hat{X}$ is the many-particle (center-of-mass) position operator.

The EGP is a natural generalization of the geometric Zak phase relevant to pure quantum states. In particular, we demonstrate below that the value of the EGP for a mixed quantum state is given by the zero-temperature Zak phase of the ground state up to corrections that vanish in the thermodynamic limit. We thus find a positive answer to the question whether geometric and topological properties of pure quantum states may retain their integrity at finite temperatures or in general mixed quantum states. In fact, our construction identifies an order parameter for topological phase transitions in general mixed quantum states. By studying concrete examples in and out of equilibrium, we demonstrate that topological phase transitions do exist in mixed quantum states, revealed by nontrivial (integer) changes in the winding of the EGP along a closed parameter cycle. We emphasize that single-particle quantities like the current, which can serve as topological order parameters at $T=0$, are no longer related to a geometric phase at any finite temperature, and thus cannot provide information on the existence of topological phase transitions for mixed states. We will give the EGP topological order parameter a concrete physical meaning by specifying an experimental detection protocol based on many-body interferometry.

In more detail, in Sec. III, we investigate the geometric nature of the EGP in the context of noninteracting systems, focusing on mixed states generically described by a Gaussian density matrix $\rho\sim e^{-G}$. (We also consider translationally invariant lattice systems, for convenience.) All information about fermionic Gaussian states is encoded in the Hermitian matrix $G$ or, equivalently, in the covariance matrix $\Gamma$ collecting all expectation values of bilinears of fermion creation and annihilation operators. It will be convenient to think of $G$ as a “fictitious Hamiltonian”. The spectrum of this Hamiltonian — the “purity spectrum” — describes the occupation probabilities of individual fermionic modes, and we will show that the existence of a gap in this spectrum — a “purity gap” — is required for the robustness of our construction [see Fig. 1(b)].

In Sec. III.1 and III.2, we show that the EGP defines a geometric phase for mixed states in the thermodynamic limit, in the sense that

[TABLE]

where $\varphi_{\text{Z}}$ is the Zak phase of the “ground state” — or lowest “purity band” — of the fictitious Hamiltonian $G$, and $\Delta(N)$ is a correction that vanishes in the limit of a large system size $N$ (number of unit cells) [see Fig. 1(a)]. More explicitly, $\varphi_{\text{Z}}$ is defined as the loop integral over the Brillouin zone of the $U(1)$ gauge field or Berry connection of the lowest purity band $|u_{k,0}\rangle$, i.e., $\varphi_{\text{Z}}=\oint dkA_{k}$, where $A_{k}\equiv i\langle u_{k,0}|\partial_{k}u_{k,0}\rangle$. Equation (2) establishes the emergence, in the thermodynamic limit, of a bona fide $U(1)$ geometric phase for mixed quantum states. The underlying mechanism is discussed in Sec. III.1 [and illustrated in Fig. 1(c),(d)]. The EGP is defined up to integer multiples of $2\pi$ (i.e., it is a “phase”), and its actual value on the unit circle is observable.

The above findings have an important consequence for EGP differences $\Delta\varphi_{\text{E}}$ accumulated over closed cycles in parameter space $(k,\phi)$ (where $\phi$ is an external parameter as above). Namely, we show that

[TABLE]

where $C$ is the Chern number associated with the lowest purity band $|u_{k,0}\rangle\equiv|u_{0}(k,\phi)\rangle$. This relation holds irrespective of the system size $N$, as the finite-size correction $\Delta(N)$ in Eq. (2) cannot contribute to the integer winding of $\partial_{\phi}\varphi_{\text{E}}$. This demonstrates the existence of an exactly quantized observable for mixed quantum states, with an explicit connection to microscopic parameters of the system.

For pure states in the usual zero-temperature ground-state scenario, the correction $\Delta(N)$ in Eq. (2) vanishes for any system size $N$, such that $\varphi_{\text{E}}=\varphi_{\text{Z}}$. Only in that case does the temporal variation $\partial_{\phi}\varphi_{\text{E}}$ of the EGP with respect to adiabatic changes of $\phi\equiv\phi(t)$ coincide with the physical charge current, and the difference $\Delta\varphi_{\text{E}}$ accumulated per adiabatic cycle correspond to the (quantized) transported charge. At finite temperature or in the nonequilibrium setting, in contrast, this connection breaks down. This can be understood from the fact that charge transport is related to the expectation value of a single-particle operator (the current), while $\Delta\varphi_{\text{E}}$ is related to a many-particle correlator [arbitrary powers of the single-particle operator $\hat{X}$ contribute to Eq. (1)]. This underlines the different phenomenology found in both cases: in charge transport, corrections to the current are intensive and finite in the thermodynamic limit, leading to the breakdown of exact quantization Wang et al. (2013); Nakajima et al. (2016). In contrast, corrections to the EGP of thermal or nonequilibrium mixed states vanish in the thermodynamic limit, allowing for the strict quantization of $\Delta\varphi_{\text{E}}$ (to nontrivial values).

In Sec. IV, we illustrate our general findings in minimal two-band models corresponding to the finite-temperature Rice-Mele model and to a nonequilibrium analog introduced in Ref. Linzner et al. (2016). In the thermal case, we demonstrate the existence of a topological phase transition (signaled by $\Delta\varphi_{\text{E}}$) at infinite temperature where the purity gap closes. In the other setting, we illustrate a chief feature of nonequilibrium dynamics, namely, the possibility that the purity gap closes at points where the “damping gap” of the Liouvillian describing the dynamics does not close Bardyn et al. (2013). Such singular points give rise to an observable nonzero $\Delta\varphi_{\text{E}}$ when encircled in parameter space, but are not associated with more conventional signatures such as divergent length and time scales in correlation functions.

Here, the requirement of adiabaticity in the conventional zero-temperature setting — the smallness, as compared to some spectral gap (the Hamiltonian gap), of the rate of parameter changes along some path — is replaced by a “purity adiabaticity” requirement: as detailed in Sec. III.5, the number of points at which the EGP must be measured or “sampled” along the relevant closed path in parameter space increases with decreasing purity gap. This reveals another analogy — and important fundamental difference — to the conventional zero-temperature setting: the adiabaticity condition comparing dynamical (energy or damping) scales is replaced, here, by a condition comparing dimensionless numbers.

The many-body nature of the correlator corresponding to the EGP [Eq. (1)] may shed doubts on its practical observability. In Sec. V, however, we propose a scheme harnessing the tools available in current experiments to measure this phase in mesoscopic ensembles of ultracold atoms, via Mach-Zehnder interferometry using photons.

At this point, we would like to mention some more related works: Refs. Avron et al. (2011, 2012) develop an adiabatic response theory for nonequilibrium (Liouvillian) dynamics for systems with few degrees of freedom, with recent generalization to a many-body context Albert et al. (2016), where the connection to Hamiltonian ground-state responses is elucidated. This construction yields geometric phases and quantized invariants only in cases where the (instantaneous) stationary state is pure, or mixed in a specific fine-tuned way, in contrast to our situation. Furthermore, observable geometric phases have been identified in pumping protocols for open quantum dots in the high-temperature regime Pluecker et al. (2017). This construction builds up on ideas for geometric phases in classical dissipative systems Landsberg (1992), which however do not relate to the geometric structure of the underlying quantum state.

II Resta polarization and its generalization

In this section, our goal is to construct the ensemble geometric phase $\varphi_{\text{E}}$, satisfying the following criteria:

1. (i)

The EGP is defined up to integer multiples of $2\pi$, i.e., it is a “phase”. 2. (ii)

Differential changes $\partial_{\phi}\varphi_{\text{E}}(\phi)$ with respect to a parameter $\phi$ are well defined and observable. 3. (iii)

As a direct consequence of (i), the normalized change $\frac{1}{2\pi}\Delta\varphi_{\text{E}}\equiv\frac{1}{2\pi}\oint d\phi\,\partial_{\phi}\varphi_{\text{E}}(\phi)$ accumulated over a closed loop in parameter space is integer quantized. 4. (iv)

The EGP has a simple enough physical meaning for differences $\Delta\varphi_{\text{E}}$ to be measurable in a realistic setup. 5. (v)

In limits where the relevant density matrix reduces to a projector onto a ground state, the EGP coincides with a conventional geometric phase (the Zak phase) characterizing ground-state band structures. A similar situation occurs when the EGP reduces to a projector onto an arbitrary pure state.

In the following, we construct a quantity satisfying these conditions, and show how it defines a geometric phase for mixed states. We then identify the topological nature of its quantization property, and define a notion of quantized adiabatic pump for mixed states.

II.1 Resta polarization

Key to our construction is a formulation of the electronic polarization of periodic quantum systems pioneered in an insightful paper by Resta Resta (1998). Resta argues that the textbook expression $P=\langle\hat{X}\rangle/L=\langle\psi_{0}|\hat{X}|\psi_{0}\rangle$ (setting $e=\hbar=1$) for the ground-state polarization of a 1D system with size $L$ in terms of the expectation value of the many-body position operator $\hat{X}\equiv\sum_{j}\hat{x}_{j}$ (where $\hat{x}_{j}$ is the position operator of individual particles $j$) is not applicable when the system is periodic as $\hat{X}$ is not a proper operator in the space of wavefunctions obeying periodic boundary conditions $\psi(L)=\psi(0)$. Instead, Resta suggests the alternative formula

[TABLE]

where $\delta k=2\pi/L$. In this form, the polarization is defined modulo an integer (as required from the periodicity of the host system), and is expressed in terms of a phase. Measurable incremental changes $\Delta P$ are more relevant than the value of $P$ itself. In particular, the introduction of an adiabatically slow time-dependent parameter $\phi(t)$ leads to the observable $\partial_{t}P(t)=I(t)$ corresponding to the electronic current. Importantly, the charge $Q\equiv\Delta P=\oint dt\,I(t)$ “pumped” during a cyclic protocol $t\in[0,\tau]$ with $\phi(\tau)=\phi(0)$ is integer quantized. This fact follows in full generality from the observation that, for ground states $\left\lvert\psi_{0}(t)\right\rangle$ that are weakly time dependent, the derivative $\partial_{t}P(t)$ is equal to an expression derived by Thouless and Niu Niu and Thouless (1984) for the current flowing in linear response to adiabatic changes. A more specific construction providing a connection to the topological band theory of noninteracting lattice systems considers $\left\lvert\psi_{0}\right\rangle$ as a Slater determinant formed from single-particle Bloch states $\left\lvert\psi_{k,s}\right\rangle$, where $k\in 2\pi\mathbb{Z}/N$ (for a system of length $L=Na$, with $N$ unit cells and lattice constant $a=1$), and $s=0,\ldots,n-1$ are band indices. Noticing that the many-body operator $\hat{T}$ introduced in Eq. (4) acts by translating all single-particle momenta by $\delta k$ (such that $k\to k-\delta k$), it is then straightforward to verify Resta (1998) that

[TABLE]

where $S_{k}$ is a $n_{\text{p}}\times n_{\text{p}}$ matrix (where $n_{\text{p}}$ is the number of particles in the system) formed by momentum-shifted ground-state wavefunction overlaps, namely, $(S_{k})_{s,s^{\prime}}\equiv\langle\psi_{k,s}|\hat{T}|\psi_{k,s^{\prime}}\rangle=\langle\psi_{k,s}|\psi_{k-\delta k,s^{\prime}}\rangle\simeq\delta_{s,s^{\prime}}+i\delta k(A_{k})_{s,s^{\prime}}$, with $(A_{k})_{s,s^{\prime}}\equiv i\langle\psi_{k,s}|\partial_{k}\psi_{k,s^{\prime}}\rangle$. The first expression in Eq. (5) identifies $P$ as the (normalized) phase of a “Wilson loop” corresponding to the product of overlap determinants $\det(S_{k})$ across the Brillouin zone. The second expression represents the same quantity as a discretized Berry phase, and the third provides a continuum approximation in terms of the Zak phase $\varphi_{\text{Z}}\equiv\oint dk\,\mathop{\mathrm{Tr}}(A_{k})$ (Ref. Zak (1989)), which corresponds to the loop integral of the multiband $U(n)$ Berry connection $A_{k}$ (where $n$ is the number of bands).

The introduction of a time-dependent parameter $\phi\equiv\phi(t)$ leads to variations $\Delta P=\int dt\,\partial_{t}P=\oint d\phi\,\partial_{\phi}P$. When performing a full adiabatic cycle, one finds

[TABLE]

where the expression in the integral is a trace over the Berry curvature corresponding to the Berry connection $A_{k}$ (and its analog $A_{\phi}$ for variations in $\phi$). Equation (6) measures the integer homotopy invariant of the map $(\phi,k)\to|\psi_{0}\rangle\equiv|\psi_{0}(k,\phi)\rangle$ from the torus defined by the two cyclic parameters $(\phi,k)$ to the ground-state manifold, i.e., $\Delta P$ corresponds to the number of times the ground-state wavefunction $|\psi_{0}(k,\phi)\rangle$ fully “covers” the torus in the process of a full parametric variation.

II.2 Generalization to mixed states: the EGP

Taking advantage of the above formulation of the Zak phase of a pure state, we now turn to mixed states and construct a generalization of Eq. (4) (see also Ref. Linzner et al. (2016)) designed to preserve all properties (i)–(v) above. We consider mixed states that arise as the unique stationary state of a gapped equilibrium or nonequilibrium quantum evolution — the direct analog of gapped nondegenerate pure states Bardyn et al. (2013). The corresponding density matrix $\rho$ can be decomposed in the generic form $\rho=\sum_{m}p_{m}\left\lvert\psi_{m}\right\rangle\left\langle\psi_{m}\right\lvert$, where $p_{m}>0$ is the probability of finding the system in state $\left\lvert\psi_{m}\right\rangle$. A natural generalization of the Zak phase $\varphi_{\text{Z}}$ [given by Eqs. (4) and (5)] would be the average phase $\bar{\varphi}_{\text{Z}}=\sum_{m}p_{m}\varphi_{\text{Z},m}$, where $\varphi_{\text{Z},m}$ is the Zak phase of each individual pure state $\left\lvert\psi_{m}\right\rangle$. This choice, however, trivially breaks property (i) above: the statistical average of phases defined modulo $2\pi$ is generally not defined modulo $2\pi$. Here, instead, we consider the phase of the statistical average $\sum_{m}p_{m}\langle\psi_{m}|\hat{T}|\psi_{m}\rangle$, i.e., we consider the “ensemble geometric phase”

[TABLE]

where $\langle\ldots\rangle\equiv\mathop{\mathrm{Tr}}(\rho\ldots)$.

Equation (7) is designed to satisfy the benchmark criteria (i)–(v) above. In particular, we note that $\varphi_{\text{E}}$ reduces to a Zak phase $\varphi_{\text{Z}}$ in the limit of pure states, which hints at the geometric nature of the EGP for mixed states. Differential changes $\partial_{\phi}\varphi_{\text{E}}(\phi)\sim\mathop{\mathrm{Im}}(\langle\hat{T}\rangle^{-1}\partial_{\phi}\langle\hat{T}\rangle)$ are physically observable, as we will discuss in Sec. V, and changes $\frac{1}{2\pi}\Delta\varphi_{\text{E}}\equiv\frac{1}{2\pi}\oint d\phi\,\partial_{\phi}\varphi_{\text{E}}(\phi)$ accumulated over parameter cycles are by construction integer quantized. For mixed states, $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ does not have the meaning of a pumped electric charge as in the case of pure states where it reduces to $\Delta P\equiv Q$ [Eqs. (5) and (6)]. Nevertheless, we will show that it can be measured in many-body interferometric protocols.

In the following, we will demonstrate for a wide class of mixed states — namely, Gaussian mixed states — that $\varphi_{\text{E}}$ is related to the Zak phase of a pure state up to corrections that vanish in the thermodynamic limit — thereby establishing the topological nature of the quantization of $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$. Gaussian mixed states can be represented by a quadratic Hermitian operator $\hat{G}$, and $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ corresponds, as we will show, to the ground-state topological invariant of the latter. In the specific case of thermal Gaussian states $\rho\propto e^{-\beta\hat{H}}$ — the finite-temperature extensions of zero-temperature ground states — the relevant operator is $\hat{G}=\beta\hat{H}$ (where $\beta$ is the inverse temperature), such that $\Delta\varphi_{\text{E}}$ reflects the topology of the ground state of $H$. Remarkably, for a gapped Hamiltonian $\hat{H}$, the quantity $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ remains quantized and coincides with the ground-state topological invariant $\Delta P$ of $\hat{H}$ at finite temperature, for as long as $\hat{G}=\beta\hat{H}$ is gapped (i.e., up to infinite temperature where $\beta=0$). In this regard, the observable $\Delta\varphi_{\text{E}}$ reflects the zero-temperature ground-state invariant in a more robust way than single-particle observables (such as linear-response conductances) whose quantization is affected by intensive (system-size independent) ratios of band gaps over temperature (see, e.g., Refs. Wang et al. (2013); Nakajima et al. (2016)).

III Geometric phase and topological invariant for mixed states

In this section, we explicitly relate the EGP defined in Eq. (7) to a geometric phase for mixed Gaussian states. To this end, we consider a set $\{\hat{a}^{\dagger}_{i},\hat{a}_{i}\}$ of fermionic creation and annihilation operators, where $i\equiv(r,s)$ is a composite index where $r=0,\ldots,N-1$ labels unit cells, and the band index $s=0,\ldots,n-1$ indexes fermionic sites in the unit cell (with lattice constant $a=1$). We consider Gaussian states defined by a density operator of the form

[TABLE]

with $G$ is a Hermitian matrix and $\mathcal{Z}$ is a normalization constant ensuring that $\mathop{\mathrm{Tr}}(\rho)=1$. The matrix $G$, which we call the state “fictitious Hamiltonian”, plays a key role in what follows: it uniquely identifies the state, and its spectrum defines what we call the state “purity spectrum”. The latter essentially corresponds to the spectrum of $-\ln\rho$, and its eigenvalues indicate the purity of the state in the corresponding eigenspaces Bardyn et al. (2013). Note that $G=\beta H$ for a thermal state $\rho\propto e^{-\beta\hat{H}}$ (where $H$ is the matrix representation of $\hat{H}$), as anticipated above.

The Gaussian density matrices that we focus on, defined in Eq. (8), describe states that are fully characterized by single-particle correlations of the form $\langle\hat{a}^{\dagger}_{i}\hat{a}_{j}\rangle$, with vanishing “anomalous” correlations $\langle\hat{a}^{\dagger}_{i}\hat{a}^{\dagger}_{j}\rangle$ or $\langle\hat{a}_{i}\hat{a}_{j}\rangle$. Such states are typically found in equilibrium or nonequilibrium systems of noninteracting fermions without particle-number fluctuations. It is straightforward to check that operator expectation values calculated with respect to the density matrix in Eq. (8) are given by

[TABLE]

where $f(G)=(e^{G}+\mathbb{1})^{-1}$, with vanishing anomalous correlations. The same information is often collected in the so-called “covariance matrix” of the distribution, with matrix elements

[TABLE]

Higher-order correlation functions can be calculated using Wick’s theorem. Alternatively, and more efficiently here, one may compute the operator correlation functions of Gaussian states via a Grassmann integral, namely,

[TABLE]

where $\hat{O}(\hat{a}^{\dagger},\hat{a})$ is a (normal-ordered) operator defined in terms of the creation and annihilation operators, $\bar{\psi}_{i}$ and $\psi_{i}$ are Grassmann variables, $\hat{O}(\bar{\psi},\psi)$ is obtained by the formal replacement $a_{i}^{\dagger}\to\bar{\psi}_{i}$, $a_{i}\to\psi_{i}$, and $\mathcal{N}=\mathrm{det}[-f(G)]$ normalizes the integral. We note that $\langle\hat{a}^{\dagger}_{i}\hat{a}_{j}\rangle\to\langle\bar{\psi}_{i}\psi_{j}\rangle=[f(G)]_{ij}$ readily follows from the Gaussian form of the integral, and the combinatorial signs in $\langle\hat{a}^{\dagger}_{i}\hat{a}^{\dagger}_{k}\hat{a}_{j}\hat{a}_{l}\rangle=-f(G)_{ij}f(G)_{kl}+f(G)_{il}f(G)_{kj}$ are faithfully reproduced by the Grassmann anticommutation $\psi_{i}\psi_{j}=-\psi_{j}\psi_{i}$.

We now use the above Grassmann representation to calculate the operator expectation value $\langle\hat{T}\rangle$ in Eq. (7). To this end, we first express the operator in the form

[TABLE]

where we have defined $t_{i}=\exp(i\delta kx_{i})$. The normal ordering of this expression is a straightforward operation as the indices $i$ carried by the factors entering the product $\prod_{i}$ are all different and, hence, nontrivial commutators do not appear. Next, we substitute $\hat{T}(\bar{\psi},\psi)=\exp[-\sum_{i}\bar{\psi}_{i}(1-t_{i})\psi_{i}]\equiv 1-\bar{\psi}(\mathbb{1}-T)\psi$, where $T\equiv\mathrm{diag}(t_{i})$, into the Grassmann integral to obtain

[TABLE]

which is still a formal expression at this point, since $f(G)$ is a highly nondiagonal matrix in position space. A more tangible representation can be obtained by assuming that the state is translation invariant 222This requires both the dynamics and the initial state to be translation invariant Rivas et al. (2013)., in which case the matrix $G$ can be prediagonalized in a “Bloch basis” as

[TABLE]

where $\langle x|k\rangle=N^{-1/2}e^{ixk}$. Here, $G_{k}$ is a $n\times n$ Hermitian matrix defined in band space, with elements $(G_{k})_{s,s^{\prime}}$, which we can be cast in diagonal form

[TABLE]

where $U_{k}$ is a unitary matrix collecting the “purity eigenstates” (analogous to Bloch vectors) of the fictitious Hamiltonian $G_{k}$ representing the state, and $B_{k}$ is a diagonal matrix containing the “purity eigenvalues” $\beta_{k,s}\in\mathbb{R}$ of the latter, ordered in increasing order $\beta_{k,1}\leq\ldots\leq\beta_{k,n}$, for convenience [see Fig. 1]. The “Fermi functions” appearing in Eq. (13) can then be expressed as

[TABLE]

where $f(B_{k})=(e^{B_{k}}+\mathbb{1})^{-1}$. The matrix $T$ in Eq. (13) (originating from the many-body translation operator $\hat{T}$) takes an intuitive form in Bloch basis: writing $x_{i}\equiv x_{r,s}\equiv x_{r}+x_{s}$, where $x_{s}$ is the position shift of the fermionic site $s$ located in the unit cell $r$ with position $x_{r}$, we obtain $\langle k|T|k^{\prime}\rangle=\delta_{k,k^{\prime}+1}$ 333More precisely, one finds $\langle k|T|k^{\prime}\rangle=\delta_{k,k^{\prime}+1}S$ with $S=\text{diag}_{s}(e^{i\delta kx_{s}})$. However, the matrix $S$ describing momentum shifts in the unit cell plays a very minor role here: it is diagonal in band space and can essentially be removed by a gauge transformation of the Bloch eigenstates of $G_{k}$, $\left\lvert u_{k,s}\right\rangle\to(S^{\dagger})^{k}\left\lvert u_{k,s}\right\rangle$. This transformation leads to a modified boundary condition $\left\lvert u_{N,s}\right\rangle=S^{N}\left\lvert u_{k,0}\right\rangle$, which is irrelevant for the geometric and topological properties of the EGP. We thus set $S=\mathbb{1}$, for simplicity.. [Note that we use $k$ interchangeably as a momentum index ($k=0,\ldots,N-1$) and as the momentum itself (with a factor $2\pi/N$).] Consequently, the matrix $T$ in Eq. (13) has the Bloch representation

[TABLE]

i.e., $T$ can be considered as a unit matrix on the diagonal next to the principal diagonal.

The EGP defined in Eq. (7) corresponds to the complex phase of $\langle\hat{T}\rangle$ in Eq. (13). Since we are only interested in this phase, we can write

[TABLE]

where we have observed that $\ln\det[\mathbb{1}-f(G)]=\ln\det[\mathbb{1}-f(B)]$ is real and thus and does not contribute to the imaginary part of $\varphi_{\text{E}}$. Note that this assumes that the matrix $\mathbb{1}-f(G)$ is invertible, which is generically true for mixed states 444In the limit of pure states, one can always choose a proper regularization of $\mathbb{1}-f(G)$ to make it invertible.. To reduce $\varphi_{\text{E}}$ to a more transparent form, we use Eqs. (16) and (17) along with the identity $\det(\ldots)=\exp\mathop{\mathrm{Tr}}\ln(\ldots)$ to write

[TABLE]

where the trace and determinant act in band space, in the last two lines, $U\equiv\{\delta_{k,k^{\prime}}U_{k}\}$, $B\equiv\{\delta_{k,k^{\prime}}B_{k}\}$, and $M_{T}$ is a path-ordered matrix product (“transfer matrix”)

[TABLE]

with “link matrices” defined as $U_{k+1,k}=U^{\dagger}_{k+1}U_{k}$. The matrices $U\equiv\{\delta_{k,k^{\prime}}U_{k}\}$ and $B\equiv\{\delta_{k,k^{\prime}}B_{k}\}$ appearing in Eq. (19) are block diagonal. In the crucial third equality in Eq. (19), we have noted that $U^{\dagger}TU=\sum_{k}U^{\dagger}_{k+1}U_{k}|k+1\rangle\langle k|$ is a matrix with blocks on the next-to-leading diagonal. Matrices of this form, viz. $A=\sum A_{k}|k+1\rangle\langle k|$, have the property that $A^{N}=(\prod_{k}A_{k})\,\mathbb{1}$, while powers of $A$ different from multiples of $N$ do not contain diagonal matrix elements, and, hence, do not contribute to the expansion of expressions of the form $\mathop{\mathrm{Tr}}\ln(\mathbb{1}+A)$. This feature was used to arrive at the final expression in Eq. (19) containing the transfer matrix $M_{T}$. For convenience, we will assume an odd number of sites throughout, such that the factor $(-1)^{N+1}$ in Eq. (20) can be omitted.

Equations (19) and (20) express the EGP as a path-ordered product of link matrices $U_{k+1,k}$ with points $k=0,\ldots,N-1$ along a closed loop corresponding to the Brillouin zone. This structure is reminiscent of discretized Wilson loop Fukui et al. (2005), where $U_{k+1,k}$ plays the role of a discrete $U(n)$ gauge connection. In the following, we will build on this observation to show how the EGP reduces to a more simple $U(1)$ gauge structure in the limit of large system sizes.

III.1 Gauge-reduction mechanism

Equations (19) and (20) directly relate the EGP to a path-ordered product of two types of matrices: the link matrices $U_{k+1,k}$ and the “weight factors” $e^{-B_{k}}=\text{diag}_{s}(e^{-\beta_{k,s}})$. The link matrices describe the “geometry” underlying the band structure of the mixed state $\rho$. This can be understood from the fact that the $n\times n$ matrices $U_{k}\equiv(\left\lvert u_{k,0}\right\rangle,\ldots,\left\lvert u_{k,n}\right\rangle)$ contain the $k$-dependent Bloch eigenvectors diagonalizing the fictitious Hamiltonian $G_{k}$ representing the state [see Eq. (8)]. In the limit of large $N$ where points $k$ and $k+1$ are infinitesimally close, the link matrices take the unitary form $(U_{k+1,k})_{ss^{\prime}}=\langle u_{k+1,s}|u_{k,s^{\prime}}\rangle\simeq 1-\delta k\left\langle u_{k,s}\right\lvert\partial_{k}u_{k,s^{\prime}}\rangle\simeq\exp[i\delta k(A_{k})_{ss^{\prime}}]$, where $A_{k}$ is the (non-Abelian) Berry connection

[TABLE]

Therefore, the link matrices $U_{k+1,k}$ describe the geometric “twist” of the state band structure (or “purity” bands) when moving from $k$ to $k+1$ in the Brillouin zone.

The weight factors $e^{-B_{k}}=\text{diag}_{s}(e^{-\beta_{k,s}})$, on the other hand, are purely real and determine the statistical weight with which a given purity band $s$ contributes to the EGP. It is at this point that the many-body nature of the correlator $\langle\hat{T}\rangle$ in Eq. (7) really kicks in: while conventional few-body expectation values take the symbolic form $\sim\mathop{\mathrm{Tr}}[e^{-G}(\ldots)]$, i.e., a structure where weight factors appear once, here we have a much more rigid structure of the form $\mathop{\mathrm{Tr}}[\ldots(\ldots)e^{-G_{k+1}}(\ldots)e^{-G_{k}}\ldots]$, in which selection through weight factors occurs a thermodynamically large number of times $N$ [see Fig. 1]. For mixed states with a purity gap, weight factors efficiently select the lowest band $s=0$ with weight $\beta_{k,0}$ (for thermal states $\sim e^{-\beta\hat{H}}$, we recall that $\beta_{k,s}=\beta\epsilon_{k,s}$, where $\epsilon_{k,s}$ is the energy spectrum of the underlying Hamiltonian $\hat{H}$). In that case, Eq. (20) leads to a crucial gauge reduction to a single $U(1)$ component:

[TABLE]

As anticipated above, the relevant Berry connection for the EGP is now the $U(1)$ Berry connection $A_{k,0}$ describing geometric properties of the lowest purity band.

As we will demonstrate below, corrections to Eq. (22) generally vanish in the thermodynamic limit $N\to\infty$. In addition, differences $\frac{1}{2\pi}\Delta\varphi_{\text{E}}\equiv\frac{1}{2\pi}\oint d\phi\,\partial_{\phi}\varphi_{\text{E}}(\phi)$ accumulated per parameter cycle are topologically quantized irrespective of the system size $N$. The underlying gauge-reduction mechanism applies, in particular, to the thermal density matrices $\rho\sim\exp(-\beta\hat{H})$ of topological band insulators: while the probing of topological invariants via conventional response coefficients generally leads to compromised results for temperatures $\beta\Delta\epsilon\gtrsim 1$ exceeding the gap $\Delta\epsilon$ of $\hat{H}$, the observable $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ considered here remains topologically quantized at finite temperature. The only requirement is that the purity gap remains finite, i.e., $\beta\Delta\epsilon>0$. In the next section, we will illustrate the reduction of the EGP to a geometric phase and the scaling of corrections in a simple example.

III.2 Two-band example

In the following, we discuss the EGP in the illustrative case of a two-band model with gapped purity spectrum $\pm\beta_{k}\not=0$. As a warmup, we first examine the case in which the density matrix reduces to a projector onto a pure state, which corresponds to the limit $\beta_{k}\equiv\beta\to\infty$. In this setting, the weight factors $e^{-(-\beta_{k})}=e^{\beta}\to\infty$ of the lower purity band dominate over those of the upper band (with weights $e^{-\beta}$), and Eq. (22) reduces to

[TABLE]

where all approximate equalities become exact in the limit $\beta\to\infty$. Note that the accumulated weight factors $\prod_{k}e^{\beta}=e^{\beta N}$ drop out as they do not contribute to the imaginary part. Therefore, in the pure-state limit, the EGP reduces to the Zak phase of the lower purity band, as expected 555Recall that, for thermal states $\rho\sim e^{-\beta H}$, the lower purity band corresponds to the ground state of $H$..

We now examine the case of mixed states, anticipating that Eq. (23) will be reproduced up to corrections that vanish in the thermodynamic limit. The EGP is given by the path-ordered product defined by Eqs. (19) and (20). We first parameterize the $2\times 2$ unitary link matrices:

[TABLE]

where $\mathcal{A}_{k}\equiv(\mathcal{A}_{k}^{0},\dots,\mathcal{A}_{k}^{3})^{T}$ is a real four-component vector, and $\mathcal{A}_{k}\cdot\sigma\equiv\sum_{i=0}^{3}\mathcal{A}_{k}^{i}\sigma_{i}$ is a shorthand for the expansion of the $2\times 2$ matrix $A_{k}$ in terms of Pauli matrices (with $\sigma_{0}\equiv\mathbb{1}$). In the two-band setting examined here, the diagonal matrix $e^{-B_{k}}\equiv\text{diag}_{s}(e^{-\beta_{k,s}})$ in Eq (20) takes the form $e^{-\beta_{k}\sigma_{3}}$. For convenience, we decompose the link matrices in the form

[TABLE]

where $U^{\text{diag}}_{k+1,k}\equiv e^{i\delta k(\mathcal{A}_{k}^{0}\sigma_{0}+\mathcal{A}_{k}^{3}\sigma_{3})}$ is diagonal in band space, and $V_{k}\equiv\delta k(\mathcal{A}_{k}^{1}\sigma_{1}+\mathcal{A}_{k}^{2}\sigma_{2})$ causes transitions between bands. Defining $W_{k_{1},k_{2}}=\prod_{k_{2}\leq k\leq k_{1}}U^{\text{diag}}_{k+1,k}e^{-\beta_{k}\sigma_{3}}$, we may then expand the “$\ln\det$” in Eq. (19) perturbatively in transition matrix elements as

[TABLE]

One may read these expressions in the spirit of time-dependent perturbation theory, where $k$ plays the role of time, and $W_{k+1,k}$ corresponds to the unperturbed discrete time-dependent propagator. The presence of weight matrices in the “Green’s function” $G$ implies that

[TABLE]

where we have noticed that $\mathcal{A}_{k}^{0}-\mathcal{A}_{k}^{3}=i\langle u_{k,0}|\partial_{k}u_{k,0}\rangle\equiv A_{k,0}$ is the $U(1)$ gauge field of the lowest purity band, and have defined the “average purity” of the lowest band $\bar{\beta}=N^{-1}\sum_{k}\beta_{k}$. The key point here is the exponential enhancement $\sim e^{\bar{\beta}N}$ with system size $N$ of the lower-band contribution relative to that of the upper band (with weight $\sim e^{-\bar{\beta}N}$). Inserting the above expression for $G^{-1}$ into Eq. (26) leads to a zeroth-order contribution which reproduces the result $\varphi_{\text{E}}=\varphi_{\text{Z}}$ obtained in Eq. (23). While the exponential factor $e^{\bar{\beta}N}$ drops out upon taking the imaginary part of the logarithm, it does play a role in the experimental detection of the EGP. We will return to this point in Sec. V.

Next, we examine perturbative corrections: in Eq. (26), contributions with an odd number of interband-transition matrices $V_{k}$ vanish, and we focus on the leading second-order contribution. Referring to Appendix B for details, we find that corrections to the zeroth-order result $\varphi_{\text{E}}=\varphi_{\text{Z}}$ scale as

[TABLE]

where $\Delta\beta=2\min_{k}{\beta_{k}}$ is the purity gap, and $c$ is a real constant independent of $N$ and $\Delta\beta$. Specifically, $\Delta(N)\sim(T/N)^{2}$ for thermal equilibrium states, implying that the EGP retains its zero-temperature value $\varphi_{\text{E}}=\varphi_{\text{Z}}$ even for finite $N$. Higher-order corrections vanish with higher powers of $N$ and, hence, are of negligible relevance. To understand the power-law suppression in an intuitive way, note that a thermodynamically nonvanishing (independent of $N$) correction can only arise if the smallness of the weight $\delta k\sim N^{-1}$ multiplying the $k$-local action of the transition operators $V_{k}$ gets compensated by an unconstrained summation over $k$ [since $\sum_{k}\delta k=\mathcal{O}(1)$]. However, excitations from $k$ to $k^{\prime}$ of the (symbolic) form $\sum_{k,k^{\prime}}(G_{0})_{1,k}V_{k}(G_{1})_{k,k^{\prime}}V_{k^{\prime}}G_{k^{\prime},N}$ from the lower band “[math]” with “propagator” $G_{0}$ to the excited band “$1$” with propagator $G_{1}$ get weighed by a factor $\sim\exp(-\Delta\beta|k-k^{\prime}|)$, on account of the spectral weights $e^{-\beta_{k}}$ of the density matrix penalizing excursions into the excited sector. This leads to “confinement” $|k-k^{\prime}|\sim 1/(\Delta\beta)$, and implies that $\delta k^{2}\sum_{k,k^{\prime}}e^{-\Delta\beta|k-k^{\prime}|}\sim 1/(\Delta\beta N)$. The fact that the correction actually scales as $1/(\Delta\beta N)^{2}$ with power two has to do with the fact that the leading-order perturbative expression comes out real, such that an additional factor $1/(\Delta\beta N)$ must be paid to obtain an imaginary contribution. The above mechanism applies regardless of the order of perturbation theory, and establishes the strong robustness of the geometric phase for mixed states defined by the EGP. Our perturbative calculations detailed in Appendix B are supported by numerical simulations for various equilibrium and non-equilibrium models presented in Sec. IV.

To summarize, we have shown that the EGP of purity-gapped fermionic Gaussian states with Bloch matrix representation $G_{k}$ satisfies $\varphi_{\text{E}}=\varphi_{\text{Z}}+\Delta(N)$, where $\varphi_{\text{Z}}$ is the Zak phase of the lowest purity band [given by Eq. (23)], and $\Delta(N)$ is a correction that vanishes in the thermodynamic limit. For equilibrium thermal states $\rho\propto e^{-\hat{H}_{k}/T}$, the Zak phase is equal to $2\pi$ times the zero-temperature ground-state polarization of $\hat{H}_{k}$, and the observable $\varphi_{\text{E}}$ probes this value even at temperatures $T\sim\Delta\epsilon$ of the order or higher than the characteristic band gaps $\Delta\epsilon$ in the system.

We emphasize that the temperature dependence $\sim(T/N)^{2}$ of corrections in $\varphi_{\text{E}}$ is fundamentally different from that in single-particle observables probing topological quantization. In the latter case, corrections generally scale exponentially ($\sim e^{-\Delta\epsilon/T}$) with temperature $T\lesssim\Delta\epsilon$, independently of system size (see, e.g., Ref. Rivas et al. (2013)), and crucially do not approach zero in the thermodynamic limit. The general mechanism identified here is also different from previous approaches focusing on Uhlmann-type phases Viyuela et al. (2014a, b); Huang and Arovas (2014); Budich and Diehl (2015): the latter are based on the construction of a system-size-independent geometric phase for density matrices, in contrast to the present construction where a gauge structure emerges only in the thermodynamic limit. Finally, we note that the scaling of corrections $\Delta(N)$ to $\varphi_{\text{E}}=\varphi_{\text{Z}}$ may be even more favorable in the presence of specific symmetries. An example is provided by the thermal density matrix of an SSH chain Su et al. (1979): in that case, the sublattice (chiral) symmetry of the system leads to a Berry connection $\mathcal{A}_{k}$ in Eq. (24) where one Pauli-matrix component is symmetry forbidden. Without going into detail, we mention that this symmetry leads to a correction $\Delta(N)\sim\exp(-\bar{\beta}N)$, exponential with system size.

III.3 Topological nature of the quantized pumping

We have argued in Sec. II that the EGP difference $\frac{1}{2\pi}\Delta\varphi_{\text{E}}=\oint d\phi\,\partial_{\phi}\varphi_{\text{E}}(\phi)$ per cycle in some parameter $\phi$ is quantized in integer multiples of $2\pi$, which is not a priori obvious to reconcile with the finding that $\varphi_{\text{E}}$ is given by a geometric (Zak) phase plus a perturbative correction $\Delta(N)$. Since the EGP is defined modulo $2\pi$, however, the loop integral $\oint d\phi\,\partial_{\phi}\varphi_{\text{E}}(\phi)=\varphi_{\text{E}}(\phi_{f})-\varphi_{\text{E}}(\phi_{i})$ (where $\phi_{i}$ and $\phi_{f}\equiv\phi_{i}$ are the start and end points of the parameter cycle) must indeed be quantized in units of $2\pi$ irrespective of $N$. To reveal the topological nature and, hence, the robustness of this quantization, we first recall the above result that

[TABLE]

where $A_{0}(k,\phi)$ is the Berry connection of the lowest purity band (which here depends on $\phi$). Since $\Delta\varphi_{\text{E}}$ is quantized irrespective of $N$, the correction $\Delta(N)$ cannot contribute to its value. Instead, $\Delta\varphi_{\text{E}}$ is determined by the winding of the Zak phase $\varphi_{\text{Z}}(\phi)=\oint_{\text{BZ}}dkA_{0}(k,\phi)$ as $\phi$ is varied from $\phi_{i}$ and $\phi_{f}$. This winding formally corresponds to an integer topological invariant known as the Chern number, and we can write

[TABLE]

where $F_{0}(k,\phi)\equiv\partial_{\phi}A_{k,0}(k,\phi)-\partial_{k}A_{\phi,0}(k,\phi)$ is the Berry curvature of the lowest purity band, defined in terms of the $U(1)$ gauge potential $A_{j,0}\equiv-i\left\langle u_{0}(k,\phi)\right\lvert\partial_{j}\left\lvert u_{0}(k,\phi)\right\rangle$ (where $\left\lvert u_{0}(k,\phi)\right\rangle$ are the Bloch vectors forming the lowest purity band). Eq. (32) shows that $\Delta\varphi_{\text{E}}$ is a topologically quantized integer which coincides with the Chern number $C$ of the lowest purity band — for any system size $N$ — which is one of the key results of this work.

III.4 Discussion

So far, we have illustrated our key gauge-reduction mechanism in a simple two-band model with a symmetric gapped purity spectrum ($\pm\beta_{k}\not=0$ for all $k$). Here, we consider the more general case of $n$ bands and examine the role of the chemical potential (i.e., of band filling) in the thermal setting. We recall that the purity spectrum is given by $\beta_{k,s}=\beta\epsilon_{k,s}$ for thermal states $\sim e^{-\beta\hat{H}}$, where $\epsilon_{k,s}$ is the energy spectrum of $\hat{H}$ (and $s$ is the band index). If we work in the grand-canonical ensemble and introduce a chemical potential $\mu$, the relevant states become $\sim e^{-\beta(\hat{H}-\mu)}$, and the corresponding purity spectrum reads $\beta_{k,s}=\beta(\epsilon_{k,s}-\mu)$. Therefore, purity eigenvalues $\beta_{k,s}$ are positive (negative) for states located below (above) the chemical potential. In turn, the weight matrix $e^{-B_{k}}=\text{diag}_{s}(e^{-\beta_{k,s}})$ controlling the gauge-reduction mechanism in Eq. (20) contains exponentially decreasing (increasing) diagonal elements for states located below (above) the chemical potential. We can thus distinguish two cases depending on whether the chemical potential lies (i) within a gap, or (ii) inside a band (which would correspond to complete or partial filling, respectively, at $T=0$):

In case (i), the direct analog of Eq. (22) reads

[TABLE]

where $\prod^{\prime}_{s}$ denotes the product over bands $s$ located below the chemical potential (“filled” bands). Accordingly, the EGP becomes $\varphi_{\text{E}}=\sum^{\prime}_{s}\varphi_{\text{Z},s}+\Delta(N)$, where $\varphi_{\text{Z},s}$ is the Zak phase of the purity band $s$ and the sum runs over filled bands [as before, $\Delta(N)$ is a correction that vanishes in the thermodynamic limit]. The relevant topological invariant is then $\frac{1}{2\pi}\Delta\varphi_{\text{E}}=\sum^{\prime}_{s}C_{s}$, where $C_{s}$ is the Chern number of the band $s$ located below the chemical potential. Therefore, in the general case of multiple bands, $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ reduces to its zero-temperature analog in a similar way as in the two-band model detailed above.

In case (ii) where the chemical potential lies within a specific band $s^{\prime}$, the purity eigenvalues $\beta_{k,s^{\prime}}$ change sign at certain values of $k$. As a result, the weight factors $e^{-\beta_{k,s}}$ only partially amplify the gauge-field contribution of the band $s^{\prime}$, and $\varphi_{\text{E}}$ does not reduce to a sum of geometric (Zak) phases. In that case, as expected, $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ is not a topological invariant.

III.5 Measurement of $\Delta\varphi_{\text{E}}$ and “purity adiabaticity” requirement

In the conventional zero-temperature setting, where relevant states are pure (ground) states, topological order parameters can be determined by measuring currents integrated over a closed parameter cycle — as typically envisioned in solid-state setups Niu and Thouless (1984) — or, equivalently, by measuring the Zak-phase difference accumulated over a cycle — as done in experiments with ultracold atoms Atala et al. (2013). Such measurements rely on a dynamical notion of adiabaticity, where pump parameters must be varied slowly in time as compared to the timescale set by some relevant gap (typically, the Hamiltonian gap).

Here we show that the requirements for measuring the mixed-state topological order parameter defined by the accumulated EGP difference $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ is more naturally related to a “purity adiabaticity” criterion. To this aim, we propose to determine the topological invariant $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ from a set of $M$ independent measurements of EPG values $\varphi_{\text{E}}(\phi_{j})$ at a discrete set of points $\phi_{j}$ along some relevant parameter cycle $\phi\in[0,2\pi]$. The purity adiabaticity condition expressed in terms of the dimensionless purity gap $\Delta\beta$ and the dimensionless “sampling rate” (or inverse “mesh size”) $\Delta\phi\equiv 1/M$ along the cycle (both assumed to be constant, for simplicity) then reads

[TABLE]

to be contrasted to the usual dynamical adiabaticity criterion $\dot{\phi}\ll\Delta\epsilon$ relating the rate of parameter changes to an energy or damping gap (see also Appendix C).

To derive the above criterion, we examine how to extract the integer-quantized topological invariant $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ from a generically imperfect set of $M$ distinct EGP measurements $\varphi_{\text{E}}(\phi_{j})$ along the relevant cycle in $\phi$ (where $j=1,\ldots,M$). Following the approach of Ref. Fukui et al. (2005), we discretize the integral $\frac{1}{2\pi}\Delta\varphi_{\text{E}}=1/(2\pi)\oint d\phi\,\partial_{\phi}\varphi_{\text{E}}$ in a way that crucially preserves two key properties: (i) the gauge invariance of $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$, and (ii) its integer quantization. Specifically, we define the $U(1)$ “link variables” $U(\phi_{j})\equiv\exp(i[\varphi_{\text{E}}(\phi_{j+1})-\varphi_{\text{E}}(\phi_{j})])$ and the corresponding “lattice field strengths” $F(\phi_{j})\equiv\mathop{\mathrm{Ln}}U(\phi_{j})$, where “$\mathop{\mathrm{Ln}}$” denotes the principal branch of the logarithm defined such that $-\pi<F(\phi_{j})/i\leq\pi$. We then estimate the topological invariant of interest as the sum $\frac{1}{2\pi}\Delta\varphi_{\text{E}}^{\prime}\equiv 1/(2\pi i)\sum_{j}F(\phi_{j})$. Clearly, this quantity is invariant under gauge transformations $\varphi_{\text{E}}(\phi_{j})\to\varphi_{\text{E}}(\phi_{j})+2\pi n_{j}$ (where $n_{j}$ is an arbitrary integer), and $\Delta\varphi_{\text{E}}^{\prime}\to\Delta\varphi_{\text{E}}$ as $M\to\infty$. Remarkably, $\frac{1}{2\pi}\Delta\varphi_{\text{E}}^{\prime}$ is additionally restricted, by construction, to integer values Fukui et al. (2005). As a result, one finds that $\Delta\varphi_{\text{E}}^{\prime}=\Delta\varphi_{\text{E}}$ as long as the mesh size $M$ (or number of sampling points in parameter space) is larger than a critical size $M_{c}$. In fact, $\frac{1}{2\pi}\Delta\varphi_{\text{E}}^{\prime}$ can only change (i.e., jump by an integer value) when $|F(\phi_{j})|=\pi$ at some point $j$ in parameter space, which corresponds to a large discontinuity $|\varphi_{\text{E}}(\phi_{j+1})-\varphi_{\text{E}}(\phi_{j})|=\pi$ [modulo $2\pi$, as $2\pi$ jumps do not contribute to $U(\phi_{j})]$. Accordingly, the critical mesh size can be estimated as the size below which the “admissibility condition” $|F(\phi_{j})|<\pi$ (for all $j$) breaks down.

In summary, the value $\frac{1}{2\pi}\Delta\varphi_{\text{E}}^{\prime}$ extracted from independent EGP measurements via the above procedure exactly coincides with the integer topological invariant $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ provided that $\varphi_{\text{E}}(\phi)$ is measured at a sufficiently large number of points $M>M_{c}$. In general, the critical mesh size $M_{c}$ is controlled by the proximity of the cyclic path $\phi\in[0,2\pi]$ to gap closing points (see, e.g., Ref. Dauphin et al. (2014)): the latter can be seen as sources of Berry-type curvature, in the sense that the field strength $F(\phi)$ is concentrated at such points in the limit of an infinitesimal mesh $M\to\infty$. Here, the relevant gap is the purity gap. Indeed, as we have demonstrated in Sec. IV, the EGP winding $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ vanishes for parameter cycles that do not encircle one (or more) purity-gap-closing point(s). This allows us to define the above notion of “purity adiabaticity” unique to thermal and nonequilibrium systems: to be able to observe the topological invariant $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$, one must sample a number $M>M_{c}$ of EGP values which gets larger and larger as one approaches purity-gap-closing points — diverging exactly at such points 666The critical mesh size $M_{c}$ increases with the absolute value of the topological invariant $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ [i.e., with the amplitude of the curvature concentrated at the gap-closing point(s)]. The critical value $M_{c}$ is typically not large in systems where $\Delta\varphi_{\text{E}}\sim\mathcal{O}(1)$ Dauphin et al. (2014).. This leads to the criterion presented in Eq. (34).

We emphasize that measurement errors on the discrete values $\varphi_{\text{E}}(\phi_{j})$ are irrelevant as long as the admissibility condition $|F(\phi_{j})|<\pi$ (for all $j$) remains satisfied (recall that $\Delta\varphi_{\text{E}}^{\prime}$ cannot change without breaking this condition). Therefore, errors can generically be compensated for by (i) using a finer mesh, or/and (ii) choosing parameter cycles further away from purity-gap-closing point(s).

Unlike usual measurements of accumulated Zak phase differences in the zero-temperature setting, the above procedure for measuring $\Delta\varphi_{\text{E}}$ does not rely on any dynamical protocol. This provides intuition as to why the purity adiabaticity criterion in Eq. (34) involves a comparison of dimensionless numbers instead of dynamical scales. Since the values $\varphi_{\text{E}}(\phi_{j})$ can be determined via completely independent measurements, the system can always be prepared with fixed parameters and measured after the time required for reaching its stationary state (controlled by possibly complex thermalization processes, in the thermal Hamiltonian case, or by a given damping gap, in the nonequilibrium Liouvillian case) — with otherwise no requirement for adiabaticity under dynamical changes of parameters. For completeness, however, we present in Appendix C a detailed analysis of dynamical quasi-adiabatic measurements of $\Delta\varphi_{\text{E}}$, where parameters are varied continuously in time. The advantage of such measurements as opposed to independent ones as above is that the state of the system follows the quasi-adiabatic evolution of parameters, which naturally fixes the gauge and leads to continuous changes in $\varphi_{\text{E}}(\phi)$. At the end of the parameter cycle, the topological invariant $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ is simply given by $\frac{1}{2\pi}|\varphi_{\text{E}}(\phi=2\pi)-\varphi_{\text{E}}(\phi=0)|$. The downside, however, is that this value generically deviates from an integer, due to the dynamical errors that come into play (see Appendix C).

IV Equilibrium (thermal) and non-equilibrium examples

In this section, we demonstrate our analytical results numerically in two illustrative examples: (i) the Rice-Mele model in thermal equilibrium, and (ii) its nonequilibrium driven-dissipative analog introduced in Ref. Linzner et al. (2016). Both models are noninteracting and translationally invariant. They exhibit Gaussian states $\rho\sim e^{-G}$ [Eq. (8)], described by a “fictitious Hamiltonian” $G$ (or $G_{k}$, in momentum space), and the EGP can be computed, e.g., using the path-ordered formula found in Eq. (19) and (20). We will illustrate three key features: (i) the convergence, in limit of large system sizes, of the EGP $\varphi_{\text{E}}$ to the Zak phase $\varphi_{\text{Z}}$ of the lower band of $G_{k}$ (the lower purity band), (ii) the quantization of the EGP difference $\Delta\varphi_{\text{E}}$ accumulated over a closed cycle in parameter space, and the coincidence of $\tfrac{1}{2\pi}\Delta\varphi_{\text{E}}$ with the Chern number of the lower purity band, and (iii) the direct connection between purity-gap-closing and topological transitions in $\tfrac{1}{2\pi}\Delta\varphi_{\text{E}}$.

We first examine the Rice-Mele model Rice and Mele (1982), defined by the Hamiltonian

[TABLE]

where $r=0,\ldots,N-1$ indexes unit cells and $s=0,1$ indexes fermionic sites in the unit cell. The first line describes the hopping of fermions on a 1D lattice with alternating hopping amplitudes $t_{1}$ and $t_{2}$, and the second line describes a staggered potential. At $\Delta=0$, the model reduces to the SSH model Su et al. (1979). It exhibits chiral symmetry, which promotes the Zak phase to a topological invariant. In that case, two topologically distinct phases (protected by chiral symmetry) can be distinguished for $t_{1}>t_{2}$ and $t_{1}<t_{2}$, respectively (separated by a gapless point at $t_{1}=t_{2}$). The corresponding quantized values of the Zak phase are $\varphi_{\text{Z}}=0$ and $\pi$ (modulo $2\pi$), respectively, which corresponds to ground-state polarizations $P=\varphi_{\text{Z}}/(2\pi)=0$ and $1/2$. In the Rice-Mele model, the parameter $\Delta$ provides a way to break chiral symmetry and, hence, to adiabatically connect the two phases originating from the SSH model and induce quantized polarization changes $\Delta P$ (note that the two phases are not symmetry protected anymore when $\Delta\neq 0$). In particular, adiabatic cycles in parameter space $(t_{1}-t_{2},\Delta)$ lead to an integer-quantized polarization difference (i.e., a pumped charge) $\Delta P=1$ whenever the gapless point $t_{1}=t_{2}$, $\Delta=0$ is encircled. This process corresponds to a topological (Thouless) pump.

We now examine the behavior of the Rice-Mele model at finite temperature where the EGP replaces the Zak phase as the relevant probe for topology. The fictitious Hamiltonian representing the thermal state $\rho\sim e^{-\beta H_{k}}\equiv e^{-G_{k}}$ of the system is given by $G_{k}=\beta H_{k}$, where $H_{k}$ is the momentum-space Hamiltonian matrix $H$ corresponding to Eq. (35). It can be expressed in the form

[TABLE]

where $\bm{\sigma}=(\sigma_{1},\sigma_{2},\sigma_{3})^{T}$ is a vector of Pauli matrices. The matrix $G_{k}$ is diagonalized by the unitary matrices $U_{k}$, and its spectrum (the purity spectrum) takes the form $\pm\beta_{k}$ with $\beta_{k}=\beta\epsilon_{k}$, where $\epsilon_{k}$ is the energy spectrum of the underlying Hamiltonian $H_{k}$ [using similar notations as before, as in Eq. (15)].

We have used the above representation as input for the numerical evaluation of Eqs. (19) and (20) [which provide an exact reformulation of Eq. (7) for the EGP]. In Fig. 2, we plot the difference $|\varphi_{\text{E}}-\varphi_{\text{Z}}|$ between the EGP and the Zak phase for fixed system parameters $t_{1}-t_{2},\Delta\neq 0$, over a wide range of $\beta$ including temperatures much larger than the Hamiltonian gap (of order $1$ for the chosen parameters). For all but the largest values of $T$, the data confirms the scaling $\sim 1/N^{2}$ predicted by perturbation theory [see Eq. (30)].

Next, we consider the EGP difference $\Delta\varphi_{\text{E}}$ accumulated over a closed path in parameter space $(t_{1}-t_{2},\Delta)$ encircling the origin (gapless point of $G_{k}=\beta H_{k}$). As argued in previous sections, we expect $\tfrac{1}{2\pi}\Delta\varphi_{\text{E}}$ to be an integer equal to the topologically quantized change (Chern number) $\tfrac{1}{2\pi}\Delta\varphi_{\text{Z}}$ of the Zak phase of the lowest band of $G_{k}=\beta H_{k}$ over the same parameter cycle. This equality must hold regardless of the system size $N$ and temperature $T$. The inset of Fig. 2 confirms this behavior: as the temperature is increased away from zero ($\beta$ decreased away from $\infty$), the difference $\tfrac{1}{2\pi}\Delta\varphi_{\text{E}}$ remains quantized for all system sizes $N$ accessible numerically. Moreover, its value indeed coincides with the value $\tfrac{1}{2\pi}\Delta\varphi_{\text{Z}}=1$ corresponding to the quantized charge $\Delta P=1$ that would be pumped through the system at $T=0$. The topological quantization of $\tfrac{1}{2\pi}\Delta\varphi_{\text{E}}$ requires the spectral gap of $G_{k}=\beta H_{k}$ (the purity gap) to be finite all along the chosen cycle in parameter space (as required for “purity adiabaticity”; see Sec. III.5). More importantly, the nontrivial value of $\tfrac{1}{2\pi}\Delta\varphi_{\text{E}}$ crucially depends on the existence of purity-gap-closing points encircled by the parameter cycle. In the thermal setting of interest here, the spectral gap of $G_{k}=\beta H_{k}$ can close either (i) via the closure of the energy gap $\Delta\epsilon$ of the underlying Hamiltonian $H_{k}$, or (ii) at infinite temperature where $\beta\to 0$. This leads to two possibilities for topological phase transitions. In the inset of Fig. 2, the value of $\tfrac{1}{2\pi}\Delta\varphi_{\text{E}}$ computed for a range of negative to positive temperatures 777We consider negative temperatures for the sake of this argument. illustrates these two possibilities: (i) when $|\beta|\neq 0$, the value of $\tfrac{1}{2\pi}\Delta\varphi_{\text{E}}$ is nontrivial because the parameter cycle that we consider encircles the purity-gap-closing point at the origin in parameter space $(t_{1}-t_{2},\Delta)$. It would be zero otherwise. (ii) When going from positive to negative temperatures, a topological transition occurs where $\tfrac{1}{2\pi}\Delta\varphi_{\text{E}}$ changes sign. Intuitively, the reason for this “jump” is that at positive/negative temperatures, the lower/upper band is predominantly occupied. At $\beta=0$, occupation inversion occurs, and the sign of $\varphi_{\text{E}}$ changes. We note that negative temperatures strictly speaking do not correspond to an equilibrium situation anymore, although the generator of dynamics still is a Hamiltonian operator alone.

We now turn to a second illustrative example provided by the nonequilibrium analog of the Rice-Mele model introduced in Ref. Linzner et al. (2016). In that case, the relevant dynamics is governed by a gapped Liouvillian (Markovian quantum master equation) of the generic form

[TABLE]

where $\rho$ is the density matrix of the system and $L_{r,s}$ are so-called Lindblad operators. This type of time evolution can be realized, e.g., in quantum systems with an engineered system-bath coupling. Provided that the bath energies are separated from those of the system by a large energy scale, Born-Markov and rotating-wave approximations become applicable, and lead to Lindblad master equations as above (see, e.g., Ref. Bardyn et al. (2013) for details). Here we assume, for simplicity, that the right-hand side of Eq. (37) does not contain a coherent contribution $\sim-i[\hat{H},\rho]$.

The nonequilibrium Rice-Model model analog of interest is defined by the set of Lindblad operators

[TABLE]

where $L_{r,0}$ and $L_{r,1}$ act inside and between unit cells, respectively, and in this regard are analogous to the two hopping terms in Eq. (35). Referring to Ref. Linzner et al. (2016) for details, we note that these operators are defined such that, on a timescale $\sim 1/\Delta_{\text{d}}$ set by the so-called “damping gap” $\Delta_{\text{d}}$ of the Liouvillian $\mathcal{L}$ Bardyn et al. (2013), the dynamics drives an arbitrary initial Gaussian state to a specific stationary Gaussian state $\rho$ satisfying $\mathcal{L}(\rho)=0$. The latter does not obey strict particle number conservation. Its mean particle number $\langle\hat{n}\rangle$, however, is stationary, and fluctuations $\langle\delta\hat{n}\rangle$ around it are intensive, i.e., $\langle\delta\hat{n}\rangle/\langle\hat{n}\rangle\to 0$ in the thermodynamic limit. For all practical purposes, the Gaussian state $\rho\sim e^{-G_{k}}$ can then be described by a number-conserving fictitious Hamiltonian $G_{k}$. Due to the specific form of the Lindblad operators, the matrix $G_{k}$ exhibits the same structure as the Hamiltonian matrix $H_{k}$ of the Rice-Mele model. Specifically, $G_{k}$ is given by Eq. (36) with $\beta=1$ (no notion of temperature here), $t_{1}=\tfrac{1}{4}(1+\epsilon)(\lambda^{2}-1)/(\lambda^{2}+1)$, $t_{2}=\tfrac{1}{4}(1-\epsilon)(\lambda^{2}-1)/(\lambda^{2}+1)$, and $\Delta=\lambda/(\lambda^{2}+1)$. Since $t_{1}-t_{2}\propto\epsilon$ and $\Delta\propto\lambda$, the real parameters $(\epsilon,\lambda)$ play a similar role as the parameters $(t_{1}-t_{2},\Delta)$ in the Rice-Mele model. In particular, the origin $\epsilon=\lambda=0$ corresponds to a purity-gap-closing point. As in the thermal case, the EGP of the stationary state can be computed using Eqs. (19) and (20).

In Fig. 3, we show the computed difference $|\varphi_{\text{E}}-\varphi_{\text{Z}}|$ between the EGP and the Zak phase of the lowest purity band (lowest band of $G_{k}$). As in the thermal case, we verify that $\varphi_{\text{E}}\to\varphi_{\text{Z}}$ in the thermodynamic limit, with power-law scaling $\sim N^{-2}$. Purity-gap-closing points play the same key role here: in particular, we observe a topological transition where the EGP difference $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ accumulated over a complete cycle in $(\epsilon,\lambda)$ parameter space “jumps” from trivial ([math]) to nontrivial ($1$) when varying a parameter $\delta\equiv\delta(\epsilon,\lambda)$ controlling whether the purity-gap-closing point at $\epsilon=\lambda=0$ is encircled [see insets of Fig. 3]. In equilibrium thermal systems, purity-gap-closing points necessarily coincide with points at which the gap of the system Hamiltonian closes, due to the tight correspondence $G_{k}=\beta H_{k}$. In nonequilibrium systems, in contrast, purity-gap-closing points need not coincide with points where the gap of the system Liouvillian (nonequilibrium analog of a Hamiltonian) closes. In fact, here, the gap of the Liouvillian (the damping gap) is given by $\Delta_{\text{d}}=4[1+\lambda^{2}(2\epsilon^{2}+\lambda^{2}+2(\epsilon^{2}-1)\cos k)]^{1/2}$ 888We typically consider parameters with $|\lambda|,|\epsilon|<1$, to avoid gap-closing points $|\lambda|=\pm 1$, $\epsilon=0$.. Clearly, it does not close at the purity-gap-closing point $\epsilon=\lambda=0$. The damping gap plays a similar role as a Hamiltonian gap, in the sense that it ensures the exponential decay of spatial correlations Bardyn et al. (2013). Therefore, here, the fact that it remains open at the point $\epsilon=\lambda=0$ where the purity gap closes and a topological transition occurs exemplifies a remarkable possibility unique to nonequilibrium systems: the fact that topological transitions can occur without the appearance of divergent length or time scales (see Ref. Bardyn et al. (2013) for a detailed discussion).

V Measurement of the EGP

As shown above, the EGP of a thermal or nonequilibrium state is not related to particle currents, which implies that it cannot be measured via particle transport. In the following, we propose an interferometric scheme to detect it, building on an idea by Sjöqvist et al. Sjöqvist et al. (2000). To this end, we recall that $\varphi_{\text{E}}=\mathop{\mathrm{Im}}\ln\langle\hat{T}\rangle$ is nothing but the argument of the complex-valued observable $\langle\hat{T}\rangle$, i.e., $\varphi_{\text{E}}=\mathrm{arg}\langle\hat{T}\rangle$, where $\hat{T}=\exp(i\delta k\sum_{i}x_{i}\hat{a}_{i}^{\dagger}\hat{a}_{i})$ [recall that $x_{i}$ denotes the position of fermions on site $i\equiv(r,s)$, with creation and annihilation operators $\hat{a}_{i}^{\dagger},\hat{a}_{i}$].

We consider a Mach-Zehnder interferometer whose lower arm contains the system to be probed [Fig. 4(a)] — fermions in a 1D lattice of length $L$ corresponding to a few tens of sites, as in typical setups with ultracold atoms. Photons moving along the two directions defined by the interferometer geometry can be represented by a two-state wavefunction with upper and lower components describing photons propagating in the “vertical” and “horizontal” directions, respectively. The action of mirrors and beam splitters is then described by unitary $2\times 2$ matrices

[TABLE]

respectively. We assume that each lattice site $i$ consists of two internal fermionic levels (“ground” and “excited”) corresponding to annihilation operators $\hat{a}_{i}$ and $\hat{b}_{i}$, respectively. As illustrated in Fig. 4b, fermions are coupled to a photonic mode with carrier frequency $\omega_{0}=k_{0}z$ (where $z$ is the position along the propagation path of interest), described by the annihilation operator $\hat{c}(z)$ satisfying the commutation relation $[\hat{c}(z),\hat{c}^{\dagger}(z^{\prime})]=\delta(z-z^{\prime})$. We assume that this mode couples the internal ground and excited states with coupling constant $g$ and a large detuning $\Delta$. When $\Delta$ is larger than all other relevant energy scales, internal excited states can be adiabatically eliminated, leading to the effective Hamiltonian

[TABLE]

where $f_{j}$ is the complex mode function of the photonic mode $\hat{c}(z)$ at site $j$. This Hamiltonian describes forward Brillouin scattering.

Next, we assume that the amplitude of the mode $f_{j}$ exhibits a spatial gradient along the axis of the probed system (realized, e.g., by a TEM01 Gauss-Hermite mode), such that $|f_{j}|^{2}\propto x_{j}$ (recall that $x_{j}$ is the position of the fermionic site $j$). Denoting $|f_{j}|^{2}\equiv\eta\,x_{j}/L$, Eq. (41) then describes the coupling between photons and the center-of-mass position operator $\hat{X}=\sum_{j}x_{j}\hat{a}_{j}^{\dagger}\hat{a}_{j}$ of the probed system, as desired. The propagation of the photonic mode along the lower arm of the interferometer is described by the equation

[TABLE]

(setting the speed of light $c=1$), with solution

[TABLE]

where $z$ lies beyond the region where photons interact with the probed system. By adjusting the detuning $\Delta$ such that $g^{2}\eta/(L\Delta)=2\pi/L\equiv\delta k$, photons in the lower arm of the interferometer pick up a phase proportional to the center of mass of fermions in the probed system — described by the unitary transformation $\hat{T}=\exp(-i\delta k\hat{X})$, as desired.

We remark that photons additionally experience a momentum “kick” in the direction of the lattice of the probed system due to the spatial gradient in their mode function. The kick imparted to each photon, however, is less than their initial momentum $k_{0}$. This effect thus only leads to a small diffraction which we neglect here, for simplicity (though it should to be taken into account when designing an actual experiment). We note that the prefactor $g^{2}\eta/(\Delta L)$ can be increased, at fixed detuning $\Delta$, by adding a build-up cavity to the lower arm of the interferometer to make photons bounce back and forth through the probed system before propagating further — thereby enhancing the effective interaction between photons and fermions.

The total unitary matrix describing the propagation of a photon through the interferometer reads

[TABLE]

where $\mathbf{U}_{B}=U_{B}\otimes\mathbb{1}$ and $\mathbf{U}_{M}=U_{M}\otimes\mathbb{1}$ ($\mathbb{1}$ being the unit operator in the Hilbert space of the fermions), and

[TABLE]

where $\chi$ is a tunable phase in the upper arm of the interferometer [Fig. 4(a)]. Overall, the interferometer transforms the input state $\varrho_{\text{in}}=\rho_{\text{ph}}\otimes\rho$, where $\rho_{\text{ph}}$ is the input state of photons and $\rho$ is the mixed state of the probed system, into $\mathbf{U}\varrho_{\text{in}}\mathbf{U}^{\dagger}$. The resulting intensity at the two output ports (“$+$” and “$-$”) is then given by

[TABLE]

where $\hat{n}$ denotes photon number operators and the expectation value $\langle\hat{T}\rangle$ is only over the fermionic state of the probed system. Therefore, in the above setup, photons injected into the lower arm of the interferometer pick up a phase which crucially corresponds to the EGP $\varphi_{\text{E}}=\mathrm{arg}\langle\hat{T}\rangle$. Most importantly, this phase can be measured by monitoring directly the intensity difference \Delta\hat{n}\equiv\hat{n}_{\text{out}}^{+}-\hat{n}_{\text{out}}^{-}=|\langle\hat{T}\rangle|\cos\bigl{[}\chi-\mathrm{arg}\langle\hat{T}\rangle\bigr{]}\hat{n}_{\mathrm{in}} between outputs (balanced detection) as a function of the reference phase $\chi$ set in the upper arm. To extract the accumulated phase difference $\Delta\varphi_{\text{E}}=\oint d\phi\,\partial_{\phi}\varphi_{\text{E}}$, one can then (i) repeat the experiment for multiple parameter values $\phi$ along a cyclic path, and (ii) follow the procedure discussed in Sec. III.5 to extract the exact integer value of $\Delta\varphi_{\text{E}}$ in a way that is gauge-invariant and, most importantly, robust against (small) measurement errors on $\varphi_{\text{E}}$.

The visibility of the EGP signal resulting from the above balanced detection scheme is unity. Since $|\langle\hat{T}\rangle|$ is typically small, however, the characteristic number of output photons per input photon is small and one may have to accumulate photons for a longer time to reach the desired sensitivity in the measurement of $\varphi_{\text{E}}$. The minimum detectable phase is set by shot noise \varphi_{\text{E}}\bigr{|}_{\mathrm{min}}=\Delta\varphi_{\mathrm{shot}}\sim 1/\sqrt{P_{\mathrm{out}}\tau}, where $P_{\mathrm{out}}$ is the maximum output flux of photons per unit time and $\tau$ is the overall measurement time.

In the illustrative two-band model with purity spectrum $\pm\beta_{k}$ examined above, one finds $|\langle\hat{T}\rangle|=\exp[-\tfrac{N}{2\pi}\int dk\ln(1+e^{-\beta_{k}})+\mathcal{O}(1/N^{2})]\approx\exp(-\tfrac{N}{2\pi}e^{-\Delta\beta/2})$ (with an additional factor $\sin\theta$), where $\Delta\beta=2\min_{k}\beta_{k}$ is the purity gap 999Note that the real factor $\ln\det[1-f(G)]$ neglected in Eq. (19) must be taken into account in $|\langle\hat{T}\rangle|$.. Note that also $|\langle T\rangle|\approx 1$ can be achieved e.g., in a mesoscopic thermal system of $N\approx 50$ sites with a purity gap $\Delta\beta=\Delta\epsilon/T\approx 5$ (where $\Delta\epsilon$ is the energy gap of the underlying Hamiltonian).

We finally comment on the effects of boundary conditions: although we have always considered a 1D lattice system of fermions with periodic boundary conditions, for pedagogical purposes, we emphasize that the EGP is a bulk quantity (reducing to the polarization, in the zero-temperature limit) which is, hence, essentially unaffected by boundary conditions. In particular, open boundary conditions can also be used, as implicitly assumed in the above measurement setup. In that case, correlations $\langle\hat{a}^{\dagger}_{i}\hat{a}_{j}\rangle=[f(G)]_{ij}$ between sites $i$ and $j$ located at opposite ends of the system — corresponding to the corner elements of the matrix $f(G)$ — are essentially removed 101010Correlations correspond to the state defined in Eq. (8), which can be seen as the thermal density matrix of a fictitious Hamiltonian $G$ (with temperature $\beta=1$). Accordingly, correlations must decay exponentially with system size over a length scale proportional to the inverse spectral gap of $G$ (the purity gap) Hastings (2004). As in the rest of this work, we assume that the purity gap is finite., and the matrix $f(G)(T-1)$ appearing in Eq. (18) for the EGP (namely, $\varphi_{\text{E}}=\mathop{\mathrm{Im}}\ln\det[\mathbb{1}+f(G)(T-\mathbb{1})]$) is similarly quasi-diagonal in position space, with no corner elements. Since $\det[\mathbb{1}+f(G)(T-\mathbb{1})]=\exp\mathop{\mathrm{Tr}}\ln[\mathbb{1}+f(G)(T-\mathbb{1})]$, elements of $f(G)(T-\mathbb{1})$ contribute to $\varphi_{\text{E}}$ to order $m$ (in a series expansion of the logarithm) whenever they contribute to diagonal elements in $[f(G)(T-\mathbb{1})]^{m}$. Here, quasi-diagonal elements in $f(G)(T-1)$ contribute at order $\sim 1$, while corner elements contribute at order $\sim N$ (where $N$ is the total number of sites in the system). As a result, corner elements determining boundary conditions lead to corrections to the EGP that are exponentially small with system size, as we have verified numerically. In particular, the plots presented in Sec. IV are visually unchanged when considering open instead of periodic boundary conditions.

VI Conclusions and outlook

We have shown that density matrices describing mixed fermionic Gaussian states in one dimension encode topological information in a way that enables a direct interpretation in terms of a physical observable. The connection to observables is provided by the ensemble geometric phase (EGP), which is defined for arbitrary density matrices — equilibrium and non-equilibrium states alike. The EGP is constructed from the expectation value of a many-particle operator: the operator of translations in momentum space by the smallest possible step $\delta k=2\pi/L$, which crucially equips us with a geometric notion of parallel transport for state vectors in that space. The mechanism underlying the robustness of the EGP as a geometric phase for mixed states is based on the statistical selection of the most strongly occupied Bloch band, i.e., typically the lowest one. This selection is a many-body effect related to the presence of $N$ modes in each Bloch band, leading to an effective physical purification of the selected Bloch band in the thermodynamic limit $N\to\infty$.

Although the price to pay to see fingerprints of topology in mixed states is the many-body character of the EGP, we have shown that the latter can be detected in interferometric measurements, e.g., in current setups based on cold atomic gases. These results have two important physical implications: conceptually, they demonstrate that topological phase transitions persist to finite temperatures, or, more generally, in mixed quantum states. More precisely, the degeneracy in the spectrum of the density matrix, measured by the purity-gap closing, is associated with a singularity in the EGP, i.e., with a jump in the accumulated EGP differences upon enclosing such a point in parameter space. Practically, the EGP provides a viable in situ alternative to detecting topological order in fermionic systems of ultracold atoms in low-temperature states and at finite density. So far, the Zak phase has been determined in the single-particle limit only, where particle statistics is irrelevant, by propagating a test particle through an otherwise empty band structure Atala et al. (2013) (see Ref. Jotzu et al. (2014) for a related strategy in a two-dimensional system). Cooling fermions to extremely low temperatures to access topological ground-state properties remains an outstanding challenge at finite fermion density. Due to its robustness towards finite temperatures, however, the EGP examined here could be used as a direct detection tool in such systems.

Our construction based on a unitary translation operator may lend itself to generalizations to interacting systems Shapourian et al. (2017); Shiozaki et al. (2017) in mixed states. Here we have focused on noninteracting, translation-invariant systems, and the crystal momentum may have seemed to play a very fundamental role. We emphasize, however, that the truly relevant object is the many-body translation operator $\hat{T}=e^{i\delta k\hat{X}}$ used to define the EGP. This operator describes a shift $\delta k=2\pi/L$ of the physical momentum of all particles, irrespective of interactions or disorder. It acts as a canonical transformation $\hat{p}_{j}\to\hat{T}^{\dagger}\hat{p}_{j}\hat{T}=\hat{p}_{j}+\delta k$ (where $\hat{p}_{j}$ is the momentum operator of particle $j$), which can conveniently be seen as the insertion of a magnetic flux $\Phi=2\pi$ through the periodic system Niu and Thouless (1984). Therefore, in the general case where the crystal momentum $k$ is not a good quantum number, $\Phi$ simply replaces $k$ as the relevant parameter for state vectors (as in many-body generalizations of the Zak phase and of the Chern number Niu and Thouless (1984)). In that case the operator $\hat{T}$ defines, for vectors $\left\lvert\psi_{\Phi}\right\rangle$, a similar notion of parallel transport in $\Phi$ space as it does in $k$ space for vectors $\left\lvert\psi_{k}\right\rangle$ in noninteracting systems with translation invariance. As a result, we expect the gauge-reduction mechanism identified in this work to hold in interacting or disordered systems, as long as the density matrix $\rho$ of the system (or, more precisely, the corresponding “fictitious Hamiltonian” $-\ln\rho$) has a gapped ground state. In particular, we expect the EGP to remain nontrivial when weak interactions and/or disorder preserving the purity gap are added to the examples examined in Sec. IV.

Beyond interacting and disordered systems, several directions will be exciting to explore: First, the gauge-reduction mechanism identified here seems very generic. In particular, we anticipate that one would obtain other well-defined geometric phases for mixed states by replacing the operator $\hat{T}$ of translations in momentum space by other unitary operators, generating translations in a different space (at the expense, however, of possibly losing the direct connection to physical observables which is a key feature of the EGP). Second, and more broadly, an intriguing direction for future research will be to examine whether other types of gauge structures can be extracted from many-body correlators Shapourian et al. (2017); Shiozaki et al. (2017), to construct topological classifications of mixed states. Extensions to systems with more than one spatial dimension will also be interesting to explore. Third, it will be exciting to ask whether the gauge-reduction mechanism identified in fermionic systems here can also play a role in bosonic ones. While bosons at equilibrium and low temperature tend to condense with a strong occupation of low-momentum modes, experiments based, e.g., on ultracold atoms in modulated potentials Aidelsburger et al. (2015) could lead to topological signatures in the EGP. Exciton polaritons are also potential candidates, as they routinely produce driven open quantum states with occupation properties reminiscent of fermionic systems Baboux et al. (2017).

Finally, our work may pave the way towards other probes of mixed-state topology. One promising candidate are Loschmidt amplitudes — another many-body observable related to the expectation value of a unitary matrix, namely, the time-evolution operator of Hamiltonian quantum dynamics. Cases which could perhaps be related to our mechanism were recently pointed out in Refs. Bhattacharya et al. (2017); Heyl and Budich (2017). Fingerprints of bulk topological properties at the edges of insulators and superconductors have been shown to persist at finite temperatures Quelle et al. (2016); Kempkes et al. (2016), and establishing a connection to the results presented here provides another challenge for future research.

Acknowledgements

We would like to thank D. Linzner for invaluable input and support, and J. C. Budich and M. Heyl for useful discussions. C.-E. B. gratefully acknowledges support from the DQMP at the University of Geneva and from the Swiss National Science Foundation under Division II. S. D. acknowledges support by the German Research Foundation (DFG) through the Institutional Strategy of the University of Cologne within the German Excellence Initiative (ZUK 81), as well as support by the European Research Council via ERC Grant Agreement n. 647434 (DOQS). L. W. and M. F. acknowledge support by the German Research Foundation (DFG) through the SFB TRR 185, and A. A. and S. D. (DFG) through the CRC 183 (project B02). C.-E. B., M. F. and S. D. would like to extend specials thanks to the KITP at UCSB for hospitality. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1125915.

Appendices

Appendix A EGP vs. adiabatic pump currents

In this appendix, we discuss the connection between the EGP and the current flow induced in an adiabatic pump protocol. The main conclusion will be that while $\tfrac{1}{2\pi}\Delta\varphi_{\text{E}}$ is quantized and physically observable, it does not correspond, for mixed states, to a quantized charge transfer. This is surprising inasmuch as the EGP does coincide with the adiabatic current in the zero-temperature (ground-state) limit. In the following, we start by reviewing the connection between Resta’s formula for the ground-state polarization of periodic quantum systems and the current flow.

A.1 Resta polarization and current flow

Resta formula for the ground-state polarization $P$ of a periodic quantum system [Eq. (4)] is constructed in such a way that differential changes $\partial_{t}P\equiv\dot{P}$ induced, e.g., by an external parameter are equal to the physical current $I$ (as they should). Here, we review how the connection between Eq. (4) and the current explicitly arises, in line with Resta’s arguments in Ref. Resta (1998).

The starting point of the construction is the observation that $\hat{T}|\psi_{0}\rangle=e^{i\delta k\hat{X}}|\psi_{0}\rangle$ corresponds, to first order in $\delta k=2\pi/L$, to the ground state of the momentum-shifted Hamiltonian $\hat{H}(q=-\delta k)\equiv\hat{T}\hat{H}\hat{T}^{\dagger}\simeq\hat{H}+i\delta k[\hat{X},\hat{H}]=\hat{H}-\delta k\hat{I}$, where $\hat{I}\equiv\partial_{q}\hat{H}(q)$ is the usual current operator. The action of $\hat{T}$ can be understood as a momentum shift of all particle momenta $\hat{p}_{j}\to\hat{T}\hat{p}_{j}\hat{T}^{\dagger}=\hat{p}_{j}-\delta k$. Equivalently, one can see $\hat{T}\hat{H}\hat{T}^{\dagger}$ as the Hamiltonian describing the periodic 1D system of interest after the adiabatic insertion of a single quantum of magnetic flux through the latter. To first order in $\delta k$, we obtain

[TABLE]

where $\varphi_{\text{Z}}$ is the Zak phase accumulated during the insertion of the flux quantum. In accordance with Resta’s formula [Eq. (4)], the geometric phase $\varphi_{\text{Z}}$ determines the instantaneous value of the polarization, namely, $P=\tfrac{1}{2\pi}\mathop{\mathrm{Im}}\ln\langle\psi_{0}|\hat{T}|\psi_{0}\rangle=\tfrac{1}{2\pi}\mathop{\mathrm{Im}}\ln\langle\psi_{0}|e^{i\varphi_{\text{Z}}}|\psi_{0}\rangle=\varphi_{\text{Z}}/(2\pi)$. The Zak phase, however, crucially does not enter the expression of the current flow, which involves excitations out of the ground state. Indeed, starting from Resta’s formula and using Eq. (47), the derivative $\dot{P}$ takes the form

[TABLE]

This expression was identified in Ref. Niu and Thouless (1984) as the pump current flowing in response to the adiabatic variation of external parameters. For the sake of completeness, we review this connection in the next subsection in a manner somewhat different from the original derivation, tailored to the present discussion. Readers wishing to proceed directly to the generalization to mixed states may skip this discussion.

A.2 For ground states, the time derivative of the EGP [Eq. (48)] coincides with the pump current

Let $\left\lvert\psi_{0}(t)\right\rangle\equiv\left\lvert\psi_{0}(\phi)\right\rangle$ be the instantaneous ground state at a value $\phi\equiv\phi(t)$ of a varied external parameter (such that $\hat{H}(\phi)\left\lvert\psi_{0}(\phi)\right\rangle=E_{0}(\phi)\left\lvert\psi_{0}(\phi)\right\rangle$). We assume that $\left\lvert\psi_{0}(\phi)\right\rangle$ is current-free, i.e., $\left\langle\psi_{0}(t)\right\lvert\hat{I}\left\lvert\psi_{0}(t)\right\rangle=0$. When the parameter $\phi$ is adiabatically varied, the time evolution of the system initially prepared in its ground state is described by the time-dependent Schrödinger equation $i\partial_{t}\left\lvert\psi\right\rangle=\hat{H}\left\lvert\psi\right\rangle$ (omitting explicit time dependences, for simplicity). A convenient ansatz consists in writing $\left\lvert\psi\right\rangle\equiv e^{-i(E_{0}t-\varphi)}\left\lvert\psi_{0}\right\rangle+\left\lvert\delta\psi\right\rangle\equiv\left\lvert\psi^{\prime}_{0}\right\rangle+\left\lvert\delta\psi\right\rangle$, with dynamical phase $E_{0}t$, geometric phase $\varphi\equiv\varphi(t)$, and out-of-ground-state excitations described by $\left\lvert\delta\psi\right\rangle$. Since the parameter $\phi$ is varied with frequency $\omega$ smaller than the excitation gap $\Delta\equiv\min_{n}(E_{n}-E_{0})$ (adiabacity condition), $\left\lvert\delta\psi\right\rangle$ will be small, but nevertheless important, as it is the wavefunction component responsible for the current flow, $\langle\hat{I}\rangle(t)\simeq\left\langle\psi^{\prime}_{0}(t)\right\lvert\hat{I}\left\lvert\delta\psi(t)\right\rangle+(\text{c.c.})$. Substitution of the ansatz into the Schrödinger equation yields

[TABLE]

To make progress, we expand $\left\lvert\delta\psi\right\rangle\equiv\sum_{n\neq 0}c_{n}e^{-iE_{n}t}\left\lvert\psi_{n}\right\rangle$ in terms of the instantaneous excited eigenstates of the Hamiltonian, where $c_{n}$ are the time-dependent coefficients to be solved for, with initial condition $c_{n}(0)=0$. In view of the expected smallness $c_{n}\sim\omega$, the temporal variation $|\dot{\psi}_{n}\rangle\sim\omega$ can be neglected as higher order, such that $(i\partial_{t}-H)\left\lvert\delta\psi\right\rangle\simeq\sum_{n}i\dot{c}_{n}e^{-iE_{n}t}\left\lvert\psi_{n}\right\rangle$. Substitution of this expression into Eq. (49) followed by a projection onto $\left\lvert\psi_{0}\right\rangle$ then identifies $\dot{\varphi}=i\langle\psi_{0}|\dot{\psi}_{0}\rangle$ as the Berry connection of the ground state (as expected). Projection onto excited states $\left\lvert\psi_{n}\right\rangle$, on the other hand, leads to the equation $\dot{c}_{n}=e^{i(E_{n}-E_{0})t}\langle\psi_{n}|\dot{\psi}_{0}\rangle$, where we have neglected a Berry-phase factor as higher order. Likewise neglecting the slow variation of the matrix elements $\langle\psi_{n}|\psi_{0}\rangle$ in comparison to the dynamical factor $e^{i(E_{n}-E_{0})t}$, we obtain $c_{n}\simeq i(E_{n}-E_{0})^{-1}\langle\psi_{n}|\dot{\psi}_{0}\rangle[1-e^{i(E_{n}-E_{0})t}]$. Inserting this solution (neglecting rapidly oscillatory factors) into the spectral decomposition of the current expectation value $\langle\hat{I}\rangle(t)\simeq\left\langle\psi^{\prime}_{0}(t)\right\lvert\hat{I}\left\lvert\delta\psi(t)\right\rangle+(\text{c.c.})$ finally leads to Eq. (48).

For completeness, we note that the current accumulated during a periodic pump cycle, $\Delta P=\oint dt\,\partial_{t}P=\oint d\phi\partial_{\phi}P$, is integer quantized. In Ref. Niu and Thouless (1984), this quantization was established in a three-step argument: first, the above derivation was generalized to include the presence of a general magnetic flux $\Phi\in[0,2\pi/L]=[0,\delta k]$ threading the system. Second, it was shown that the current expectation value is the same at any value of the flux, i.e., $\langle\hat{I}\rangle(\Phi)\simeq\langle\hat{I}\rangle(0)$ (up to corrections that vanish in the thermodynamic limit), such that the relevant current can be expressed as the average $\langle\hat{I}\rangle\equiv\delta k^{-1}\oint d\Phi\,\langle\hat{I}\rangle(\Phi)$. Third, it was argued that the presence of the flux leads to a modification of the Hamiltonian $\hat{H}\to\hat{H}+\Phi\hat{I}$, such that perturbation theory similar to the one outlined above yields $(E_{0}-E_{n})^{-1}\langle\psi_{0}|\hat{I}|\psi_{n}\rangle=\langle\psi_{0}|\partial_{\Phi}\psi_{n}\rangle$ to first order in $\delta k$ (where wavefunctions now depend on $\phi$ and $\Phi$). Using the additional identity $|\dot{\psi}_{0}\rangle=\dot{\phi}|\partial_{\phi}\psi\rangle$, the charge transported per adiabatic cycle becomes

[TABLE]

We recognize here the standard expression for the first Chern number of the Berry connection $A_{j}=i\langle\psi_{0}|\partial_{j}\psi_{0}\rangle$ (with $j=\phi,\Phi$), which shows that $\Delta P$ is indeed a topologically quantized integer.

A.3 For generic states, the time derivative of the EGP differs from the pump current

We have shown in the main text that the EGP is a physical observable and that its integrated change over a complete cycle in parameter space is integer quantized. We have also argued that the EGP reduces to the Zak phase ($2\pi$ times the polarization) in cases where the density matrix reduces to a ground-state projector. Nonetheless, its topological quantization cannot be interpreted in terms of current flow, as we demonstrate now.

For the sake of concreteness, consider a thermal density matrix with general spectral decomposition $\rho=\sum_{n}p_{n}|\psi_{n}\rangle\langle\psi_{n}|$, where $p_{n}=e^{-\beta E_{n}}/\mathcal{Z}$ and $\mathcal{Z}=\sum_{n}e^{-\beta E_{n}}$. When introducing a time-dependent parameter $\phi\equiv\phi(t)$, both $E_{n}(\phi)$ and the states $|\psi_{n}(\phi)\rangle$ become time dependent. The time derivative of the EGP defined by Eq. (7) then becomes

[TABLE]

This expression is nonlinear in the occupation numbers $p_{n}$ and, hence, does not lend itself to analytical simplifications. Assuming that parameter changes are slow enough for an adiabacity principle to hold for individual states with weight $p_{n}$, the matrix elements appearing in this expression will carry geometric phases [see Eq. (47)] which generally do not cancel out [except in the specific case of ground states $\rho=\left\lvert\psi_{0}\right\rangle\left\langle\psi_{0}\right\lvert$ where the above expression reduces to Eq. (48)]. Even for time-independent states $(\dot{p}_{n}=0)$, the presence of a nontrivial sum of operator expectation values in the denominator makes the above expression formally different from linear-response expectation values describing current flows.

Appendix B Perturbative corrections to the EGP

In this appendix, we provide additional details regarding the second-order perturbative expansion of the EGP [Eq. (26)] and the scaling of second-order corrections [Eq. (30)] in the same context as in Sec. III.2, i.e., in an illustrative two-band model with purity spectrum $\pm\beta_{k}$. For convenience, we focus on the limit of large system sizes $N$ where sums over momenta can be approximated as continuous integrals. Second-order corrections then read

[TABLE]

Here, upper and lower bands are identified by “$+$” and “$-$” indices, respectively, with Bloch states $|u_{k}^{\pm}\rangle$. The operators $W^{\pm}_{k_{1},k_{2}}=\exp[\int_{k_{1}}^{k_{2}}dk\,(iA_{k}^{\pm}\pm\delta k^{-1}\beta_{k})]\,\mathbb{1}$ describe the “evolution” in individual bands under the influence of the Berry connection $A_{k}^{\pm}=\mathcal{A}_{k}^{0}\pm\mathcal{A}_{k}^{3}=i\langle u_{k}^{\pm}|\partial_{k}u_{k}^{\pm}\rangle$ and the weight factors $\pm\beta_{k}$. The operators $V_{k}^{-+}=(\mathcal{A}_{k}^{1}-i\mathcal{A}_{k}^{2})\sigma^{+}$ and $V_{k}^{+-}=(\mathcal{A}_{k}^{1}+i\mathcal{A}_{k}^{2})\sigma^{-}$, on the other hand, are “jump operators” causing transitions between bands [with $\sigma^{\pm}\equiv(\sigma_{1}\pm i\sigma_{2})/2$]. The operator $G$ is the “Green’s function” defined in Eq. (26), which acts trivially in the upper band, and via the inverse evolution operator in the lower band [where $P^{\pm}=(\sigma_{0}\mp\sigma_{3})/2$ projects onto individual bands]. Using the above expressions, the second-order corrections can be cast in the form

[TABLE]

where the first and second lines are the contributions of $V^{(2)}$ and $(V^{(1)})^{2}$, respectively, and where we have used the properties $(W^{\pm}_{k_{1},k_{2}})^{-1}=W^{\pm}_{k_{2},k_{1}}$ and $W^{\pm}_{k_{1},k_{2}}W^{\pm}_{k_{2},k_{3}}=W^{\pm}_{k_{1},k_{3}}$. In equation (52), the second and fourth terms cancel out, and the third term is massively suppressed due to the presence of $W^{+}_{2\pi,0}\sim\exp(-\bar{\beta}N)$ with $\bar{\beta}\equiv(2\pi)^{-1}\int_{0}^{2\pi}dk\,\beta_{k}$. We are thus left with the first term, and noting that $(W^{-}_{k,k^{\prime}})^{-1}W^{+}_{k,k^{\prime}}=\exp[2\int_{k^{\prime}}^{k}dq(i\mathcal{A}_{q}^{3}-\delta k^{-1}\beta_{k})]$, we obtain

[TABLE]

This expression makes the essence of the EGP gauge-reduction mechanism manifest: due to the global presence of weight factors $\sim\exp(-\delta k^{-1}\beta_{k})\sim\exp(-N\beta_{k})$, excursions from the ground state to the higher band are exponentially costly and effectively confined to short intervals $\sim(N\Delta\beta)^{-1}$ in the Brillouin zone, where $\Delta\beta=2\min_{k}{\beta_{k}}$ is the purity gap. A quantitative estimate of the suppression factor may be obtained by noting that, for $k=k^{\prime}$, the integrand is real. A straightforward Taylor expansion to first order in $k-k^{\prime}$ then yields the leading-order contribution

[TABLE]

where the integral over gauge-potential components yields a nonextensive and $\beta$-independent factor. This is the result quoted in Eq. (30).

Appendix C Dynamical adiabaticity

If one chooses to measure the topological invariant $\frac{1}{2\pi}\Delta\varphi_{\text{E}}=1/(2\pi)\oint d\phi\,\partial_{\phi}\varphi_{\text{E}}$ by observing the EGP $\varphi_{\text{E}}(\phi)$ continuously in time while varying some system parameter $\phi$ along a closed loop, one must do so slowly enough to ensure that the state of the system follows the expected stationary state.

If one considers a nonequilibrium driven-dissipative system with dissipative evolution governed by a Liouvillian, the rate at which the latter reaches its stationary state after an arbitrary change of parameters is controlled by the damping gap $\Delta_{\text{d}}$ (see Sec. IV). In that case, the rate at which the EGP changes as $\phi$ is varied should be less than $\Delta_{\text{d}}$, i.e., $|\partial_{t}\varphi_{\text{E}}|\ll\Delta_{\text{d}}$. This translates as the “dynamical adiabaticity” condition

[TABLE]

which should hold all along the path $\phi\in[0,2\pi]$. Note that the “rate” of change $|\partial_{\phi}\varphi_{\text{E}}|$ generically increases when moving closer to purity-gap-closing points (recall that $\phi$ typically encircles such points). Therefore, the requirements for dynamical adiabaticity are determined by both the damping and the purity gaps.

If one considers, instead, a system at thermal equilibrium with unitary evolution governed by a Hamiltonian, the above discussion cannot be applied directly. Unless one explicitly takes into account the coupling between the system and the reservoir(s) which make(s) the latter thermalize, the damping gap is not defined and, due to the finite temperature, there is no many-body energy gap either — even if the chemical potential lies in the gap between two bands. As we have shown in the main text, however, the temperature-induced population of single-particle energy states above the chemical potential does not affect the EGP winding $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$, at any finite temperature. Therefore, in the thermal case, we expect dynamical adiabaticity to be controlled by the energy gap $\Delta\epsilon$ between single-particle energy bands below and above the chemical potential, rather than by the damping gap as in Eq. (55). The fact that this is indeed correct is exemplified in Fig. 5 for the Rice-Mele model (see Sec. IV) at half filling with $8$ sites and thermal initial state (temperature $T=10$). The hopping amplitudes $t_{1},t_{2}$ and the staggered potential $\Delta$ [see Eq. (35)] are varied continuously in time so as to encircle the purity-gap-closing point $t_{2}=t_{1}$, $\Delta=0$. Namely, we parameterize $t_{1,2}=2\pm\sin(2\pi\phi/T_{\phi})/4$ and $\Delta=-\cos(2\pi\phi/T_{\phi})/2$, and vary $\phi$ linearly in time from [math] to $2\pi$ over the time period $T_{\phi}$, i.e., $\phi=2\pi t/T_{\phi}$. Figure 5 shows the EGP as function of $\phi(t)/(2\pi)$ for different values of $T_{\phi}$. As expected, the EGP difference $\frac{1}{2\pi}\Delta\varphi_{\text{E}}$ accumulated over one cycle is approximately quantized provided that $T_{\phi}$ is large as compared to the inverse energy gap $\Delta\epsilon=\mathcal{O}(1)$ (inset of Fig. 5).