One-particle-density-matrix occupation spectrum of many-body localized states after a global quench

TL;DR

This paper investigates the occupation spectrum of the one-particle density matrix in many-body localized states after a global quench, revealing a nonthermal occupation profile and partial quasiparticle occupation influenced by initial conditions.

Contribution

It demonstrates the absence of a Fermi-liquid-like discontinuity in the steady state occupation spectrum post-quench, linking it to the local structure of the initial density-wave state.

Findings

Discontinuity in occupation spectrum is absent after a global quench.

Occupation function remains strongly nonthermal despite Fermi-liquid-like features.

Partial quasiparticle occupation can be tuned by initial state and described by an effective temperature.

Abstract

The emergent integrability of the many-body localized phase is naturally understood in terms of localized quasiparticles. As a result, the occupations of the one-particle density matrix in eigenstates show a Fermi-liquid-like discontinuity. Here we show that in the steady state reached at long times after a global quench from a perfect density-wave state, this occupation discontinuity is absent, reminiscent of a Fermi liquid at a finite temperature, while the full occupation function remains strongly nonthermal. We discuss how one can understand this as a consequence of the local structure of the density-wave state and the resulting partial occupation of quasiparticles. This partial occupation can be controlled by tuning the initial state and can be described by an effective temperature.

Equations

$$\displaystyle H=J\sum_{i=1}^{L}\Big{[}$$

$$\displaystyle+V\Big{(}n_{i}-\frac{1}{2}\Big{)}\Big{(}n_{i+1}-\frac{1}{2}\big{)}\Big{]},$$

$$|{\Psi }_0⟩=\prod _{i=1}^{L/2}{c}_{2i}^{†}|0⟩,$$

$${\rho }_{ij}(t)=⟨\Psi (t)|c_i^{†}{c}_j|\Psi (t)⟩$$

$${\rho }_{ij}(t)=\sum _{n,m}{e}^{-i(E_n-E_m)t}{a}_m^{\ast }{a}_n⟨m|c_i^{†}{c}_j|n⟩.$$

$$\overline{\rho }$$

$${\bar{\rho }}_{ij}$$

$${\overline{\rho }}_{ij}\approx \sum _n|a_n{|}^2⟨n|c_i^{†}{c}_j|n⟩\equiv {\rho }_{ij}^{\mathrm{diag}}.$$

$$⟨\Delta n^{\mathrm{diag}}⟩\equiv ⟨n_N^{\mathrm{diag}}⟩-⟨n_{N+1}^{\mathrm{diag}}⟩\approx ⟨|a_0{|}^2\Delta n^{(0)}⟩.$$

$$|{\Psi }_{\gamma }⟩=\prod _{i=1}^{L/2}[\cos (\gamma )c_{2i}^{†}+\sin (\gamma )c_{2i+1}^{†}]|0⟩.$$

$${\mathcal{I}}_{\text{OPDM}}(t)=\frac{⟨N_{+}(t)⟩-⟨N_{-}(t)⟩}{N},$$

$$n_{\alpha }^{\mathrm{diag}}=|a_0{|}^2{n}_{\alpha }^{(0)}+\sum _{n>0}|a_n{|}^2{n}_{\alpha }^{(n)}.$$

Full text

Introduction.—One of the basic notions of condensed matter physics is adiabatic continuity Anderson (1984). The Fermi liquid, a major example, is adiabatically connected to the Fermi gas by slowly turning on interactions: the ground state of the Fermi gas evolves into that of the Fermi liquid and low-energy excited states evolve into excited quasiparticle states with identical quantum numbers Nozieres and Pines (1999); Coleman (2015). The quasiparticle density operators represent conserved quantities; the Fermi-liquid Hamiltonian is diagonal in the quasiparticle basis but is not quadratic as it contains quasiparticle density-density interactions that represent the back-action of all the other particles on excitations. The Fermi liquid fundamentally relies on reduced scattering due to limited phase space, resulting from the Fermi-sphere structure of the ground state, and is therefore only a valid description at low temperatures compared with the Fermi temperature Nozieres and Pines (1999); Coleman (2015).

Closed interacting disordered quantum systems can exhibit many-body localization (MBL) Basko et al. (2006); Gornyi et al. (2005), resulting in an ideal insulator with vanishing conductivities at finite, or even infinite Oganesyan and Huse (2007); Pal and Huse (2010), temperatures. This MBL insulator is adiabatically connected to the Anderson insulator and therefore shares many features with the Fermi liquid Bera et al. (2015, 2017). In contrast to the Fermi liquid, where only the ground state and the lowly excited states are adiabatically connected to the free Fermi gas, in an MBL insulator every eigenstate is adiabatically connected to some eigenstate of the Anderson insulator—for example in fully many-body-localized systems Huse et al. (2014) (with exceptions Parameswaran and Gopalakrishnan ) where the relation can be provided by a finite-depth quantum circuit Bauer and Nayak (2013). The MBL phase is thus an emergent integrable phase Serbyn et al. (2013a); Huse et al. (2014); Vasseur and Moore (2016) characterized by conserved quasiparticle densities, which are the density operators of Anderson orbitals locally dressed by particle-hole excitations Ros et al. (2015); Imbrie (2016a, b); Rademaker and Ortuño (2016). The eigenstates are product states of these quasiparticles and therefore satisfy an area law of entanglement Bauer and Nayak (2013); Kjäll et al. (2014); Luitz et al. (2015), necessarily violating the eigenstate thermalization hypothesis Deutsch (1991); Srednicki (1994); Rigol et al. (2008). The construction of the conserved quantities (commonly referred to as $l$-bits) is extensively studied Chandran et al. (2015); Kim et al. ; Rademaker and Ortuño (2016); Monthus (2016); He and Lu ; Inglis and Pollet (2016); O’Brien et al. (2016); Pekker and Clark (2017); Pekker et al. ; Ros et al. (2015); Imbrie (2016a, b). As in the Fermi liquid, the MBL Hamiltonian is diagonal in the quasiparticle basis, but contains quasiparticle density-density interaction terms, which are absent in the Anderson insulator Serbyn et al. (2013a); Huse et al. (2014). These interaction terms give rise to dephasing in dynamics that results in a logarithmic growth of entanglement entropy Žnidarič et al. (2008); Bardarson et al. (2012); Serbyn et al. (2013b) (for examples of other quantum information measures, see Goold et al. (2015); Devakul and Singh (2015); Geraedts et al. (2016); Singh et al. (2016); Bera and Lakshminarayan (2016); Iemini et al. (2016); De Tomasi et al. (2017); Campbell et al. ), and a slow relaxation of observables towards nonergodic stationary states at long times Serbyn et al. (2014). The adiabatic connectivity of the MBL phase to the Anderson insulator relies on the stability of Anderson localization against interactions Basko et al. (2006); Gornyi et al. (2005); Ros et al. (2015); Imbrie (2016a, b).

This formal analogy between MBL and Fermi liquids was further developed in Refs. Bera et al., 2015, 2017, which evinced a Fermi-liquid-like discontinuity in the eigenvalues of the one-particle density matrix (OPDM) in many-body eigenstates, analogous to a finite quasiparticle weight (see also Ref. Rademaker et al., 2017). The discontinuity signals Fock-space localization, while the eigenvectors of the OPDM give localized orbitals, the natural orbitals, that can be used to construct an optimized single-particle approximation to the quasiparticles Bera et al. (2017).

In a Fermi liquid the occupation spectrum is discontinuous only at zero temperature; any nonzero temperature leads to a smooth occupation spectrum. With the MBL eigenstates providing an analog to a zero-temperature Fermi liquid, it is natural to ask if there is also a finite-temperature analog. We limit our consideration to temperature effects on quasiparticle occupations and assume that quasiparticle lifetimes are not affected. In this phenomenological analogy, in which each MBL eigenstate is a zero-temperature reference state, this requires partial occupations of quasiparticles compared with the reference occupations in a given eigenstate. A generic combination of eigenstates, described by a mixed density matrix, does not work as this corresponds to summing over different random occupations of quasiparticles, or to mixing reference states. Instead, we propose that a global quench from a local product state provides the physics we are after. Intuitively, a local density has a large overlap with some quasiparticle density. An expansion of such a local density in the quasiparticles will therefore mainly contain the quasiparticles localized close-by, as if they were excited by a relatively small effective temperature. Initially we focus on a perfect density-wave state as the initial state. Such a product state still has systematic phase differences between different quasiparticles unlike in thermal states. During time evolution, however, this quasiparticle superposition dephases such that the initial phase relationship is scrambled in the infinite-time steady state. The main result of our work is a characterization of this steady state with one-particle density matrix occupations that indeed mimic occupation effects of temperature in a Fermi liquid. We further show how the partial occupations in the steady state can be controlled by the structure of the initial state.

An initial density-wave state is also used in the ultracold atoms experiments that observed a finite imbalance between the density on even and odd sites as a signature of the absence of thermalization Schreiber et al. (2015); Bordia et al. (2016, 2017) (see Smith et al. (2016); Choi et al. (2016); Lüschen et al. (2017, ) for further experiments). As a corollary result we therefore obtain a relation between the OPDM occupations and experiments. In particular, we introduce an OPDM occupation imbalance, which behaves similar to the density imbalance but with a slower relaxation towards the steady state, thereby capturing dephasing.

Model and methods.—We study a system of spinless fermions hopping and repulsively interacting with their nearest neighbours in a disordered 1D lattice, with Hamiltonian

[TABLE]

where $c_{i}^{\dagger}$ creates a fermion on site $i$ (among $L$ sites) and $n_{i}=c_{i}^{\dagger}c_{i}$ is the number operator. Energies are expressed in terms of the hopping constant $J$, whereas disorder and interaction strengths are denoted by the dimensionless quantities $W$ and $V$, respectively. The disorder is diagonal and taken from a box distribution $\epsilon_{i}\in[-W,W]$. We set $J=V=1$ throughout this work, in which case the localization-delocalization transition is found to be at $W_{c}=3.5\pm 1$ for energies in the middle of the spectrum Pal and Huse (2010); Bar Lev et al. (2015); Luitz et al. (2015); Bera et al. (2015); De Luca and Scardicchio (2013).

Using exact diagonalization, we study the system described in (1) for different system sizes $L$ and average over $10^{4}$ ($L=8,10,12$), $5\times 10^{3}$ ($L=14$) and $4\times 10^{3}$ ($L=16,18$) disorder realizations. We use periodic boundary conditions and fix the number of particles to half filling $N=L/2$. The symbol $\langle\cdot\rangle$ denotes the disorder average.

The initial state, unless stated otherwise, is a perfect density-wave state,

[TABLE]

which then evolves under the Hamiltonian (1) according to (we set $\hbar=1$) $|\Psi(t)\rangle=\exp(-iHt)|\Psi_{0}\rangle.$ To characterize the state $|\Psi(t)\rangle$, we calculate the instantaneous OPDM

[TABLE]

and diagonalize it. The eigenvalues $\{n_{\alpha}(t)\}$, with $\alpha=1,2,\dots,L$, are the occupations and the eigenfunctions $\{|\phi_{\alpha}\rangle\}$ are the natural orbitals. For each time, we order the OPDM eigenvalues in descending order $n_{1}(t)\geq n_{2}(t)\geq\dots\geq n_{L}(t)$, noting that the total particle number is conserved $\sum_{\alpha=1}^{L}n_{\alpha}(t)=\mathrm{tr}\,\rho(t)=N$ at all times.

Evolution of occupations.—We first address the nature of the relaxation dynamics of the occupations $\{n_{\alpha}(t)\}$. In the initial state $|\Psi_{0}\rangle$, half of the occupations are equal to one and the other half equal to zero, i.e., $n_{\alpha}(0)=1$ for $\alpha\leq N$ and $n_{\alpha}(0)=0$ for $\alpha\geq N+1$.

The time evolution of the occupation spectrum in the MBL phase is plotted in the main panel of Fig. 1. Initially, the spectrum captures a fast expansion up to the localization length, followed by a slow relaxation in which the occupations approach their saturation values as a power law $\nu\,t^{-\gamma}+\delta$, starting at times of the order of $t\sim 10^{2}$ (see Fig. 1(c)). The parameters $\delta$, $\nu$, and $\gamma$ depend non-universally on $\alpha$ with the exponent $\gamma$ ranging between $0.3$ and $0.6$. In Fig. 1(b), the time evolution of the occupation $\langle n_{N}(t)\rangle$ is shown for both phases and at the transition ($W\approx W_{c}$). In the MBL phase ($W>W_{c}$), it undergoes a slow relaxation towards a nonthermal stationary state at long times ($t\sim 10^{8}$). This slow relaxation is due to dephasing and is characteristic of the MBL phase Serbyn et al. (2014). The instantaneous natural orbitals evolve from the initial onsite densities towards localized orbitals at long times, and the instantaneous occupations $\{n_{\alpha}(t)\}$ can therefore be seen as expectation values of local observables. In this sense, their approach to their stationary values is consistent with general arguments for power-law relaxation of local observables Serbyn et al. (2014). In the ergodic phase ($W<W_{c}$), in contrast, we observe a fast relaxation towards a thermal stationary state.

Steady-state properties.—On the basis of the above, it is natural to ask about the behavior of the occupation spectrum in the steady-state limit. To this end, we explore the asymptotic behavior ($t\rightarrow\infty$) of the time-averaged density matrix, which we compare to density matrices that capture the separate effects of dephasing and thermalization. The steady-state density matrix at long times is described by the diagonal ensemble in both the MBL and the ergodic phase; only the latter is additionally reproduced by a thermal ensemble. Specifically, expanding $|\Psi_{0}\rangle=\sum_{n}a_{n}|n\rangle$ in terms of the many-body eigenstates $H|n\rangle=E_{n}|n\rangle$, the time-evolved state takes the form $|\Psi(t)\rangle=\sum_{n}e^{-iE_{n}t}a_{n}|n\rangle$, and the density matrix is

[TABLE]

Taking the time average (denoted by $\bar{\cdot}$) yields

[TABLE]

We take $t_{0}=10^{5}$ at which point the time evolution has reached a steady state. For a nondegenerate system, the phases in (6) are random and sum to zero if $n\neq m$; therefore

[TABLE]

It is important to contrast the time-averaged density matrix with the instantaneous occupations in Fig. 1. The ordering of eigenvalues does not generally commute with time averaging, in particular if there is an interchange of occupations in the time evolution. This can be expected to occur in the MBL phase, where the eigenvalues correspond to local quantities separated in space and consequently do not couple. From now on, we therefore first time average the density matrix as in (5) and only then determine and order its eigenvalues, denoted by $\bar{n}_{\alpha}$, in descending order. The occupations $\langle\bar{n}_{\alpha}\rangle$ are plotted in Fig. 2 as a function of $L$ and for three different values of $W$. In particular, we compare $\langle\bar{n}_{\alpha}\rangle$ with the ordered eigenvalues obtained directly from the diagonal ensemble (7), denoted by $\langle{n}_{\alpha}^{\mathrm{diag}}\rangle$; both are plotted in Fig. 2(a) as a function of $L$ for $W=8$, with excellent agreement. In the ergodic phase, we further find good agreement with the eigenvalues of the thermal OPDM, $\langle{n}_{\alpha}^{\mathrm{TH}}\rangle$, obtained from $\rho_{\text{TH}}=\mathrm{tr}\left(\rho_{c}c_{i}^{\dagger}c_{j}\right)$, where we use the density matrix of the canonical ensemble $\rho_{c}=e^{-\beta H}/\mathrm{tr}\,(e^{-\beta H})$, with inverse temperature $\beta$ set by the requirement that the energy of the state be $E=\langle\Psi_{0}|H|\Psi_{0}\rangle=\mathrm{tr}\left(\rho_{c}H\right)$ (see Fig. 2(b)). The occupations obtained in the MBL phase are, in contrast, highly nonthermal as revealed by the comparison with the thermal occupation spectrum, plotted as stars in Fig. 2(a). The OPDM occupations tend to exhaust the full range of values between 0 and 1, similar to the occupations in eigenstates Bera et al. (2015, 2017), but with a discontinuity that vanishes exponentially as $L\to\infty$, see the inset in Fig. 2(a). This main result of our work suggests that a global quantum quench from a product state of local densities results in partial quasiparticle occupations and thus a continuous occupation spectrum, similar to the effect of a finite temperature in a Fermi liquid.

The absence of the discontinuity is best understood in the diagonal ensemble. With the initial state being a product state of single-site occupations, it has a large overlap with the eigenstates that have a large weight on these sites. The quenched state will, to first order in perturbation theory, inherit the step function from this eigenstate, while all the other states provide smearing of the step. To support this argument, we order the many-body eigenstates $|n\rangle$ according to their overlap with the initial state such that $|a_{0}|\geq|a_{1}|\geq\cdots\geq|a_{n}|$, and in Fig. 2(c), we plot the disorder-averaged overlap $\langle|a_{n}|^{2}\rangle$ as a function of $n$. This function decays quickly with $n$ and the largest overlap also decays exponentially with $L$. This is consistent with eigenstates built from $N$ quasiparticles each with an overlap with absolute value $c<1$ with a given site density, and therefore, a total overlap that scales like $|a_{0}|\propto c^{N}$. The maximum-overlap eigenstate $|0\rangle$ has an OPDM $\rho^{(0)}$ with a zero-temperature Fermi-liquid-like step function. The unitary transformation that diagonalizes $\rho^{(0)}$ approximately diagonalizes the OPDM $\rho^{(n)}$ of the higher eigenstates $|n\rangle$, but with random ordering such that the disorder average $\langle n_{\alpha}^{(n)}\rangle$ becomes a smooth function without any discontinuity (see Appendix A) for detailed calculations supporting this picture). The resulting prediction of the diagonal ensemble for the discontinuity is then

[TABLE]

This indeed goes to zero exponentially with $L$ since $\langle|a_{0}|^{2}\rangle\sim e^{-bL}$, consistent with the inset to Fig. 2(a).

Tuning quasiparticle occupations.—We can systematically tune the distance to a reference MBL eigenstate by generalizing the density-wave inital state to

[TABLE]

In Fig. 3 we plot the steady-state occupation spectrum for several values of $\gamma$. With increasing $\gamma$ the energy variance of the state $|\Psi_{\gamma}\rangle$ increases, corresponding to exciting quasiparticles and an increased effective temperature, which is reflected in the long-time occupation spectrum deviating more and more from the step function. This provides a systematic, and experimentally feasible Folling et al. (2007), procedure to tune the quasiparticle occupations.

While the energy variance provides a proxy for the effective temperature, by being a measure of the amount of quasiparticle excitations, it is important that these excitations are with respect to a definite reference eigenstate that acts as the ground state. This is guaranteed by the local structure of the initial state $|\Psi_{\gamma}\rangle$. If we instead take an initial state without such a reference, for example a highly excited eigenstate $|\Psi_{\mathrm{free}}\rangle$ of the clean Hamiltonian (1) with $W=0$, we obtain a practically flat distribution $\langle n_{\alpha}^{\mathrm{diag}}\rangle\approx 1/2\,\,\forall\alpha$), see diamonds in Fig. 3. This is the case even though the energy variance of this state is similar to that of $\Psi_{\gamma}$ with $\gamma=\pi/4$, and results from the mixing of reference states (see Appendix B).

Imbalance and connection to experiments.–The density imbalance $\mathcal{I}=(N_{\text{e}}-N_{\text{o}})/N$ between the number of physical particles $N_{\text{e}}$ on even sites and $N_{\text{o}}$ on odd sites is experimentally seen to relax to zero in the ergodic phase, whereas in the many-body localized phase it exhibits a fast relaxation towards a nonzero value, reflecting the absence of thermalization Schreiber et al. (2015); similar conclusions were obtained numerically in Ref. Luitz et al., 2016. From the occupation spectrum, we define a related imbalance between the occupied and unoccupied halves of the spectrum as

[TABLE]

where N_{+}(t)=\sum_{\alpha=1}^{N}\big{(}n_{\alpha}(t)-n_{\alpha}^{\mathrm{TH}}\big{)} and N_{-}(t)=\sum_{\alpha=N+1}^{L}\big{(}n_{\alpha}(t)-n_{\alpha}^{\mathrm{TH}}\big{)}. We can view $\mathcal{I}_{\mathrm{OPDM}}(t)$ as a measure of how close a state is to a step-function occupation spectrum where the imbalance is maximal, or to a thermal occupation spectrum where its value is zero. This imbalance is plotted for both phases in Fig. 4(a); it saturates to zero in the ergodic phase (dotted lines) and to a finite value in the MBL phase (solid lines). This conclusion is supported by the finite-size scaling given in Fig. 4(c). In comparison with the density imbalance $\mathcal{I}(t)$ in the MBL phase, plotted in Fig. 4(b), the relaxation of the OPDM imbalance is much slower. The reason for this is that $\mathcal{I}(t)$ only captures the ballistic expansion part of the relaxation, while $\mathcal{I}_{\mathrm{OPDM}}(t)$ also captures the dephasing mechanism originating from interactions between quasiparticles. As a second main result of our work, we have thus demonstrated that the slow relaxation of the OPDM occupation spectrum can be directly connected to lack of ergodicity and experimentally accessible quantities.

Conclusions.—We demonstrated that in the many-body localized phase, the steady state reached at long times after a quench from a local product state has a smooth OPDM occupation spectrum with a highly nonthermal shape, in contrast to that obtained in the ergodic phase. This is consistent with the picture of local conserved quantities which have a significant overlap with the initial state. The approach towards the steady state is consequently a power law, reflecting dephasing via interactions between quasiparticles. We have also defined an occupation imbalance, similar to the density imbalance used in experiments, that captures the main effect of dephasing and absence of thermalization.

Our discussion suggests that the continuous occupation spectrum is phenomenologically similar to that of a finite-temperature Fermi liquid. The finite temperature is provided by the energy difference between the initial state and the closest eigenstate, which serve as reference states with a (zero-temperature) Fermi-liquid-like occupation spectrum. This effective temperature can be tuned by changing the local structure of the initial state. This is not thermalization in the conventional sense since the many-body localized phase is manifestly nonergodic. Nevertheless, the observation that the steady state OPDM spectrum is continuous hints at the possibility of describing it with an emergent temperature, capturing the energy variance in the system. It remains an interesting future research direction to establish whether such an emergent temperature corresponds to some thermal-like ensemble, necessarily different from eigenstate thermalization, and then how one can characterize it.

Appendix A Time average of the OPDM and smearing of its occupation spectrum in the MBL phase.

In this appendix we provide more details on some technical points mentioned in the main text, starting with the time average of the density matrix. First, in Fig. 5 we plot the instantaneous OPDM spectrum at large times ($t=10^{8}$), obtained by first ordering the occupations by size and then taking the disorder average. This should be compared with the time-averaged OPDM spectrum in Fig. 2 of the main text. We notice that the instantaneous occupations are systematically larger, although as the system size becomes larger, the distribution in the thermal phase approaches its infinite-time average, whereas the distribution in the MBL phase does not. As already mentioned in the main text, this has mainly to do with the fact that the ordering of eigenvalues does not commute with the disorder average when occupations are interchanged as a function of time.

Second, we give more details on the arguments leading to Eq. (8) in the main text. The maximum-overlap eigenstate $|0\rangle$ has an OPDM $\rho^{(0)}$ with a zero-temperature Fermi-liquid-like step function, as we confirm in Fig. 6, where we plot its eigenvalues $n_{\alpha}^{(0)}$ (the same is the case for the other eigenstates, see lines). The unitary transformation $U_{0}$ that diagonalizes $\rho^{(0)}$ approximately diagonalizes the OPDM $\rho^{(n)}$ of the higher eigenstates $|n\rangle$, such that the diagonal elements $n_{\alpha}^{(n)}=(U_{0}^{-1}\rho^{(n)}U_{0})_{\alpha\alpha}$, if ordered in decreasing order before taking the disorder average, also have an OPDM occupation discontinuity (see filled symbols in Fig. 6). The transformation $U_{0}$, however, randomizes the ordering. To demonstrate this, the occupations in the zeroth, first and second state $\{n_{\alpha}^{(n)}:n=0,1,2\}$ are plotted for three disorder configurations in Fig. 7. As a consequence, the disorder average $\langle n_{\alpha}^{(n)}\rangle$, plotted with unfilled symbols in Fig. 6, becomes a smooth function without any discontinuity.

Assuming then that the diagonal elements of the diagonal density matrix are obtained by applying $U_{0}$ term by term in Eq. (6), we write

[TABLE]

Here, as in the main text, the $\{|n\rangle\}$ are ordered in order of decreasing $|a_{n}|$. Given that all the higher order terms in Eq. (S1) are smooth, even at finite system sizes, they do not contribute to the discontinuity $\Delta n$, which is therefore given by Eq. (8) in the main text. In Fig. 9 we compare the discontinuity $\Delta n$ obtained from the time average with that obtained from the approximation of Eq. (8). The agreement between the two quantities is rather good, given the approximations used in the argument leading up to the prediction of Eq. (8), and corroborates the conclusion that in the thermodynamic limit, the discontinuity indeed goes to zero.

Appendix B Reference initial states and their overlap with the many-body eigenstates.

In the main text, we stated the importance of selecting initial states that resemble zero-temperature reference states (many-body eigenstates). We also provided a way to produce such a class of initial states $|\Psi_{\gamma}\rangle$, and compared them with an initial state $|\Psi_{\mathrm{free}}\rangle$ that is not representative of such reference states. Here, we show in Fig. 9, how the local structure of the initial states is reflected in their overlap $\langle|a_{n}|^{2}\rangle$ with the many-body eigenstates $|n\rangle$ (in line with the analysis in Fig. 2(c) in the main text). As we increase $\gamma$, the largest overlap $\langle|a_{0}|^{2}\rangle$ between the initial states $|\Psi_{\gamma}\rangle$ and the eigenstates decreases, but even for large $\gamma$, there is an eigenstate with a significant weight that serves as a reference. In contrast, the overlaps of $|\Psi_{\mathrm{free}}\rangle$ with $|n\rangle$ are essentially constant.