Probing infinite many-body quantum systems with finite-size quantum simulators

Experimental studies of synthetic quantum matter are necessarily restricted to approximate ground states prepared on finite-size quantum simulators. In general, this limits their reliability for strongly correlated systems, for instance, in the vicinity of a quantum phase transition (QPT). Here, we propose a protocol that makes optimal use of a given finite-size simulator by directly preparing, on its bulk region, a mixed state representing the reduced density operator of the translation-invariant infinite-sized system of interest. This protocol is based on coherent evolution with a local deformation of the system Hamiltonian. For systems of free fermions in one and two spatial dimensions, we illustrate and explain the underlying physics, which consists of quasi-particle transport towards the system's boundaries while retaining the bulk"vacuum". For the example of a non-integrable extended Su-Schrieffer-Heeger model, we demonstrate that our protocol enables a more accurate study of QPTs. In addition, we demonstrate the protocol for an interacting spinful Fermi-Hubbard model with doping for 1D chains and a small two-leg ladder, where the initial state is a random superposition of energetically low-lying states.


I. INTRODUCTION
A major goal of quantum simulation is the study of ground states (GSs) of quantum many-body systems with phase transitions emerging in the thermodynamic limit. Aiming for accurate characterization of GS phase diagrams raises the question of how to directly probe -on a finite-size quantum simulator -the GS |Ω ∞ of an infinite system described by a translation-invariant Hamiltonian H ∞ , especially in the vicinity of and at critical points [see Fig. 1(a)]. In practice, any quantum simulator will necessarily be limited to a finite size, say N d sites for a d-dimensional lattice system. In this context, we note the remarkable recent experimental progress [1,2] culminating in the realization of almost perfectly isolated synthetic quantum systems consisting of tens or hundreds of atoms as Hubbard or spin lattice models in various spatial dimensions [3][4][5][6][7][8][9][10][11][12][13][14][15][16]. The restriction to finite size is less relevant as long as the intrinsic correlation length ξ of the target state |Ω ∞ is sufficiently small compared to the available system size, i.e., ξ N (with ξ in units of the lattice spacing). In this case, an approximation of the finite-size GS |Ω N of H N [system S, see Fig. 1(b)], which might have been prepared via adiabatic state preparation or variational techniques [4,7], allows for a faithful study of |Ω ∞ . However, if the correlation length diverges, ξ → ∞, such as at a quantum phase transition (QPT), significant finite-size effects destroy a faithful representation of the infinite-size GS on a finite system with N ≤ ξ.
In this paper we resolve this situation by considering pure many-body states |ψ , living on a finite-size quantum simulator S of linear dimension N , such that its reduced density operator ρ n for a subsystem A of size n < N equals the mixed state represented by the reduced (1) Here the traces Tr ¬A are taken over the complement of A in the finite and infinite systems, respectively. On a formal level, we note that quantum information theoretical considerations guarantee the existence of a purification |ψ of ρ n for the case N = 2n [17], with n lattice sites representing the system of interest A, and the remaining n sites playing the role of an auxiliary system, or 'reservoir'. The challenge to be addressed in a quantum simulation setting is to find specific, experimentally realizable protocols which allow the preparation of the purifications |ψ for properly chosen n for given N . Below we will refer to such protocols as Purification Preparation (PP).
The central result of this paper is that Purification Preparation can be achieved in analog quantum simulation setups as engineered quench dynamics with spatially deformed system Hamiltonians, something which is readily implemented in today's experiments. We consider a situation where the analog quantum simulator provides us not only with an implementation of the (in bulk translationally invariant) Hamiltonian of the form H N = j∈S h j , but also allows the realization of spatially deformed Hamiltonians H (def.) N = j∈S g j h j with {g j } a given spatial pattern of local couplings [18]. In its simplest form, the protocol starts with the preparation of a pure state approximating the GS of the finite system, |ψ(0) ≈ |Ω N , with H N |Ω N = E |Ω N . Our claim is that the coherently evolved state |ψ(t) = e −iH (def.) N t |ψ(0) provides an approximate purification |ψ for a proper choice of deformation parameters {g j } (in particular in form of a parabolic deformation) and evolution time t = t PP , i.e. |ψ = |ψ(t PP ) . , parametrized by gj, on the quantum simulator (obtained in the continuum limit of freefermion system for evolution with a parabolically deformed Hamiltonian, see Appendix A). At time tPP the initial state |ΩN is transformed into a purification |ψ . (d) The final state |ψ with subsystem A that corresponds to the part of |Ω∞ . White to red colors indicate the density of excitations (colorbar at the bottom).
While the present paper will provide detailed reasoning behind the functioning of the protocol from various perspectives, which includes connections with recent studies of quench dynamics within conformal field theory (CFT) [19][20][21][22][23] and numerous numerical illustrations, we outline in Fig. 1 a simple physical picture behind PP as engineered quasi-particle dynamics. Here we interpret the initial state |ψ(0) ≈ |Ω N in terms of excitations above the target state, and the physical essence of the dynamical evolution with the parabolically deformed Hamiltonian consists of "cleaning" the bulk A from these excitations by moving them to, and accumulating them at the edge S A [ Fig. 1(c)], and thus effectively "cooling" the bulk. In this way the PP protocol essentially eliminates boundary effects and restores translation invariance in the bulk A. We demonstrate our protocol for non-interacting fermions in one and two spatial dimensions as well as an interacting spin chain, a 1D Fermi-Hubbard model and preliminary results for a Fermi-Hubbard model on a two-leg ladder. These examples suggest broad applicability of our approach to both interacting and non-interacting systems in different spatial dimensions. While the present work focuses on initial pure states in PP, we emphasize that the approach is readily applied to mixed thermal states, where again the bulk is "cooled" by accumulating thermal excitations dominantly at the edge.
Finally, in connecting the present theoretical work to a feasible experimental protocol, tools for verification, quantifying the approach to the target state ρ n , must be developed. This can be achieved via recently developed techniques for Entanglement Hamiltonian (EH) learning [18,[24][25][26][27][28]. Writing ρ n ∼ exp(−H n ) withH n defining EH, the Bisognano-Wichmann theorem of CFT [29][30][31][32] applied to lattice models (see, e.g., [18,25,33]) predicts that the EH for ground states has the simple structure of a parabolically deformed system Hamiltonian. EH learning thus provides the techniques to monitor and thus verify the approach to the target state ρ n in PP. Again we will illustrate this with examples.
The structure of this paper is as follows. We first provide an overview of our results in Sec. II where we introduce our state preparation protocol and present the results of its application to an interacting spin chain, demonstrating the potential of our approach for mapping out the GS phase diagram. In the following Sec. III, we discuss the physics underlying our protocol in more detail. For the examples of non-interacting fermions in one and two spatial dimensions, we provide extensive numerical evidence supplemented by analytical predictions to explain how the desired purification is realized. Subsequently, we return to the interacting spin chain in Sec. IV and demonstrate that our protocol enables a more accurate determination of universal critical exponents that characterize a QPT and the underlying CFT. In Sec. V we demonstrate the protocol for an interacting spinful Fermi-Hubbard model with doping for 1D chains, and further results for two-leg ladders are presented in Appendix E. In this study the initial state are taken as random superposition of energetically low-lying states, i.e. effectively mixed initial states. We conclude our work with a brief discussion in Sec. VI. The appendices contain additional details about our protocol, its analysis, and our numerical and analytical calculations.

II. OVERVIEW OF CONCEPTS AND RESULTS
This section will provide an overview of the concepts and application of the Purification Preparation protocol. We will first describe and motivate the protocol per se, followed by the illustrative example of the study of the phase diagram of the extended Su-Schrieffer-Heeger model, comparing the case of finite vs. infinite system.

A. Purification Preparation Protocol
We consider a quantum simulator with finite linear size N , and we assume here pure initial states |ψ(0) , which are close to the GS |Ω N of a translation-invariant finitesize local Hamiltonian H N . In a typical quantum simulation experiment |Ω N corresponds to a GS with open boundary conditions (OBC). For simplicity of writing, we describe here the case of a finite one-dimensional (1D) chain with OBC (see Fig. 1 and h j,j+1 couple neighboring lattice sites in the total system S with N sites. Our goal is to realize the reduced state ρ n [Eq. (1)] on n < N sites in the middle region A.
Our protocol to prepare the purification |ψ in Eq. (1) as |ψ(t PP ) consists of two main steps: 1. Evolution of the initial state |ψ(0) with a deformed Hamiltonian H N t |ψ(0) , for some time t. As we argue below, a good choice for models with dynamical critical exponent z = 1 is corresponding to a parabolic deformation of the system Hamiltonian 2. Finding the optimal time t PP to stop the evolution as the time that corresponds to the first maximum of the second Rényi entropy (purity) of |ψ(t) measured for the edge region E ⊂ S A of the simulator S.
Let us briefly comment on the rationale behind these steps (for more details, see Sec. III). The main idea relies on interpreting the initial state |ψ(0) as a subsystem of the infinite system in some excited state. More formally, this can be formulated as follows: According to Kraus' theorem [34], the initial state can be written as where the operatorsÔ α act non-trivially only on the finite system S, viewed as a subsystem on an infinite system, creating there excitations above the ground state |Ω ∞ . The first step of our protocol implements an evolution to in such a way that the support of the evolved operatorŝ in the subsystem A ⊂ S decreases with time t [35]. This describes transport of the excitations from the bulk towards the edge of the simulator (see Sec. III B for an explicit demonstration of this transport). Ideally, we want that at time t PPÔα (t) has no more effect on A, such that If additionally Tr i.e., if the "bulk vacuum" is preserved by the evolution, we obtain the desired result Tr This preservation is guaranteed if the deformed Hamiltonian is proportional to the Entanglement Hamiltonian (EH)H In this context, the Bisognano-Wichmann theorem [29,30] for a CFT [31,32], which states that the EH is given by a parabolic deformation, suggests our choice of Eq. (2).
The criterion used in the second step to determine the optimal time t PP is motivated by the spreading of entanglement due to the propagation of entangled quasiparticle pairs [36]. We expect the (von Neumann) entanglement entropy S S A of the boundary region S A, which equals the entropy S A of the bulk for an initial pure state, to grow while excitations are transported out of A. Eventually, these excitations will reflect at the boundary of S and return into the bulk A. We use this to estimate the optimal evolution time t PP as the first maximum of the entropy of an edge region E ⊂ S A, which for a spatially symmetric system we take as a few sites on one side of the chain, described by state ρ E . Further, we choose to measure the second Rényi entropy as an experimentally accessible [37] alternative for the von Neumann entropy of the edge region E. A more accurate way to monitor the build-up of the entanglement structure corresponding to the target reduced state ρ n is discussed further below in section III C. Before showing the first example of our protocol, we would like to briefly discuss its experimental requirements. Both state preparation, e.g., adiabatically or variationally, as well as quench dynamics, are nowadays routinely implemented in quantum simulation experiments. The main challenge of our protocol thus consists of realizing a deformation of the system Hamiltonian. This requires a certain degree of programmability, i.e., we explicitly assume the capability of tuning local Hamiltonian parameters. Alternatively, the corresponding time evolution can also be realized stroboscopically. For completeness, we show an explicit example of such a "Trotterized" approach in Appendix D.
Before entering a detailed technical discussion of our protocol, we find it worthwhile to present results illustrating the power and potential of our approach. We focus on obtaining a quantum phase diagram obtained from observables measured in A for the purification |ψ(t PP ) . Specifically, we study an extended Su-Schrieffer-Heeger [38,39] (SSH) model on a system S with an even number of sites N and OBC, described by the 1D Hamiltonian where {σ x j , σ y j , σ z j } denote the Pauli operators at sites j. The phase diagram of this model, which is for δ = {0, 1} is generally non-integrable, is shown in Fig. 2(a). In the GS |Ω N , the symmetry-protected topological (SPT) phase is detected by measuring highly non-local many-body correlations, such as the partial reflection invariant [39] in the bulk subsystem A of n < N sites (n is even) centered in the middle of the system. Here, ρ A (and analogously ρ A L/R ) is the reduced density matrix of subsystem A = A L ∪ A R with A L and A R the corresponding left and right partitions of A with n/2 sites. The operator R A interchanges A L and A R with respect to the reflection center. In the thermodynamic limit, sending first N → ∞ and then n → ∞ with n being a multiple of four, Z n,N approaches quantized values: Z n,N → 1 for the trivial phase, Z n,N → −1 for the SPT, and Z n,N → 0 for the symmetry-broken antiferromagnetic (AF) phase.
Since the typical subsystem size n required to achieve convergence for Z n,N is determined by the correlation length ξ [39], increasingly large n becomes necessary when approaching a phase transition. At the same time, in order to avoid boundary effects in the OBC GS |Ω N , n must be taken relatively small compared to the system size N . We, therefore, expect a significant advantage in determining the phase boundaries on a finite system when using, instead of |Ω N , the state |ψ(t PP ) prepared by our protocol.
This expectation is confirmed in Fig. 2, where we present the results for the phase diagram obtained using a numerical simulation of our protocol: the order parameter Z (PP) n,N for the state |ψ(t PP ) with |ψ(0) = |Ω N is shown in comparison to Z (OBC) n,N for the OBC GS |Ω N , and Z n,∞ for the infinite system. The finite systems are simulated for N = 36 using MPS techniques [40]: DMRG for finding |Ω N and Trotter gates for the time evolution resulting in |ψ(t PP ) . The infinite case is treated with iDMRG [40,41] for details). The dependencies of Z (PP) n,N and Z (OBC) n,N on n for δ = 1.4 and δ = 2.5 are shown in Figs. 2(b and (c). We observe that the convergence for the OBC GS [(b)] stops for n N/3 due to boundary effects, resulting in a large uncertainty in identifying the transition points. In contrast, due to suppression of the boundary effects, our protocol allows for larger values of n and provides the results [(c)] which are much closer to those of the infinite system. This enables much more precise identification of the phase transitions and leads to a quantitatively improved phase diagram with sharp and accurate phase boundaries, which are in very good agreement with the thermodynamic limit results, as demonstrated in Fig. 2(a) for n = 28. Further below in section IV, we will return to the present model and show that also critical exponents, which characterize the phase transitions, can be extracted more accurately using the PP protocol.

III. DETAILED ANALYSIS FOR NON-INTERACTING FERMIONS
We now turn to an in-depth discussion of our protocol. We start with a one-dimensional chain of non-interacting fermions at half filling described by the Hamiltonian (throughout this article, we employ units where H N is dimensionless) where c ( †) j denote fermionic annihilation(creation) operators at site j. In the thermodynamic limit (N → ∞), H N becomes a gapless model described by a low-energy effective theory that exhibits a relativistic dispersion with dynamical critical exponent z = 1, namely the CFT of a free fermion. The goal of our protocol is to recover the properties of the GS |Ω ∞ of the infinite-size Hamiltonian, which has a translation-invariant energy density h j,j+1 and two-point correlators that decay algebraically with the spatial distance |2l + 1|.
In the following, we consider a finite chain with N = 36, for which Fig. 3(a) demonstrates the spatially modulated relative deviation of h j,j+1 from the exact value in the infinite GS, during evolution with a parabolically deformed H  are strongly decreased at t PP . Similarly, the initial correlators |C j,j+2l+1 | decay too fast or too slow depending on j being odd or even [see Fig. 3(c)], while translation invariance in the bulk and the target algebraic decay of correlators is restored at t PP . Below we explain how the PP protocol achieves these results for subsystems A ⊂ S with size n < N .
A. Physical explanation of the protocol As outlined in the previous section, we regard the state |Ω N as a part of a mixture of excited infinite-size stateŝ O α |Ω ∞ . In the considered system, these excited states differ from the ground state by a number of particle and hole excitations (quasi-particles) with momenta close to the Fermi momentum, which is p F = π/2 (or 3π/2) in our units. This situation can be described by wave packets of quasi-particles with momentum p created by the operatorsc † p ∼ j=∞ j=−∞ e ijp c † j for p > p F and by the Hermitian conjugate operatorc p for p < p F , which populate |Ω ∞ . In the infinite system, the quasi-particles have a dispersion (p) = |cos p|. For a small deviation ∆ rel h j,j+1 , the quasi-particles in the OBC GS |Ω N correspond to low-energy excitations, for which p ≈ p F , and the dispersion takes a linear relativistic form v F |p − p F | with a constant group velocity v(p) = ∂ε (p) /∂p = ±v F with v F being the Fermi velocity. The idea of our protocol is to realize a unitary evolution U = e −it PP H (def.) N which moves these quasi-particles to the boundaries of S, without creating new excitations in the middle of the system. In the absence of quasi-particles, the bulk A approximates a subsystem of the true vacuum, and ideally This evolution can be achieved with a local deformation of the finite-size Hamiltonian, H with an appropriately chosen deformation g j . A sufficiently smooth deformation locally modifies the velocity, v j ∝ g j , such that the quasi-particles quickly leave the bulk and spend a long time at the edge of the finite system before going back into the bulk (see Appendix A for a detailed derivation). At the same time, g j should be chosen such that the properties of the bulk vacuum are retained. Intuitively, any deformation g j that smoothly changes from unity at the center of the bulk to zero at the system's edge should provide such dynamics (see also the studies [20,21,23] of Floquet CFTs with deformed Hamiltonians). Since the bulk vacuum is conserved for evolution with a Hamiltonian proportional to the EH (or commuting with it), we expect that a quasi-local approximation of the EH will fulfill the desired conditions and provide the desired purification.
Returning to the chain of free fermions [Eq. (7)], a parabolic deformation, g j = (N − j)j/(N/2) 2 , is a quasi-local approximation [42] of the EH Hamiltonian that also exactly commutes [43] with the reduced den- Moreover, under evolution with this parabolic deformation any low-energy quasi-particle in the bulk reaches the edge of the system exponentially fast with a time-scale t ∝ N [Eq. (A10) in Appendix A gives an analytical expression for the quasi-particle's trajectories, with several trajectories demonstrated in Fig. 1(c)]. We note that a deformation function g j that connects more smoothly to zero at the edge of the system, specifically the so-called sine-square deformation that has been studied for CFTs [19,21], leads to an algebraic rather than an exponential approach to the edge (see Appendix A).
To test our predictions, we simulated the evolution of the initial state |Ω N under the deformed Hamiltonian for N = 36, (see Appendix B for details on the numerics). The resulting motion of quasi-particles can be tracked indirectly by monitoring the system's energy density during the evolution. After a certain evolution time t PP , the initial oscillations of the energy density are strongly suppressed, and the bulk energy density exhibits translational invariance with a value close to the target one of the true vacuum, as shown in Fig. 3(a) and Fig. 3(b). At the same time, the two-point correlators |C j,j+2l+1 | in the bulk, plotted in Fig. 3(c), demonstrate the target algebraic decay.
More complex correlations in the bulk state ρ n (t PP ) = Tr also approach their target values, which follows from the decrease of the bulk infidelity with respect to the target reduced state ρ n = Tr , grows approximately as n k with k changing to a larger value at n ≈ 28 indicating a "concentration" of the quasi-particles within a relatively small edge region, ≈ 4 sites at each side. The same conclusion can be deduced from the faster growth of the deviation of the energy density at these edge sites [ Fig. 3(b)] from its target value. Since the infidelity is related to the density of the quasi-particles [44], these numerical observations strongly support our arguments about the evolution with the parabolically deformed Hamiltonian.
In an experiment, it is hardly possible to measure the infidelity with respect to the true (Fermi sea) vacuum. Instead, we estimate the time t PP to stop the evolution by measuring the second Rényi entropy [37] on a small subsystem at the system boundary. This is illustrated in Fig. 3(d), where we plotted the entropy S (2) L of the two left-most sites of the chain. The maximum of the entropy closely matches the minima of the infidelities, which we attribute to the boundary region S A "storing" the necessary entanglement required for the global pure state to produce the correct mixed state in the bulk. For any finite lattice spacing, quasi-particles eventually reflect at the system's boundaries and return into the bulk, thereby decreasing both the edge entropy as well as the bulk fidelity for times t > t PP . We also note that the edge entropy provides experimentally more favorable criterion for determining the optimal time t PP than measuring small deviations of the energy density in the bulk from its target homogeneous value, see Fig. 3(b). This is because the values of the entropy, S L ∼ 1, and its rapid decrease for t > t PP provide much better visibility in the presence of experimental shot noise. Moreover, the entropy S (2) LR corresponding to two edge sites on each end of the chain [also shown in Fig. 3(d)] indicates a similar maximum as S (2) L with deviations that result from correlation among the two boundaries of the chain. In view of experimental feasibility (in particular for higher dimensions), we focus on the simpler edge entropy S (2) L throughout the rest of this paper.
Finally, Fig. 3(f) presents the numerically obtained values of the time t PP and the corresponding infidelity 1 − F n versus the total size N for the fixed subsystem sizes n = 10 and n = 30 before and after PP protocol. As expected from the quasi-particle picture, we obtain t PP ∼ N . The infidelities in both cases decrease approximately as N −k , where we find k to be independent on the bulk size n. For the OBC GS we retrieve k (OBC) ≈ 2.0 while for the PP states k (PP) ≈ 5.4. Therefore, our protocol on a simulator of size N results in the same infidelity in the bulk as a simulator of the size ∼ N 2.7 when using the OBC GS. We can further estimate the time required to prepare the corresponding states in both strategies: Using a non-linear adiabatic sweep [45], one can prepare the OBC GS of size N in a time ∼ N . With the time t PP ∼ N required for the PP protocol, this yields the total preparation time ∼ 2N for the PP state of size N . To prepare the OBC GS which results in the same infidelity, on the other hand, one needs time ∼ N 2.7 . These estimates demonstrate the advantages of our protocol for critical system studies.
B. Interpretation of quasi-particle motion as mode mapping In the following, we give another view on the physics that underlies the PP protocol for free fermions [Eq. (7)]. Specifically, we study how the PP dynamics maps finitesize fermionic modes in momentum spacec p , to the modes in the position space c j In the Heisenberg picture, evolution with the deformed Hamiltonian H (def.) N = j jh j j c † j c j , withh described by a Hermitian N × N matrix, acts on the fermion operators in position space as Together with Eq. (11) we obtaiñ For any sufficiently smooth deformation g j , we can study the evolution in the continuum limit, N → ∞ while keeping the length of the chain constant (see Appendix A for details). We show that, in the long-time limit, the PP evolution maps all "low-energy" momentum modesc p (t), i.e, modes with p in the vicinity of the Fermi points p = π/2 and p = 3π/2, to position space modes c j localized at the system's edges, i.e., with j close to j = 1 or j = N . In particular, for the parabolic deformation [Eq. (2)], the sum m j=1 |χ p,j | 2 for some m < N and p ≈ π/2 approaches unity exponentially [see Eq. (A15)], that is, the low-energy momentum modes quickly "collapse" to the edges. Once again, we note that a similar effect occurs for the SSD, albeit with an algebraic approach [see Eq. (A16)]. The "mode mapping" provides further insight into the PP mechanism when one considers the dynamics of the correlator c † pcq starting from the initial OBC GS |Ω N . The deviation of the correlator for the finite system of the size N = 36 from the one for the infinite system, ∆ c † pcq = | c † pcq − c † pcq ∞ |, at times t = 0 and t = t PP is shown in Fig. 4(a) in momentum space and in Fig. 4(b) in position space, respectively. We see that initially the deviation comes mainly from low-energy excitations with p and q being around the Fermi momenta π/2 and 3π/2, see in Fig. 4(a). Since the parabolically deformed Hamiltonian does not perturb the true vacuum, the evolution only maps the initial deviation of the cor- relator in the momentum space to the deviation of the spatial correlators at the system boundaries, as shown in Fig. 4(b). This is illustrated in Fig. 4(c), where we show the product X p X q with X p = j∈A |χ p,j |, which quantifies the contributions of quasi-particles with initial momenta q, p to the two-point correlator. The bulk is thus "cleaned" of excitations and provides a highly improved approximation of the true vacuum.
In practice, any initial state with relatively high energy density also exhibits excitations with momenta far from the Fermi momentum. Since the equation of motion for these quasi-particles is more complicated to solve analytically in the continuum limit (see Appendix A), we stud-ied numerically the momentum-position mapping function χ p,j (t) [Eq. (14)]. The left column of Fig. 5 shows |χ p,j (t)| versus time t and sites j for fixed values of p demonstrating the evolution of the position-space distribution for the selected modes p. We find that modes p away from the Fermi momentum still collapse to a small region in position space, however, shifted towards the center of the system. We conclude that in the presence of corresponding excitations at higher energies, the efficiency of the PP degrades in the sense that the maximum size of an effectively cooled bulk decreases.
The right column of Fig. 5 shows the matrix |χ p,j (t)| at fixed times where the corresponding modesc p shown in the left columns of the same figure are maximally localized in position space. These plots illustrate the spatial distribution of all momentum modes at the selected moments of time. We observe that the optimal focusing time for the modes decreases for p further away from the Fermi momentum. Thus, for a finite amount of initially populated momentum modes in the vicinity of the Fermi momentum, one can, in principle, optimize the time, where most of these modes are approximately localised in a finite edge region. This is demonstrated, i.e., by the middle plot of the right column of Fig. 5.
As a summary, the mode mapping provides another interpretation of the PP dynamics, similar to the behaviour of light passing through an optical lens. While the above considerations have been explicitly discussed for a 1D chain, we present evidence in favor of the applicability of our PP protocol in two dimensions in Sec. III D and Appendix E. Since any reduced state ρ n ∝ exp(−H n ) is determined by the corresponding EHH n , monitoringH n during the PP evolution allows us track the approach to the target state ρ n . This becomes particularly feasible, e.g., at a critical point described by a CFT, where the Bisognano-Wichmann theorem predictsH n to be a parabolic deformation of the system Hamiltonian. To observe the build-up of this parabolic shape, we apply Hamiltonian learning [26,27] to find a local approxima-tionH n = j g (n) j h j,j+1 of the EH [18,25,28] during the PP dynamics. Fig. 6 shows the resulting EH parameters g (n) j , extracted from a numerical simulation of the PP protocol for two different subsystem sizes. We observe that the EH indeed becomes very close to the anticipated parabolic shape at an optimal time t PP . This provides an experimentally accessible verification of the preparation of the desired target state with PP, and provides a way to identify the optimal time t PP in an experiment (see also below). For a discussion of Hamiltonian and EH learning protocols in an experimental context, we refer to [18,24,25]). We emphasize that the results shown in Fig. 6 correspond to a thermal initial state ∝ e −H N /T at temperature T = 0.15. This is in contrast to the pure initial states |ψ(0) considered so far. Since we expect that any low-energy state evolves towards a purification of ρ n , a mixture of low-energy states, such as a low-temperature thermal state, should also evolve to a mixture of purifi-cations with the correct reduced density matrix. Our simulations support these expectations and demonstrate that the purity of the initial state is indeed not essential. Thus the present scheme achieves, with coherent quench dynamics, a 'cooling of the bulk', i.e., a preparation of the reduced density matrix ρ n associated with the ground state of the infinite system.
We further note that the minimum of the infidelity (shown on the side of Fig. 6) perfectly agrees with the time at which the reconstructed EH is closest to the parabolic shape. This is in contrast to the maximum of the Rényi entropy, which misses the optimal PP time due to the existence of excitations with higher energy in the initial thermal state. In the present case, where the structure of the expected EH is known, we can thus also estimate the optimal time t PP by minimizing the distance of the learned EH to the expected shape.

D. Fermions in 2D
We conclude our investigation of PP for free fermions with an example in two spatial dimensions, considering a square lattice at half filling. The corresponding Hamiltonian reads where c ( †) x,y denote fermionic annihilation(creation) operators at sites (x, y) of a two-dimensional (2D) lattice as sketched in Fig. 7(a). Here, we consider a middle square of size n × n as subsystem A (the 'bulk') as indicated by the green rectangle highlighted in Fig. 7(a).
For the PP dynamics, we choose the following deformation of the Hamiltonian given in Eq. (15), with f (x, y) the product of two parabolas along the x and y directions, i.e., see Fig. 7(b). We have checked that an alternative deformation of parabolic type, namely f (r, R = N/2) = (R − r) 2 /R 2 with r the radius from the subsystem's origin, yields similar results in the bulk. Once again, the parabolic deformation is suggested by the CFT version of the Bisognano-Wichmann theorem, which holds also in higher dimensions. Here, we have chosen Eq. (17) because it is better compatible with the square lattice (it is non-negative everywhere and zero at the boundary).
To detect the optimal time t PP to stop the PP, we track the second Rényi entropy S (2) [37], as in the 1D case considered before, but now at one corner of the system shown in Fig. 7(a) as the red dashed square. By analogy with the 1D case, we expect that a small corner (scaling sub-linearly with the full system size N ) will suffice to detect t PP . For our simulations, we take N = 16 and a 4 × 4 corner.
Our results for the simulated PP protocol are presented in Figure 7(c), where we plot the evolution of the bulk infidelity starting from the ground states of the homogeneous finite-size Hamiltonian [Eq. (15)]. Similar to the 1D case, we find that the bulk infidelity for subsystems of different sizes n with respect to the target reduced states of the 2D infinite system decreases significantly, demonstrating the success of PP. Moreover, the minima of the infidelity for different subsystem sizes approximately coincides with the first maximum of the second Rényi entropy, measured in the small corner. At this point, the correlators of the subsystems approach the ones in the infinite system. Figure 7(d) shows two-point correlators with r i = {x i , y i }, that recover the target polynomial decay in the PP state in comparison with the initial ground state.

IV. EXTRACTING CRITICAL EXPONENTS
We now turn to application of PP to QPTs. Quantum critical phenomena can be grouped in universality classes, which are characterized by a set of critical exponents and universal scaling functions [45]. Below we demonstrate that the usage of the PP bulk states for determining critical exponents provides more accurate results than the direct usage of OBC GSs.
For this demonstration, we reconsider the interacting spin chain of Eq.(5) and focus on the transition between the trivial and AF phases along δ = 3 [see Fig. 2(a)]. Since the AF phase spontaneously breaks a Z 2 symmetry (σ z → −σ z ), we expect critical behaviour of the 2D Ising universality class. In this case, the staggered magnetization serves as an order parameter, and for a finite system of size N we consider its root-mean-square (RMS) in a subsystem A of size n ≤ N . At criticality, we assume the standard fine-size scaling (FSS) hypothesis with critical exponents β,ν, scaling functionM , and critical point J c . Note that in our scaling analysis, we use the size of the subsystem (bulk) n rather than the total system size N . For simplicity, in the following, we also use the exact value of ν = 1.  We also investigate how the average values of J c and β, together with their uncertainties, change with increasing the maximum size N . The results are presented in Figs. 9(a) and (b), showing that the usage of the PP mainly improves the accuracy and precision of the critical exponent β, while the results for the critical value J c remains close to those obtained with the OBS GSs. However, we emphasize that in our analysis, we use the exact OBC GSs that are practically impossible to prepare with a quantum simulator. Realistically, one expects not the ground but a low-energy state for which the accuracy of the FSS analysis might degrade due to the presence of the excitations. Applying our PP protocol to such imperfectly prepared GSs will push these excitations out of the bulk to the edges, making the FSS analysis insensitive to imperfections of the initial state preparation.
Finally, we consider the anomalous dimension η, which quantifies the deviation of the scaling dimension of the field from its engineering value. The quantity η determines the polynomial decay of the two-point correlator, which, for the considered critical point, has the form We obtain the values of η by fitting σ z N/2 σ z j with j ∈ [N/2 + 1, 3N/4] in the OBC GSs and in the PP states at the average points J c retrieved from the FSS analysis with M (PP) n,N . The results for N = 100, which are presented in Fig. 9(c), show that the correlator in the OBC GS decays faster than in the infinite state at the same J c , and its functional form deviates significantly from the algebraic one. In contrast, the correlators in the PP state approach the target values practically for all distances. As a result, η obtained from the PP states with increasing N converges much faster to the exact value, as demonstrated in Fig. 9(d).
The above examples illustrate that the PP protocol improves the quantitative studies of quantum phase transitions on finite-size quantum simulators. It provides definite advantages in analyzing the scaling properties of the correlation functions and in extracting the critical indices.

V. APPLICATION: FERMI-HUBBARD MODELS
As an outlook for applications of our protocol to strongly correlated systems, we study the performance of PP for Fermi-Hubbard models (FHMs) for different lattice geometries and parameter regimes. In particular, we consider spinful fermions on a lattice, described by the Hamiltonian Here, the operators c † iσ (c iσ ) create (annihilate) a fermion on site i with spin σ ∈ {↑, ↓}. The first line in Eq. (22) describes the hopping of particles between nearest-neighbour sites ij with transition amplitude J. The second term causes a repulsive (U > 0) onsite interaction if two fermions with opposite spin occupy the same site and the chemical potential µ controls the fermion filling in the ground state. Strongly correlated quantum phases of the FHM have been realized in quantum simulation experiments based on neutral atoms in optical lattices [6,11,47]. Recent developments enable single-site addressing and allow for single-site readout of atoms (see [1]). Moreover, spatial programmability of the Hamiltonian parameters, as required for PP, can be achieved by shaping optical potentials using digital mirror devices [48]. As we will show below, PP enables to reveal the behavior of spatial correlation functions in the thermodynamic limit from experiments on very small lattice systems.
We start by studying the performance of PP for a 10site Hubbard chain [see Fig. 10 (a)] with U/J = 2 and a small amount of hole doping n = 0.8. Here, we focus on a symmetry sector of zero magnetization with N ↑ = N ↓ = 4, where N σ = i n iσ is the number of spin-σ fermions. In order to demonstrate the robustness of our protocol, we take a uniform random superposition of the 10 lowest-lying eigenstates of H FHM as an initial state, which mimics an imprecise preparation of an initial state obtained e.g. in non-adiabatic ground state preparation. The initial state is subsequently time evolved with spatially modified coupling parameters J i = Jg i , U i = U g i−1/2 and µ i = µg i−1/2 , with g i given in Eq. (2). This procedure is repeated 100 times to gather statistics for different initial state superpositions. Figure 10 summarizes the results obtained by averaging over the repetitions. Panel (b) shows the infidelity between a 4-site bulk-subsystem [ Fig. 10(a)] and a subsystem from the center of a 200-site chain. As can be seen, the infidelity decreases by more than one order of magnitude close to the optimal time t PP , which approximately corresponds to 5 hopping times. Fig. 10(c) demonstrated that the optimal time t PP is reflected in a maximum of the entanglement entropy for a 2-site subsystem at the boundary. In Fig. 10 panel (d) we plot the correlation function, analogously defined to Eq. (8), for a N B=8 site bulk subsystem of the 10-site chain at different times and for a bulk in the ground state of a 200site chain (exact). Here, the correlators exhibit a rapid decay with distance modulated by oscillations with approximate 5-site periodicity. As can be seen in Fig. 10(d), these modulations emerge during the quench with the de-formed Hamiltonian and are accurately reproduced close to t PP . Fig. 10(e) shows the absolute relative error of local energy densities in the bulk-subsystem (N B = 8) of the 10-site chain as a function of time. The h i for a single site are given by with c jσ = 0 for j ∈ {0, N + 1}, and the relative error with h i 200 the energy density in the bulk of the 200sites ground state. Note that H FHM of Eq. (22) is given by Fig. 10(e) demonstrates that close to the optimal time, t PP , ∆ rel h i in the bulk vanishes, while a high energy density is accumulated at the edges. As a first proof of principle demonstration for quasi-2D systems, we also studied the performance of PP on a 2 × 5 and 2 × 6-Hubbard ladder at U/J = 2, which is presented in Appendix E. As shown in Figs. 12 and 13 of Appendix E, we observe a decreasing bulk infidelity with respect to large ladders and a corresponding improvement of local observables, providing preliminary evidence in favor of our protocol. We leave a more thorough study of PP on larger ladders, including an investigation of the sensitivity with respect to the Hamiltonian deformation, for future work.

VI. DISCUSSION AND SUMMARY
Quantum simulators represent engineered many-body systems isolated from the environment. Traditionally, analog quantum simulation addresses questions of realizing specific many-body Hamiltonians and preparing and studying corresponding equilibrium and non-equilibrium quantum phases and dynamics. This includes the preparation of ground states or the study of quench dynamics in closed quantum systems [49]. In contrast, the present work outlines quantum protocols running on quantum simulators where the goal is to prepare reduced density matrices as mixed states of many-body systems. The specific example considered in the present paper is the preparation of the reduced density matrix ρ n of a translationinvariant ground state of an infinite system on a quantum simulator, where n is the subsystem of interest and N the size of the simulator (n < N ). In this sense, a finite-size quantum simulator can 'represent' the quantum state of an infinite system, in loose analogy to what is addressed with classical techniques such as iMPS or iPEPS [40,41].
The quantum protocols proposed and studied in the present paper are based on quantum memory of the simulator storing a pure state many-particle wavefunction |ψ , which represents a purification of the desired mixed target state ρ n . The challenge addressed consists in devising coherent dynamics to prepare the target state of interest with given Hamiltonians available on the quantum simulator. As shown in this work, for ground states, this is achieved with coherent quench dynamics generated by spatially deformed system Hamiltonians and starting from an approximate finite-size ground state, as realizable with state-of-the-art quantum simulators [1,2]. The specific choice of a parabolically deformed Hamiltonian is motivated by predictions of the Bisognano-Wichmann theorem on the structure of reduced density matrices for ground states, which then also provides a toolset to monitor the approach to the 'correct' target state in an experiment. The present paper has illustrated both robustness and broad applicability of this scheme of Purification Preparation in numerical simulations for a range of many-body models.
The physical picture underlying the 'cooling' to the target (ground) state ρ n in the 'bulk' of the simulator, with time evolution generated by a parabolically deformed system Hamiltonian, is best understood in a quasi-particle picture. Here, the coherent quench dynamics with the deformed Hamiltonian moves excitations to and accumulates them at the 'edges' of the simulator, which play the role of a quantum reservoir. While, assuming preparation of an initial pure state, the overall dynamics remains in a pure state at all times, it is the entanglement (entropy) of bulk and edges which enables the preparation of the mixed target state ρ n . As demonstrated with examples, this picture also generalizes to thermal quasi-particle excitations, i.e. weakly mixed initial states. A promising application is PP in Fermi-Hubbard models realized with cold atoms [6,8], which are currently limited by the experimentally available techniques to cool fermionic atoms.
Let us briefly comment on the general applicability of our method. Throughout this paper, we have found that a parabolic deformation of the system Hamiltonian, which is motivated by the BW theorem for a relativistic CFT with dynamical critical exponent z = 1, is efficient in preparing the desired purification. However, based on the quasi-particle picture, we expect our protocol to be less efficient for transporting, e.g., low-momentum excitations with dispersion ω(p) ∼ p z with z > 1. Similarly, it remains to evaluate the performance of our protocol in the presence of non-local excitations, e.g., in a topological system. In general, the validity of the quasi-particle picture in the presence of interactions could be tested by deriving a Boltzmann-type transport equation for general Hamiltonian deformations. We leave a thorough investigation of these issues for future work. In this context, it could be useful to include an additional step optimizing the shape of the deformation into our protocol.
We conclude by commenting on the relation of Purification Preparation (PP) to algorithmic cooling, and we refer to the review [50], and experimental demonstrations in nuclear-magnetic-resonance based and ultracold atom platforms see [51] and [52], respectively. The main idea of so-called heat-bath algorithmic cooling discussed in a quantum information context consists of extracting entropy from a subset of logical qubits and compressing it into "reset" qubits, which in turn release the entropy by coupling to a heat bath. While the compression step is superficially similar to the quasi-particle transport within our PP protocol, our goal is very different. In contrast to algorithmic cooling, which aims to prepare a fresh set of pure qubits with minimal entropy, our target is a specific mixed state corresponding to a part of the infinite-size ground state. It would be interesting to study the inclusion of dissipative elements into PP to iteratively extract the compressed quasi-particles from the boundaries of the system. The dynamics of the deformed free-fermion Hamiltonian in 1D of Sec. III, permits an illuminating continuum interpretation. To reveal this, we decompose the lattice operators as where we assume that ψ x becomes a smooth functions in the limit a → 0 for suitably chosen momentumj. From the Heisenberg equation of motion, we have i∂ t ψ n ∼ J e ij g n,n+1 ψ n+1 + e −ij g n−1,n ψ n−1 . (A3) We now expand the field operators in derivatives up to second order, and assume that an analogous expansion holds for the deformation which approximately (to leading order in a) captures the operator evolution for quasi-particles aroundj under a smooth deformation g x . Note that for g x = 1, we recover the usual left-and right-moving modes for a given fermion filling determined byj, as it should be. The above continuum approximation allows us to gain further insight for a given deformation g x . In the following, we focus on the case of half filling and consider the right-moving modes atj = 3π/2 (the left-movers only differ by some signs). Employing the method of characteristics to solve the partial differential equation for ψ x , it follows that the quasi-particles travel on trajectories from x 0 to x 1 determined by with a position-dependent velocity v x . We emphasize that the term ∝ ψ x ∂ x g x ensures a proper normalization of ψ x , but it does not affect the trajectories. Indeed, for j = 3π/2 (A) is equivalent the Dirac equation for massless fermion propagating in a optical metric [53]. Next, we compare the efficiency of the PP for a system of size L evolved with a parabolic deformation or a SSD, In these cases, the involved integrals can be evaluated analytically. For the parabolic deformation, we obtain the trajectories (x 0 = −L/2) as where the approximation is valid asymptotically for t L/(4v) log [(L − 2x 0 )/(L + 2x 0 )]. Similarly, for the SSD, we obtain where the asymptotic approximation applies for t L/(πv) tan πx0 L and t L/(π 2 v). This proves our claim made in the main text that the parabolic deformation is more efficient than the SSD in transporting the quasi-particles towards the edge, as the former shows an exponential approach ∝ exp (−4vt/L), while the latter is only algebraical ∝ L/(vt).
These results also explain the mode mapping discussed in Sec. III B, i.e., the transformation of initial lattice operators in momentum space to localized operators in position space at the edges of the system. Their reflection at times t L/v is a lattice effect that happens when the ψ x operators approach a delta-like shape ∼ δ(x ± L/2) and the continuum approximation breaks down. Neglecting this reflection, we can estimate the weight w ∆x (t) of the operator c p=j (t) contained in an edge region of size ∆x by summing up all trajectories that have reached this region, i.e., where Θ(. . . ) is the Heaviside step function and we integrate of the initial conditions of the trajectory x 1 (t).
For the above results for the parabolic and SSD, we see that the characteristic curves never cross, such that w ∆x (t) = 1/2 − x 0 (t)/L, where x 0 (t) is obtained by inverting x 1 (t) = L/2 − ∆x. Explicitly, w SSD ∆x (t) = Given w ∆x , the probability density for finding a quasiparticle at position x ∈ (−L/2, L/2) and time t is given by |χ x (t)| 2 = ∂ ∆x w ∆x | ∆x=L/2+x . For comparison with the numerical results, we estimate the contribution of c p=j (t) to a correlation function · · · c p=j (t) · · · in a sub- Repeating the calculation for p = π/2 gives an analogous result for X p=π/2 (t). For the initial OBC GS discussed in the main text, the dominant contribution comes from correlations between left-and right-movers. The error in corresponding two-point functions evaluated in the bulk A, originating from the initial correlator c † p=π/2 c q=3π/2 , is thus estimated to decrease as the product X p=π/2 (t)X p=3π/2 (t). To compare these continuum predictions with the simulation results for the lattice model, we replace L = N a, v = 2Ja and x a/b = (n a/b − N/2)a, where N is the number of lattice sites for a subsystem consisting of the sites {n a + 1, n a + 2, . . . , n b − 1, n b } ⊂ {1, 2, . . . , N − 1, N }. Specifically in Fig. 4(a) we considered N = 36 and n a = N − n b = 4. The analogous curve in Fig. 4(a) for the SSD involves integrals which we evaluated numerically. sketched in Fig. 11(a) where we consider Trotterization of the parabolic deformation [Eq. (2)]. In the plot, each box corresponds to an elementary Trotter steps with widthñ and height 2(δt)ñ. Because of the experimental limitations, we assume the shortest time step accessible in the experiment from which one can obtain the smallest sizẽ n min . To achieve a Trotter error that scales quadratically with the Trotter steps, we halve the rectangles along the time axis and reorder them as shown in the figure. We also move the first elementary step (red rectangle), having the shortest evolution time, to the end of the sequence and stacked it with the identical last step. This only affects the very first elementary step. The whole procedure gives the Trotter-Suzuki expansion of the PP at second order in the composite step ∆t(ñ min ).
In Fig. 11(b), we give the result of the Trotterized PP for a N = 36 system of free fermions starting from the OBC GS. We considerñ min = 2 andñ min = 8 cases, where we fix the time step to (δt) 2 = 0.01, where we measure time in units of the inverse coupling strength in the original Hamiltonian. As a result of the Trotterization scheme, forñ min = 2, the smallest Trotter step becomes 2(δt) 2 = 0.02 (stacked elementary step between two composite steps), and, forñ min = 8, the smallest elementary step takes time (δt) 10 = 0.09. In both cases, the minima of the infidelity determines optimal times close to the one obtained with the analog evolution, Fig. 3(c). In the considered scenarios, the optimal evolution takes 5 composite Trotter steps equivalent to 170 and 140 elementary trotter steps in the case ofñ min = 2 andñ min = 8, respectively. Notice that the latter can be obtained from the former by eliminating the smallest middle steps. By doing so, we effectively realize Trotterization of a different deformed Hamiltonian with a plateau in the bulk, which explains the degradation in the efficiency of the PP observed forñ min = 8. with and N c the number of ladder rungs. The index λ ∈ {1, 2} denotes the two legs of the ladder. Note that in the definition of the deformed Hamiltonian Eq. (E1) the index i in the first line runs from 0 to N c − 1, while it runs up to N c in the remaining lines. The coefficients g i−1/2 are thus required to generate a symmetric parabola in the x-direction. We repeat the quench simulations 50 times in order to collect statistics for different random initial states. We emphasize that by sampling over random initial superpositions we demonstrate PP for two types of typical experimental imperfections. Since a single realization corresponds to a coherent superposition with of low-lying excited states, it models an imperfect preparation of the ground state, e.g. via almost adiabatic state preparation. The collection of all such realizations realizes a mixed state, akin to a microcanonical ensemble with a small energy window above the ground state. This ensemble thus models another type of experimental errors, similar to a finite (low) temperature. In either case, our simulations indicate the PP works well for these states "close" to the ground state. In Fig. 12(b) we plot the Uhlmann-infidelity between the subsystem B depicted in Fig. 12(a) and a corresponding bulk-subsystem in the ground state of a 2 × 100ladder. We observe a continuous decrease of the infidelity, reaching a minimum at about ∼10 hopping times. Fig. 12(c) demonstrates that the time of optimal infidelity is approximately reflected by a maximum in the entanglement entropy for a boundary subsystem A shown in Fig. 12(a).
We now turn to analyzing local observables at different times during the evolution with the deformed Hamiltonian. Fig. 12(d) shows nearest-neighbor hoppingcorrelations averaged over initial states (t = 0) and over states at the optimal time (t = t PP ). As can be seen, correlations close to the system boundary are closer to the exact value at t = t PP . While the bulk-correlators are already close to the exact value at t = 0, the correlators at t = t PP have a significantly reduced error bar, implying that the exact correlations are consistently reproduced close to the optimal time.
Finally, we investigate the behaviour of nearestneighbor d-wave pair correlations during the time evolution with the deformed Hamiltonian (E1). For strongly correlated systems of Fermions, the decay of the d-wave pair correlation function can serve as a detector for pairing phases [59]. Thus, it plays an important role for characterising low-temperature phases in Fermi-Hubbard models. For a 2-leg ladder, the d-wave pair correlation function is defined as ,↑ creates a singlet on rung i [59]. In Fig. 12(e) we plot the real part of the averaged nearest-neighbor d-wave pair correlations D i = Re ∆ † i ∆ i+1 for all 4 links of the 2 × 5 ladder. For comparison, the data are divided by the averaged value of D i from a 10-site bulk region of a 2 × 100 ladder. Similar to the hopping correlators discussed above, we obtain large fluctuations for the correlators close to the boundary, while the mean value of the correlators in the center approaches the exact value with a significantly reduced error bar.
We amphasize that the results presented in this chapter where carried out for a small ladder system (2 × 5) and thus may be severely affected by finite-size effects. Future investigations will have to include a more sophisticated analysis for larger ladders in several different parameter regimes as well as an investigation on how strongly the results depend on the type of Hamiltonian deformation.

PP on a 2×6 Fermi-Hubbard ladder
In this section we repeat the simulations presented in Fig. 12 for a 2 × 6 strongly repulsive Fermi-Hubbard ladder. While the filling is chosen to be n = 10/12 [see Fig. 13(a)], the quench parameters are exactly the same as described in the previous paragraph. Again, we observe a significant drop of infidelity and corresponding nearest-neighbor correlation functions approaching the values in the N c = 100 ladder.

Fermi-Hubbard Chain at Half-Filling
In this section we provide additional results for 1D Hubbard chain in a different parameter regime. The setting is equivalent to the one discussed in section V in the main text with the difference that we study the system at half filling N ↑ = N ↓ = 5 and U/J = 1. Overall, the performance of PP is very similar to the situation discussed in Fig. 10 in the main text. Fig. 14(e) shows that the power-law for the decay of correlations is accurately reproduced close to the optimal time t = t PP which we define to be at the maximum of the entanglement entropy Fig. 14(c). Analogous plots to Fig. 12 for a 2×6 Fermi-Hubbard ladder. The initial state is a random superposition of the 2 lowest-lying energy eigenstates. The quench is repeated 10 times to gather statistics for different random initial states.