Disentangling Sources of Quantum Entanglement in Quench Dynamics

Quantum entanglement may have various origins ranging from solely interaction-driven quantum correlations to single-particle effects. Here, we explore the dependence of entanglement on (time-dependent) single-particle basis transformations in fermionic quantum many-body systems, thus aiming at isolating single-particle sources of entanglement growth in quench dynamics. Using exact diagonalization methods, for paradigmatic non-integrable models we compare to the standard real space cut various physically motivated bipartitions. Moreover, we search for a minimal entanglement basis using local optimization algorithms, which at short to intermediate post-quench times yields a significant reduction of entanglement beyond a dynamical Hartree-Fock solution. In the long-time limit, we identify an asymptotic universality of entanglement for weakly interacting systems, as well as a cross-over from dominant real-space to momentum-space entanglement in Hubbard-models undergoing an interaction quench. Finally, we discuss the relevance of our findings for the development of tensor network based algorithms for quantum dynamics.

Entanglement as one of the most fundamental traits of quantum mechanics plays a key role in various physical contexts, ranging from the characterization of complex many-body states [1, 2] to its use as a resource for quantum technologies [3][4][5][6][7]. On the other hand, strong entanglement in correlated quantum matter poses a significant challenge that limits both its analytical and numerical description. This becomes particularly relevant in nonequilibrium dynamics, where entanglement generically grows rapidly with time [8][9][10][11][12][13][14][15][16][17][18][19][20]. This growth, however, may originate from mechanisms of different inherent complexity: While it can already be found in analytically solvable free theories, at the opposite end entanglement can be purely interaction-induced. Yet, it has remained largely unexplored as to what extent a heterogeneous dynamical entanglement content may be distinguished according to its various sources on general grounds, which would substantially contribute to the understanding of complex quantum matter. Furthermore, such a distinction might be crucial for the development of numerical algorithms for quantum dynamics based on tensor network approaches [21][22][23][24][25][26][27][28][29][30][31][32][33][34], whose efficiency crucially depends on entanglement scaling.
The purpose of this work is to investigate how the single-particle content of entanglement in quantum quench dynamics can be isolated from genuine many-body complexity, thereby revealing physically distinct sources of quantum correlations. To this end, we study the dependence on generally time-dependent single-particle basis rotations (entanglement-cuts) of the half-system entanglement entropy where ρ At = tr Bt ρ denotes the reduced density matrix of subsystem A t at time t obtained from the full density matrix ρ = |ψ(t) ψ(t)| by tracing out its complement B t (see Fig. 1(a)), and |ψ(t) is the state of the system at time t. As we exemplify for several non-integrable fermionic quantum many-body systems, for quite long arXiv:1905.05782v2 [cond-mat.str-el] 16 May 2019 transient times the entanglement entropy indeed exhibits a significant basis-dependence. In particular, we compare to the real space representation various physically motivated entanglement cuts, corresponding to time-evolving Wannier functions (of a dynamical Hartree-Fock solution), and Bloch functions (see Fig. 1(b)). Based on our concrete case studies in the framework of exact diagonalization, in this context we identify several general principles for both transient and asymptotically long times, including (i) an asymptotic universality (in the sense of basis-independence) of entanglement in systems with weak to modest interactions, and (ii) a cross-over between real-space and momentum-space entanglement in Hubbard models undergoing an interaction quench. Furthermore, we search for an optimal (minimal entanglement) orbital basis using local optimization algorithms (see inset in Fig. 1(b)), in order to separate genuine many-body entanglement content, originating from quantum correlations involving several particles, from single-particle contributions. Entanglement dynamics from Krylov time propagation. For our quantitative study, we consider a system of n p interacting fermions on a lattice with sites j = 1, . . . , L, which are annihilated by the vector of field operators c = (c 1 , . . . , c L ). To quench the system out of equilibrium, the Hamiltonian is suddenly changed from H i to H f at time t = 0, assuming that the many-body state |ψ(t = 0) was prepared as the ground state of H i at t < 0. Using numerically exact Krylov time-propagation methods [35,36], for several examples of H i and H f representing quenches into non-integrable systems, we then compute the time-evolved wave function ( = 1) In order to study the dynamics of the (von Neumann) entanglement entropy S(t), we decompose the system into two subsystems, A t and B t , which each contain L/2 orbitals and on average n p /2 particles. However, we do not restrict this decomposition (entanglement cut) to the real space lattice space, but consider time-dependent basis orbitalsc t = u t c, resulting from the real space lattice representation by an arbitrary U(L)-transformation u t (with U(L) denoting the group of L×L unitary matrices). To represent the time evolved wave function in an arbitrary orbital basis, we need to express u t in many-body Hilbert space, where it is denoted by U t . Naively applying U t to the real space representation ψ(t) of the state vector |ψ(t) would generate an unfeasible computational cost ∼ N 2 , where N is the dimension of the many-body Hilbert space of n p particles. However, exploiting that U t formally may be seen as a "time-evolution" with respect to the the free HamiltonianH t = j, (i log u t ) j, c † j c yieldsψ which naturally allows for the efficient (O(N )) computation of the transformed representationψ(t), again using Krylov time-progation with respect to the fictitious HamiltonianH.
In the following, we present exact numerical data on the quench dynamics of systems with up to L = 24 sites and n p = 12 fermions, where N ≈ 2.7 · 10 6 . Concretely, we study two examples of non-integrable ergodic model systems, (i) a 1D two-band band insulator that is quenched to a topological insulator (TI) phase with weak to modest interactions, and (ii) a metallic singleband Fermi-Hubbard model with next-nearest neighbor (NNN) interactions. By comparison of these physically quite different scenarios, we will identify various generic features relating to the basis-dependence of entanglement in quantum quench dynamics.
(i) Weakly interacting topological insulator. We first discuss a one-dimensional (1D) two-banded system similar to the celebrated Su, Schrieffer, and Heeger (SSH) model [37,38]. There, the lattice index j = 1, . . . L is decomposed into a unit-cell index i = 1, . . . , L/2 and a sublattice index α = a, b such that j = 2i − 1 for site a in cell i and j = 2i for site b in cell i, respectively. In reciprocal space, the two-band Bloch Hamiltonian at lattice momentum k = 4πp/L, p = 0, . . . , L/2 − 1 of the system is defined as where σ i are the standard Pauli matrices acting on a/b sublattice space, m(t) plays the role of a Dirac mass parameter that is quenched from m i to m f at t = 0, and J is the hopping strength between neighboring unit cells. For |m| < |J|, the model is in a topological insulator phase protected by the chiral symmetry σ z h(k)σ z = −h(k), while at |m| = |J| a topological quantum phase transition to a trivial band insulator phase extending over the parameter regime |m| > |J| occurs. Here, we focus on quenches which change the topological phase of the system Hamiltonian, by choosing m i in the trivial, and m f in the non-trivial phase, respectively. For the (standard) real space entanglement cut, even the non-interacting system is found to exhibit linear entanglement growth at small times. However, this entanglement is readily seen to be entirely basis-dependent: The exact solution |ψ(t) stays a single Slater determinant at all times, leading to zero entanglement in a suitable basis obtained from time-evolved Wannier functions or Bloch functions representing the initial uncorrelated state. This provides a conceptually simple example for how strongly the entanglement of a system far from equilibrium may depend on the chosen representation.
We now address the natural question as to how an integrability breaking interaction term affects this phenomenology. To this end, we add to the free model (4)

the interaction Hamiltonian
where already a non-zero U or W individually is sufficient to make the model nonintegrable, and both terms preserve the protecting chiral symmetry. We first discuss the case on U = J, W = 0 for transient times (see Fig. 1(b)). For the Bloch and Wannier entanglement cut, in which the non-interacting solution would exhibit zero entanglement, S is lower than in the real space cut, for short to intermediate times. Interestingly, evolving the Wannier basis with respect to a dynamical Hartree-Fock Hamiltonian, thus accounting for interactions at mean field level leads to a substantial reduction of S for quite long times. Searching for an optimal entanglement cut by treating the translationinvariant Hermitian matrixH t as a variational ansatz generating the basis transformation U t (see Eq. (3)) [36], we are able to further reduce S significantly below the dynamical Hartree-Fock basis (see inset of Fig. 1(b)) by using local optimization algorithms such as Gradient Descent and ADAM [39]. We now switch on W > 0 which is generally found to speed up the process of thermalization in 1D topological insulator models [40], thus allowing us to look at the long-time (close to equilibration) behavior of our model. For weak to modest interactions, we observe a quite remarkable universal aspect of entanglement in the long-time limit, namely that S becomes largely independent of the choice of basis (see Fig. 2). This behavior exemplifies a quite generic mechanism that can be understood with the following physical picture [36]. When a system is quenched far from equilibrium by changing its gapped single particle Hamiltonian, weak interactions are expected to lead to thermalization at an effective temperature β e of the system with respect to the close to non-interacting post-quench Hamiltonian. However the thermal entropy of a free Hamiltonian (such as Eq. (4)) generically is independent (up to non-extensive boundary effects) on the basis of a bipartition with an average of n p /2 particles in each subsystem. Thus, the basis-dependence of the entanglement entropy in this scenario is a transient effect, even though it might last to long times.
(ii) Spinless NNN Fermi-Hubbard model. As a second non-integrable model system, we consider a spinless single band Fermi-Hubbard model with NNN interaction at filling ν = n p /L, defined by the Hamiltonian where J (fixed to 1 in the simulations) denotes the hopping strength, n j = c † j c j , and V (t) ≥ 0 represents the repulsive NNN interaction strength which is suddenly changed from 0 to a finite value over the quench at t = 0. We focus on ν = 1/3 in the following, where the model (6) at zero temperature is known to stay in a metallic Luttinger liquid phase up to large V > 0 [41], and has been found to exhibit strongly ergodic behavior [42]. When comparing the long-time entanglement entropy of this model in real space and momentum space, we find an interesting basis-dependence that persists even in the long-time limit: For weak to modest V , the longtime entanglement in momentum space is lower than in real space, whereas at strong correlations (V 3J), the real space cut leads to lower entanglement entropy (see Fig. 3). This behavior may be interpreted as a dynamical manifestation of the generic competition between an interaction term that is diagonal in real space and a free band structure that is diagonal in momentum space in Hubbard models. At short times, by contrast, entanglement for all V grows steeper in the momentum cut, owing to the non-local nature of the interaction in reciprocal space. In both cuts, entanglement clearly exhibits extensive (volume law) scaling for long times, as expected for an equilibrating ergodic system (see inset of Fig. 3).
Concluding Discussion. We conclude by putting our present results into a broader perspective from several angles. Generally speaking, non-integrable ergodic systems are known to exhibit a generic dynamical behavior for the entanglement entropy S in coherent quench dynamics starting from a weakly entangled (equilibrium) state: S initially grows linearly in time reflecting a ballistic spreading of quantum information at a certain velocity [8][9][10][11][12][13][14][15][16][17][18][19][20]. At large times, close to thermalization, S approaches an extensive value (volume law), i.e. S = αV A , where V A denotes the volume of subsystem A, since the entropy as an equilibrium thermodynamic potential is always extensive. While all of our present findings fit into this rough generic framework, here we identified and microscopically exemplified several general principles relating to the basis dependence and heterogeneous physical nature of dynamical entanglement growth.
First, for systems with weak to modestly strong interactions that are quenched out of equilibrium by a change in their band structure, the transient entanglement exhibits very strong basis dependence and may be reduced significantly beyond a dynamical mean field solution by means of single-particle basis rotations. By contrast, the volume-law coefficient α of the long-time average of S becomes basis independent. This is because the thermal entropy of the asymptotically equivalent nearly free ensemble typically is largely independent of the choice of bipartition, as long as each subsystem contains half of the degrees of freedom and accommodates on average half of the particles.
Second, in Hubbard models that are undergoing an interaction quench leaving the non-interacting Hamiltonian unchanged, a significant basis-dependence of entanglement persists in the long-time limit, where the considered ergodic systems exhibit basis-dependent extensive entanglement scaling: Up to intermediate post-quench interaction strengths entanglement is found to be weaker in reciprocal space while at stronger correlations the realspace cut has lower in S. This behavior dynamically reflects the competition between kinetic energy and interactions characteristic for Hubbard models.
These findings on distinguishing single-particle contributions to S from inherently more complex quantum correlations are not only of fundamental interest, but might also serve as important guidelines for the pursuit of devising new numerical algorithms for quantum dynam-ics. Within the realm of thermal equilibrium, the observation that low-temperature states typically exhibit significantly lower (area-law [43,44]) entanglement than generic states has been crucial for systematically taming the exponential complexity of correlated systems, e.g. with the advent of tensor network methods [21][22][23][24][25][26][27][28][29][30][31][32][33][34], where the basis-dependence of entanglement has been successfully exploited in the context of spin-boson models [45] and quantum chemistry settings [46][47][48]. However, far from equilibrium, generically the dynamical proliferation of entanglement eventually builds up an exponential wall even for the most powerful known computational methods. In this context, our results provide systematic insights addressing the question as to what extent both transient entanglement growth and long-time entanglement may be made computationally manageable, e.g. by extending tensor network methods to a flexible representation of the physical degrees of freedom that allows for time-dependent basis choices.
In this work, we have been concerned with the influence of single-particle basis rotations on entanglement dynamics. At a conceptual level, our search for an optimal basis in which S is lowest may be interpreted as entanglementbased dynamical mean-field approach aimed at minimizing the single-particle contributions to S. Thinking further along these lines, it is readily conceivable to consider more elaborate approximate solutions based on genuine many-body methods as a starting point for the choice of a correlated basis, which may eliminate even some higherorder quantum correlations from the system. Quantitatively exploring this approach and constructing hybrid numerical methods, where a physically motivated approximate solution serves as the starting point to be augmented by an unbiased tensor-network simulation, represents an interesting subject of future work.

Krylov Time Propagation
Here we briefly elaborate on the Krylov time propagation method [A1], used in the main text for the timeevolution as well as for the single-particle basis rotation of the many-body state |ψ . We start with the case of conventional time-evolution, where the goal is to calculate |ψ(∆t) = e −iH∆t |ψ from a known initial state |ψ , with a Hamiltonian H that is assumed to be time-independent. The idea of Krylov time propagation is based on the observation that the time-evolved state |ψ(∆t) , up to a given order N K − 1 in ∆t reads as and thus manifestly belongs to the N Kth Krylov subspace K NK (H; |ψ ) = span |ψ , H|ψ , H 2 |ψ , ... , H NK−1 |ψ .
Using this insight, the Hamiltonian H may be expressed as a tri-diagonal matrix T in the N K -dimensional Krylov subspace, where typically N K ≈ 50. Then, the matrix exponential e −iT ∆t may readily be expressed and explicitly applied to |ψ within this Krylov subspace, before transforming the resulting vector back to the full Hilbert space.
Practically, one needs to generate a basis for K NK (H; |ψ ), which can be done using the Lanczos method (see e.g. Ref. [A2]). We denote with Q the matrix whose columns are the N K orthonormal basis states |q n (with n = 0, ..., N K − 1) for the Krylov subspace, where |q 0 = |ψ . The initial state |ψ is thus represented in the Krylov subspace by the N K -component vector Q † |ψ = (1 0 ... 0) . The vector e −iT ∆t Q † |ψ therefore represents the time-evolved state |ψ(∆t) in K NK (H; |ψ ) truncated to O(∆t NK−1 ). Hence, the full time-evolved state is retrieved from at an overall computational cost of O(N K N ), where N denotes the dimension of the many-body Hilbert space.
To apply this method also to basis transformations, we note that it can efficiently perform any unitary transformation U = e −iH on the state |ψ , provided that the Hermitian many-body operatorH can be efficiently applied to a state vector. Turning to the single particle transformationsc = u t c with u t ∈ U(L) used in the main text, the Hermitian many-body operatorH t can be constructed as the following fictitious non-interacting Hamil-tonianH t = j, (i log u t ) j, c † j c , where the c j denote the annihilation operators for the j-th orbital. This reduces the intuitive cost of O(N 2 ) for a unitary basis transformation to the lower computational cost of O(N ) that is typical of Krylov methods, thus removing the computational bottleneck of our simulations from the application of basis transformations.

Entanglement Optimization
Here we provide some details on the variational optimization of the entanglement entropy S used in the main text. The goal is to find an optimal single-particle basis transformation u ∈ U(L) yielding the minimum entanglement entropy for a equal-size bipartition of the resulting single-particle modes (c 1 , . . .c L ), where we dropped the time-dependence of all involved quantities for notational brevity. For a given set of L single-particle modes corresponding to u, we denote withÃ andB the two subsystems containing L/2 modes each. The (von Neumann) entanglement entropy forÃ is then written as where ρÃ = trB|ψ ψ|.
As mentioned in the main text the representation (coordinate) vector ψ = (ψ 1 , . . . , ψ N ) of the state |ψ transforms under the basis rotation asψ where againH = j, (i log u) j, c † j c . For our variational minimization of S(ρÃ), we interpret the elements (i log u) j, ≡h j, as fictitious hopping amplitudes ofH, and use them as variational parameters. We furthermore would like to impose the constraint that each bipartition should contain on average n p /2 particles, i.e. half of the total number of particles n p . To this end we minimize the entanglement cost function where nÃ = tr(NÃρÃ) is the expectation value of the particle number operator NÃ in subsystemÃ. The second term with λ > 0 in the above equation punishes the entanglement cuts resulting in an imbalanced particle number between subsystemsÃ andB, where we choose λ ≈ 40 to basically obtain a hard constraint for practical purposes.
For the minimization we adopt gradient based methods such as simple steepest descent or the Adam optimization algorithm [A3]. Such gradient methods require the calculation of the derivatives of C(ρÃ) with respect to all the parametersh j, , which we do numerically to first order. Imposing translation invariance onh j, automatically leads to Wannier-like (i.e. shift-orthogonal) single-particle orbitals, and reduces the number of real variational parameters from O(L 2 ) to O(L). Denoting by T the translation operator on the lattice, this ansatz limits our search to single-particle basis sets satisfying T u T † = u. We point out here that this choice excludes the physically relevant basis choice of Bloch states that are eigenstates of the translation operator T , satisfying T u = u D with D a diagonal matrix of phase factors. However, our simulations using the unrestricted set of parametersh j, for small systems indicate that the momentum cut of Bloch states always constitutes an isolated local minimum for the entanglement cost function (A5). Hence, for gradient-based optimization methods, the minimization within the space of Wannier-like orbitals, where the ansatz of translation invariance applies, is more promising and indeed found to be more successful.
Finally, we point out that the translation operator T in real space is different for the two models studied in the main text. In the case of the SSH model, T translates by one unit-cell, which consists of two orbitals, hence the ansatz Th T † =h yields 2L independent variational parameters. In the case of the spinless Hubbard model with NNN interaction, T translates by one physical site with a single orbital, hence resulting in an ansatz depending on L independent variational parameters.

Asymptotic Basis-Independence of Entanglement Entropy
In this section, we briefly elaborate on the long-time basis-independence of the entanglement entropy observed Figure A4: Schematic representation of the equivalence of the spectra of real-space and momentum restrictions of H f . The single-particle spectrum of H f for the total system is shown by the thick black lines, on the right labeled by the momentum k. For the real-space cut A (region in the blue box) the spectrum {εν } of H f |A is shown by the dashed blue line. For the momentum cutÃ, the spectrum {εν } of H f |Ã is shown by the orange thick lines. The two spectra concur, up to boundary modes (the two light-blue dots in the spectrum on the left), and corrections that vanish in the thermodynamic limit.
in the post-quench dynamics of the thermalizing SSH model. As explained in the main text, when the band structure of the system is suddenly changed by quenching its single-particle Hamiltonian from H i to H f , the presence of weak interactions is generically expected to induce thermalization of the system, as intuitively the scattering processes will equilibrate the single-particle excitation created during the quench protocol [A4].
In the weakly interacting regime we can hence assume that in the long-time limit the system to good approximation thermalizes with respect to the free part of the post-quench Hamiltonian H f , in the sense that a subsystem A would be described by a reduced density operator of the form where β e denotes the inverse temperature, generally depending on the density of excitations created during the quench, and H f | A denotes the restriction to subsystem A of the free post-quench Hamiltonian H f . Since H f | A is a free fermionic Hamiltonian, the (von Neumann) entanglement entropy S(ρ A ) = −tr[ρ A log ρ A ] can be calculated as the thermal entropy of a free Fermi gas, which reads as (see e.g. Refs. [A5, A6]) S(ρ A ) = ν log 1 + e −βeεν + β e ε ν 1 + e βeεν (A7) where the ε ν denote the single-particle energies of H f | A .
Since the above expression of S(ρ A ) depends only on the single-particle eigenvalues ε ν , we argue in the fol-lowing that for a broad class of Hamiltonians H f the value of S(ρ A ) is independent of the choice of the subsystem A up to non-extensive boundary terms, as long as A has the same size and contains, on average, half of the particles. Specifically, let H f be the Hamiltonian of a system with an inversion symmetry I such that the Bloch Hamiltonian h f (k) in momentum space satisfies I h f (k) I † = h f (−k), hence its spectrum is symmetric with respect to k = 0. As an example of quite different entanglement cuts, consider the subsystem A as of consisting of half the system in real space, andÃ a momentum space partition comprising only Bloch states in the momentum interval IÃ = [−π, 0). Notice that this is the case for the SSH Hamiltonian and the entanglement cuts studied in the main text. If the system is thermalized with respect to H f , then the reduced density operator forÃ reads as ρÃ = Z −1 A e −βeH f |Ã . The crucial observation is that, under the above assumptions, the spectra of H f | A and H f |Ã are the same up to nonextensive boundary terms (e.g. topological edge states at the boundaries of A) and other corrections which vanish in the thermodynamic limit, provided that A andÃ contain the same number of single-particle modes and, on average, the same number of particles. This is schematically shown in Fig. A4. Hence the entanglement entropies S(ρ A ) and S(ρÃ) calculated from Eq. (A7) are equivalent, up to terms that do not scale with the size of the subsystems. Similar statements hold when comparing to a cut of generic Wannier-like functions. This explicates the asymptotic basis independence found in the main text for the quenched weakly interacting SSH model, and shows how it can be generalized to other weakly interacting quenched tight-binding models.