{'fraction_non_alphanumeric': 0.03289384259888541, 'fraction_numerical': 0.011213809976892756, 'mean_word_length': 4.391901795529498, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics of the entanglement spectrum in spin chains

We study the dynamics of the entanglement spectrum, that is the time evolution of the eigenvalues of the reduced density matrices after a bipartition of a one-dimensional spin chain. Starting from the ground state of an initial Hamiltonian, the state of the system is evolved in time with a new Hamiltonian. We consider both instantaneous and quasi adiabatic quenches of the system Hamiltonian across a quantum phase transition. We analyse the Ising model that can be exactly solved and the XXZ for which we employ the time-dependent density matrix renormalisation group algorithm. Our results show once more a connection between the Schmidt gap, i.e. the difference of the two largest eigenvalues of the reduced density matrix and order parameters, in this case the spontaneous magnetisation.


Introduction
In the last decade, entanglement is revealing itself an extremely useful tool in the characterisation of quantum many-body systems [1]. Since the full-tomography of a quantum state would require an exponentially large number of measurements, measures of quantum correlations, including entanglement, give a lot more information than expectation values of local operators. Many works have concentrated on bipartite entanglement including von Neumann and Renyi entropies [2,3], mutual information [4], negativity between blocks of particles [5], concurrence between two spin-1/2 particles as first studied in [6,7]. Beyond entanglement, the study of discord and its global multipartite extension in spin chains have also received recently a lot of attention [8,9,10]. The multipartite setting is far less explored because of its complexity compared to the bipartite scenario, although a few works show the emergence of genuine multipartite entanglement close to quantum phase transitions in spin chains [11,12].
In the last few years, entanglement spectrum (ES), consisting of the eigenvalues of the reduced density matrix of one of the two parties in a bipartition of the Hilbert space, for example in real or momentum space has received enormous interest in the condensed matter physics community. For strongly correlated systems, ES was studied in relation to the performance of the density matrix renormalisation group [13]. Since Li and Haldane's seminal paper [14], many authors have investigated the ES in a variety of settings including one-dimensional spin chains [15]. In particular, in critical 1D systems described by a conformal field theory, the density distribution of the eigenvalues of the reduced density matrix is universal and only depends on the central charge of the corresponding conformal theory [16]. However, it is important to stress that the discrete part of the ES, containing the largest eigenvalues, does not follow this universal distribution, and thus cannot be predicted using the von Neumann entropy. Instead this region of the ES can be connected to the energy spectrum of a boundary conformal field theory [17]. Other studies have shown that the analysis of the individual eigenvalues gives much more information than it was anticipated. In particular the so-called Schmidt gap, the difference of the two largest eigenvalues, detects very accurately the position of a critical point and scales with universal critical exponents [18,19,20,21,22,23]. The Schmidt gap can also be employed in numerical simulations to reveal Haldane phases that lack a local order parameter. As in these hidden order phases the ES must be formed of evenly degenerate multiplets [24], the Schmidt gap is bound to close [25]. Since this information is not accessible in the entanglement entropy, the study of the ES has deepen the understanding of many-body quantum states.
The dynamics of strongly correlated systems has also received particular attention recently. On one hand, one question is whether or not isolated quantum systems reach a steady state after a sudden change of the Hamiltonian (for a review check [26]) and whether this state resembles a thermal state. Ref. [27] studied the evolution of the reduced density matrix in the quantum Ising model. On the other hand, quasi adiabatic sweeps of a Hamiltonian parameter across the transition reveals universal dynamics ruled by the Kibble-Zurek mechanism (see [28] for recent reviews).
In this work we analyse the time evolution of the ES in a system of 1D spin chains after instantaneous or quasi-adiabatic quenches in one parameter of the Hamiltonian. We first study the 1D Ising model that can be exactly diagonalised and we show results on the evolution of the ES after a sudden change of the magnetic field. We show that when the magnetic field is changed from the paramagnetic phase to the ferromagnetic phase, the ES levels cross each other and the corresponding Schmidt gap oscillates. The oscillations are equivalent to the oscillations of the spontaneous magnetisation after a similar quench. We also analyse the evolution of ES after a quasi-adiabatic quench and show the appearance of an algebraic dependence of the Schmidt gap with the rate of change of the magnetic field when crossing the quantum phase transition.
We then move to the XXZ model in which we quench the value of the anisotropy parameter. The unitary evolution of the ground state in this case is calculated using the time-dependent density matrix renormalisation group (tDMRG) algorithm [29]. In this case however we do not observe crossing of the eigenvalues but rather a tendency of them to become degenerate and at the same time decreasing exponentially with time.
The paper is organised as follows: in Sec. 2 we present the models we consider in this work and review the basic properties of the ES; in Sec. 3 we present our results for the dynamics of the Ising model; in Sec. 4 we pass to the ES of the XXZ model and in Sec. 5 we summarise and conclude. In the Appendix, details of the diagonalisation of the Ising model can be found.

Preliminaries
In this paper, we will study the Hamiltonian of two spin chain models. The first is the spin-1/2 transverse field Ising model described by the quantum Hamiltonian: with L the number of spins in the chain and h proportional to the magnetic field oriented in the z direction and throughout this paper we always assume open boundary conditions. Operatorsσ x,y,z i are the Pauli spin operators of site i. As the parameter h approaches 1 (quantum critical point), the system undergoes a quantum phase transition between a ferromagnetic phase (h < 1) and a quantum paramagnetic phase (h > 1). In the paramagnetic phase the ground state is Z 2 invariant, and thus completely disordered, while for the ferromagnetic phase, in the thermodynamic limit, the ground state breaks Z 2 invariance (see for example [30]). This model is integrable and can be exactly diagonalised using the well known Jordan-Wigner and Bogoliubov transformations, as described in Appendix A.
The second model is the spin-1/2 XXZ model of a 1D chain with Hamiltonian: and anisotropy ∆. The phase diagram of this model is very well known: the ground state of Hamiltonian (2) is ferromagnetic for ∆ < −1 and Néel antiferromagnetic for ∆ > 1. In the intermediate region −1 < ∆ < 1 the system is in the critical-XY phase, characterised by absence of energy gap, divergence of the correlation length, logarithmic growth of the block entanglement entropy with the size of the block. We assume the system to be initially at time t = 0 in the ground state |Ψ 0 ofĤ 0 that is Hamiltonian Eq. (1) with magnetic field h 0 or Hamiltonian Eq. (2) with initial anisotropy ∆ 0 . The system is then evolved with a quench in the magnetic field or the anisotropy, respectively. We consider two scenarios. In the first one, the magnetic field h (anisotropy ∆) is instantaneously changed to a final value h 1 (∆ 1 ) and the system is let to evolve with the new HamiltonianĤ 1 that is Hamiltonian Eq. (1) with magnetic field h 1 or Hamiltonian Eq. (2) with final anisotropy ∆ 1 . In other words the state of the system at time t is: and we assumeh = 1 for simplicity.
In the second scenario, and only for the Ising model, we assume that the magnetic field changes linearly in time with rate γ: until the final time T = (h 0 − h 1 )/γ. Notice that γ > 0 when h 0 > h 1 and viceversa. In this case the state of the system follows Schrödinger equation: The entanglement spectrum of a state |Ψ is obtained by tracing out half of the chain and computing the eigenvalues of the reduced density matrix of one of the two halfs. In other words, we partition the chain in two parts Ω and Ω ⊥ containing the spins 1, 2, . . . , L/2 and L/2 + 1, . . . , L, respectively, assuming for simplicity even L. Then we compute the eigenvalues of the matrix: As Trρ Ω = 1, the sum of the eigenvalues λ 1 , λ 2 , . . ., sorted in decreasing order, is equal to 1. As shown in Ref. [18], for a large homogeneous chain, the exact position of the cut, as long as it is far from the chain edges, does not affect the results. In the next sections we will analyse the two models separately.

Quenches in the Ising model
In this section we discuss the results for the quenches in the Ising Hamiltonian. The details of the calculation can be found in the Appendix. We first consider instantaneous quenches and then move to quasi-adiabatic quenches.

Results for the Ising model: instantaneous quenches
In this section we discuss the results obtained for the ES after an instantaneous quench h 0 → h 1 for L = 200. The first quenches we consider are from the paramagnetic to the ferromagnetic phase, with h 0 = 1.5 and h 1 = 1.2, 1.0, 0.8, 0.5. As show in Fig. 1 the ES dynamics strongly depends on whether the system is quenched within the same phase (paramagnetic in this case) or in the other. On one hand, when the final Hamiltonian magnetic field h 1 ≥ 1 pertains to the paramagnetic phase, the Schmidt gap decreases almost monotonically (insets a-b). On the other hand when h 1 < 1, the Schmidt gap oscillates, with a series of zeros (insets c-d), inside the ferromagnetic phase. In the latter case, the ES exhibits, a series of crossings in time between the eigenvalues λ k (t). We observe that the frequency of the crossings increases lowering h 1 and going further deep in the ferromagnetic phase while it vanishes as h 1 → 1 at the quantum critical point. Such behaviour follows from the single particle eigenvalues ε k ofK, (for a definition see Eq B.14 in Appendix A) since the smallest one (which is the most important due to the exponential form of the reduced density matrix) goes periodically to zero, as plotted in Fig. 2. As we see, the Schmidt gap (black dashed curve) matches perfectly the oscillating behaviour of the lowest eigenvalue ε 1 .
To better understand the nature of these crossings we can look at the time t crossing when the first crossing between the two largest eigenvalues occurs (i.e. the first zero of the Schmidt gap), as a function of the final Hamiltonian parameter h 1 . As expected  the crossing time diverges algebraically as the post-quench regime gets closer to the quantum critical point. The result is shown in Fig. 3 where the numerical data for the crossing times are plotted against 1 − h 1 and compared to the simple scaling function: Notice that one would be tempted to associate a critical exponent to the power 1/2 that we find. However this is not the case because as we are quenching the system across the phase transition the evolved quantum state of the system is quite far from equilibrium and thus a description in terms of quasi-equilibrium critical exponents should not be adequate.
We would like to highlight a very close resemblance of the results for the Schmidt gap and those for the spontaneous magnetisation σ x (t) for a quench across the transition h 0 → h 1 [31,32]. These authors predict and observe numerically that the   crossing times go as: demonstrating that the first oscillation time t * for the magnetisation diverges indeed with the same scaling as the crossing time t crossing (h 1 ) for the Schmidt gap of Eq. (7). We emphasise, however, that, although the scaling of the time dependence of the two quantities, Schmidt gap and magnetisation is similar, here we are comparing results for quenches in the two opposite directions. We now turn to the analysis of the effects of the finiteness of the spin chain we consider. For short times we have not observed significant variation with the size of the system for sufficiently long chains. However for long times we observe revivals similar to the ones observed for the entanglement entropy. For the quench from h 0 = 1.5 to the quantum critical point h 1 = 1 the results are plotted in Fig. 4 for the six biggest eigenvalues λ k : revivals occur at the time t = L. Since the initial state is the ground state of the HamiltonianĤ 0 , when we quench the Hamiltonian, it will have a higher energy with respect to the new ground state of the Hamiltonian, and will act as a source of quasi-particle excitations along the chain [33]. The excitations travel along the chain with velocity v ≤ 1 (from the dispersion relation) and when they reach the end of the chain they are reflected backward, giving rise to such revivals at time t = L. This is the time for a quasi-particle emitted in the middle of the chain to reach the hard wall at the end of the chain, be reflected and come back at the centre.
From the large time behaviour of the entanglement spectrum we can also infer some properties of the equilibration ensemble. The Ising model should equilibrate to a Generalized Gibbs states. Generalized Gibbs states are exponential of extended Hamiltonians, that include, apart from the system Hamiltonian, also other operators encoding the local density of the conserved quantities. They thus give rise to extended phase diagrams that depend on several parameters. The values of the parameters in the generalized Hamiltonian are fixed by the values of the conserved quantities in the initial state. It has been conjectured that these generalized Hamiltonians are responsible for the appearance of dynamical quantum phase transitions [34]. This kind of transitions emerges when the combination of the initial state and quench Hamiltonian produce a generalized gapless Hamiltonian.
Our data in Fig. 2 seem to confirm such conjecture. At least in this simple case, when the system is quenched across a critical point, the generalized Hamiltonian (that is encoded in the logarithm of the reduced density matrix for large enough times), becomes gapless ( Fig. 1 panels b, c, d) while it displays a well defined gap when the system is quenched inside the same phase Fig. 1 panel a.  So far we have considered quenches from the paramagnetic phase to the ferromagnetic one and we have seen in the case of a quench connecting two different phases the occurring of crossing between the components of the ES. Such characteristic is surprisingly not found in the opposite case, where the initial state belongs to the ferromagnetic phase. As shown in Fig. 5 in both cases the degeneracy of the ground state (due to the Z 2 symmetry of the Hamiltonian in the ferromagnetic phase) is not lifted in the dynamics, and no crossing occurs. We have checked that the absence of crossing is not a consequence of the unbroken symmetry initial state. Using tDMRG and starting with a broken symmetry state with net magnetisation, we computed the ES after similar quenches. The results are reported in the lower panels of Fig. 5. As in the case with no symmetry breaking, the two largest eigenvalues never cross (although some crossing are observed for lower eigenvalues).
Despite the lack of crossings, even for the quench from the paramagnetic phase to the ferromagnetic phase, we appreciate that the generalized Hamiltonian describing the long time Generalized Gibbs State is gapped for quenches in the same phase and gapless for quenches across the phase transition.
In order to confirm that the entanglement spectrum becomes gapless once one crosses a phase transition, we have computed the entanglement spectrum of a small block that should be described at large times by the corresponding diagonal ensemble (DE) [27] ‡. Let us consider thus a small block, and compute the entanglement spectrum of a block of L/10 spins in a chain of L spins. This can be characterized by plotting the corresponding L/5 eigenvalues m k of the covariance matrix defined in Eq. (B.17), k = 1, . . . , L/5 versus k. The resulting plots are shown in Fig. 6. Any element of the entanglement spectrum is obtained by multiplying together L/10 of those m k .
The presence or absence of a gap between the first two eigenvalues of the reduced density matrix is related to the presence or absence of a gap in the plots of Fig. 6 close to L/10. A gap there means that the two largest eigenvalues are different (since they are obtained by choosing in the product one or the other, m L/10 or m L/10+1 ).
Still, first and second eigenvalues could be degenerate and a gap might exist between the second and the third. This is the typical case for Z 2 invariant theories where the eigenvalues of the reduced density matrix come in degenerate pairs. In this case the gap is now between m L/10−1 or m L/10+2 . Both cases are observed for quenches in the same phase. Gapped spectra are obtained for quenches in the paramagnetic phase while gapped degenerate spectra are obtained for quenches in the ferromagnetic phase as reported in the two panels Fig.  6. As soon as one studies quenches between two different phases, we observe that the entanglement spectrum becomes really gapless as shown by both panels of Fig. 6. ‡ We have checked that the same qualitative behavior is also observed in the half chain entanglement spectrum obtained from the same DE but in that case the DE is not expected to reproduce the correct quantitative behaviour of half spin chain since it should correctly capture the local physics only.

Results for the Ising model: adiabatic quench
Here we move to linear quasi adiabatic quenches in which the magnetic field is decreased linearly from h 0 to h 1 as in Eq. (4). Many previous works have concentrated on the dynamics of entanglement entropy and its relation to the Kibble-Zurek mechanism [35,36,37,38]. Our results of the ES dynamics for a quench from the paramagnetic phase h 0 = 1.5 to the ferromagnetic phase h 1 = 0 with rate γ = 10 −2 are shown in Fig. 7.
The system is initially in the ground state ofĤ 0 thus it is very close to a product state with all the spins pointing in the z direction. Thus initially λ 1 ≈ 1 while all the other eigenvalues are negligible. As the magnetic field decreases the first eigenvalue λ 1 decreases while all the others increase steadily until the instantaneous field h(t) passes the critical point. After this time the eigenvalues start oscillating and crossing (see Fig.7) in an analogous fashion to the instantaneous quenches. Eventually for long times and far below the critical point, the eigenvalues seem to tend to an asymptotic limit with residual oscillations as expected invoking the assumption of the generalised Gibbs ensemble. The eigenvalues distribution gives correctly the value for the von Neumann entropy predicted for slow quench (not shown).
We can see from Fig. 7 that the first crossing does not occur for h(t) = 1 but slightly later in time. The reason for this is twofold: first, close to the phase transition the reaction times of the system compete with the rate of change γ of the magnetic field, so there is a finite delay; second, for a finite system, the Schmidt gap of the ground state of the Ising model at the critical point is not zero but approaches zero with an algebraic decay § ruled by the ratio of critical exponents β/ν = 1/8 as found in [18,19].
These two effects combined together give rise to an interesting behaviour of the Schmidt gap ∆λ when the instantaneous field h(t) = 1 is at the critical point. The results for ∆λ at the critical point versus the rate of change γ of the magnetic field is shown in Fig. 8

Instantaneous quenches in the XXZ model
In this section we turn our attention to the XXZ model. The basic properties of the phase diagram have been described in Sec. 2. In contrast to the Ising model that is characterised by an isolated critical point, it exhibits an entire critical region |∆| < 1. The dynamics of the von Neumann entropy of the evolved state after a quench in the anisotropy from a value ∆ 0 to a value ∆ 1 has been first studied in Ref. [39]. Here we extend that analysis to the ES for quenches starting in a perfect Néel antiferromagnetic state. To calculate the ES, we distinguish two cases. For ∆ 1 = 0 the final Hamiltonian is the XX one and can be mapped into free fermions as for the Ising model and the ES can be computed analogously (see [39]). For ∆ 1 = 0 we employ the tDMRG method [29] to compute the time evolution of the initial state with HamiltonianĤ(∆ 1 ). In this section we fix L = 40. We consider values of ∆ 1 in or at the boundary of the critical region. The results are shown in Fig. 9. As for the Ising model the first eigenvalue starts from approximately 1 and decreases, while all other lower eigenvalues increase. This homogenisation of the eigenvalues gives rise to the linear increase of the von Neumann entropy. We do not observe dramatic crossings as for the Ising model. Instead the eigenvalues, after a fast transient, tend to decay together until a saturation. In fact, around t * ≈ 4, also the von Neumann entropy saturates because of the finiteness of the chain size. Moreover for ∆ 1 = 0 the eigenvalues become degenerate and decay approximately exponentially with time. For long times, it is possible to show that the covariance matrix M becomes a multiple of the identity matrix leading to a maximally mixed reduced density matrix after any bipartition.
In the case of the XXZ model, our data suggest the appearance of a gapless generalized Hamiltonian for quenches from the anti ferromagnetic phase to the boundary of the critical phase (∆ = 1) and in the critical phase (∆ > −1). This again suggests the appearance of a dynamical quantum phase transition for such quenches. Quenching in the ferromagnetic phase (∆ < −1), on the other hand seems to produce a gapped generalized Hamiltonian (Fig. 10). This suggests the absence of a dynamical phase transition despite our quench crossing two static quantum phase transitions (one at ∆ = 1 and another one at ∆ = −1).

Conclusions
In summary, we have presented an analysis of the entanglement spectrum dynamics in 1D spin chains. We have shown results for instantaneous quenches for the Ising and XXZ models. For the Ising model, we have shown how the dynamics of the Schmidt gap resembles very closely that of the spontaneous magnetisation, thus establishing a further connection with order parameters. Moreover, we have studied quasi-adiabatic quenches, in which the Schmidt gap scales algebraically with the rate of change of the magnetic field.
Our results demonstrate once more the ability of the ES to reveal interesting features of many-body system close to critical points. By analyzing the presence/absence of a gap in the long time evolution of the entanglement spectrum we have confirmed that in the simplest models dynamical quantum phase transitions seem to arise when quenching across a quantum critical point. More interestingly, in the XXZ model, we have observed the absence of a dynamical phase transition in a quench protocol across several quantum phase transitions. This result seems to confirm the conjecture of Ref. [34] that identifies the generalized Hamiltonian, rather than the system Hamiltonian as the one inducing dynamical phase transitions.
Many open questions remain that go beyond the scope of this work. For example it is natural to ask whether the ES dynamics across other Ising-type transitions, like the Néel-Haldane transition in spin-1 chains, behave similarly to the calculations showed here.
During the preparation of this manuscript we became aware of a related and independent work on entanglement spectrum dynamics in the Ising model [40].
first step is to map the Pauli operators for each site to a set of fermionic operators: This mapping transforms the Hamiltonian in a system of free spinless fermions: where the matrices A and B for the Ising model (1) read Notice that the real matrices A and B are symmetric and anti-symmetric, respectively. Following [41] we The single particle energies Λ kk ≥ 0 and the transformation matrices can be directly calculated by the singular-value decomposition of the sum of the matrices A and B where the column of the matrices φ = Π + Γ and ψ = Π − Γ are the eigenvectors of H. Due to the quadratic form of the Hamiltonian, the system is fully described by the one-particle Green's function G(c, c † ) Appendix B. Calculation of the entanglement spectrum dynamics for the Ising model Let us first consider the dynamics of the system after an instantaneous quantum quench on the magnetic field. The initial Hamiltonian at time t < 0 is abruptly changed tô for t ≥ 0. All the properties of the system, including the ES, can be computed directly from the time-dependent Green's function G(c, c † ; t). Within the Heisenberg picture, the time evolution is governed by the evolution operatorÛ (t, 0) = exp{−i L k=1 (Λ 1 ) kkξ † kξ k t}. It is then convenient to first transform the operators (c, c † ) into the new set (ξ, ξ † ), for which the time evolution results in a simple phase factor The operators ξ(t) can then be mapped back into the set (η, η † ), where the Green's function is trivial since the initial state |Ψ 0 is the vacuum of Hamiltonian (B.1) and therefore η |Ψ 0 = 0. If B 0 and B 1 are respectively the Bogoliubov transformation that maps the operators (c, c † ) into (η, η † ) and (ξ, ξ † ), such chain of transformations leads to where D(t) is defined as Defining the following matrices it is easy to show that the time-dependent Green's function is where we introduced the correlation matrices with C(t), F (t) ∈ Ω ⊕ Ω ⊥ . We now divide the full system into subsystem Ω and its orthogonal complement Ω ⊥ . Since we consider a fermionic system, the eigenstates of the total system Ω ⊕ Ω ⊥ are Slater determinants and thus, according to Wick theorem, the reduced density matrix for the subsystem Ω has exponential form [43] ρ Ω = Tr Ω ⊥ |Ψ Ψ| = Z exp(−K) (B.12) withK quadratic in the Fermi operatorŝ where the coefficients α ij and β ij are related to the correlation matrices C(t) and F (t).
Following an analogous procedure employed for diagonalising (A.2), the entanglement HamiltonianK can be brought to its diagonal form K = L k=1 ε kη † kη k (B.14) with the appropriate transformationB. Since the system is completely characterized by the one-particle Green's function, we can derive the time evolution of the ES λ(t) of ρ Ω (t) directly from the projection onto Ω of the correlation matrices C and F c(t) = P Ω C(t) P Ω f (t) = P Ω F (t) P Ω c(t), f (t) ∈ Ω. (B.15) If the partition is in real space as in our case, then the projection (B.15) consists in selecting the elements in C(t) and F (t) with indexes belonging to the subsystem Ω. Such correlation functions are obtained directly from Eqs. (B.10)-(B.11) but the same results must be found from which provides us the way to relate the eigenvalues ε k ofK (and thus the ES) with the Green's functions components [44,45]. The eigenvalues ofK are found from the diagonalisation of the L × L matrix M M (t) = 2c(t) − I 2f (t) −2f (t) * −(2c(t) * − I) (B.17) whose eigenvalues are m k (t) = ± tanh(ε k (t)/2).
It is important to remind here that once the m k are known the actual entanglement spectrum is obtained by a product of L of them. The presence or absence of a gap in the entanglement spectrum is then directly related to the presence or absence of degeneracies in the eigenvalues of M as discussed in the main text. The reduced density matrix is then whose eigenvalues are simply: where the numbers n i k = 0, 1 and the normalisation constant Z is found by imposing Nmax i=1 λ i = 1 with N max = 2 L/2 . As it is practically impossible to calculate the normalisation constant exactly as soon as L becomes large (we use 32 ≤ L ≤ 512), we normalise the eigenvalues restricting the sum appearing in (B.19) to 1 ≤ k ≤ k max . This limits the calculation to the first 2 kmax ≤ N max eigenvalues. In our simulations we observe that k max ≤ 8 gives accurate results to machine precision.
For non instantaneous quenches the HamiltonianĤ(t) is changing with time. This means that the fermionic operators of the instantaneous HamiltonianĤ(t) also change with time. In order to deal with this slightly more complicated case, one can use a timedependent Bogoliubov transformation as explained for example in Ref. [35] where a set of ordinary differential equations in time must be solved numerically. We avoid this step by using an explicit finite difference expansion of the evolution. In other words, in the Heisenberg picture, we evolve the Green's function G(c, c † ; t) for a small time interval dt to obtain G(c, c † ; t + dt) using the same procedure of the instantaneous quenches case assuming the HamiltonianĤ(t + dt) to be constant in the interval [t; t + dt]. In the numerical simulations we used dt = 10 −3 which is sufficiently small to obtain accurate results.