We consider the non-equilibrium dynamics after a sudden quench of the magnetic field in the transverse field Ising chain starting from excited states of the pre-quench Hamiltonian. We prove that stationary values of local correlation functions can be described by the generalised Gibbs ensemble (GGE). Then we study the full time evolution of the transverse magnetisation by means of stationary phase methods. The equal time two-point longitudinal correlation function is analytically derived for a particular class of excited states for quenches within the ferromagnetic phase, and studied numerically in general. The full time dependence of the entanglement entropy of a block of spins is also obtained analytically for the same class of states and for arbitrary quenches.
Quantum quenches from excited states
in the Ising chain
Leda Bucciantini, Márton Kormos, Pasquale Calabrese
Dipartimento di Fisica dell’Università di Pisa and INFN, 56127 Pisa, Italy
MTA-BME “Momentum” Statistical Field Theory Research Group,
1111 Budapest, Budafoki út 8, Hungary
Recent experiments in the field of ultra-cold atoms [2, 3, 4, 5, 6] triggered an intense theoretical activity aimed to understand the unitary non-equilibrium evolution of isolated many-body quantum systems. A situation which attracted a lot of interest is that of a sudden quench of a Hamiltonian parameter such as an external magnetic field . The most remarkable results that emerged from these theoretical and experimental investigations are probably the following two: (i) there is a light-cone like spreading of correlations following a quantum quench [4, 8, 9, 10, 11] and (ii) expectation values of local observables generically approach stationary values at late times in the thermodynamic limit, although the whole system is always in a pure state.
Rather amazingly, these stationary values can be predicted by a statistical ensemble without solving the complicated non-equilibrium dynamics. For non-integrable models the appropriate statistical ensemble is expected to be the standard Gibbs one with an effective temperature fixed by the value of the energy in the initial state . For integrable models, the existence of non-trivial local conservation laws strongly constrains the dynamics and the stationary values are expected to be described by a generalised Gibbs ensemble (GGE) . The most convincing evidence supporting this scenario comes from the exact solution of models that can be mapped to free fermions or bosons for which the full dynamics can be obtained analytically. Among these models a crucial role has been played by the transverse field Ising chain which is probably the most intensively investigated model in the quench literature [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31], although many other free-like models have been studied as well. [32, 33, 34, 35, 36, 37, 38, 39, 40]. For integrable interacting models, i.e. with a non-trivial scattering matrix between quasi-particle excitations, the exact solution of the quench dynamics is still an outstanding problem , but recently it has been possible in some instances to construct the GGE and compute explicitly a few observables [41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57] and also to check these predictions numerically . Clearly, for non-integrable models the evidence comes only from numerical [58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71] or experimental [2, 3] investigations.
However, most of (if not all) these studies focused on the evolution starting from the ground-state of a given pre-quench Hamiltonian. While these initial states should include all the most relevant experimental situations, it is natural to wonder how general the conclusions drawn on their basis are. Indeed, ground states of local Hamiltonians are not at all generic because their entanglement entropy (i.e. the von Neumann entropy of the reduced density matrix of a subsystem) satisfies an area law  (i.e. it scales with the area of the surface of the subsystem) or at most it has multiplicative logarithmic corrections, while generic states satisfy a volume law. This important difference is at the heart of the machinery of the so-called tensor network techniques used to effectively simulate many-body quantum systems (see e.g. Refs.  as reviews). Somewhat more general conclusions can be drawn by considering a quantum quench in which the initial state is not a ground-state, but an excited state of the pre-quench Hamiltonian, which generically follows a volume law for the entanglement entropy  rather than an area law.
In order to give a first answer to this question, in this paper we consider the quench dynamics in the prototype of the models mappable to free fermions, namely the transverse field Ising chain, which, in spite of its simplicity and the fact that it is exactly solvable, represents a crucial paradigm for quantum critical behaviour. Furthermore, while the model admits a representation in terms of free fermions, the spin variables are non-local with respect to the fermionic degrees of freedom, a property which renders it less trivial than free particle systems and, at the same time, an ideal testing ground for studying the relaxation for which locality is an essential feature.
The manuscript is organised as follows. In Sec. 2 we report all the preliminary information to set up the calculations, i.e. we introduce the model, its diagonalisation, the quench protocol and the observables we will study. In Sec. 3 we show that in the long time limit (taken after the thermodynamic one) the system can be described by a GGE. In Sec. 4 we study the time evolution of the transverse magnetisation, while in the following Sec. 5 we consider the two-point longitudinal correlation function. In Sec. 6 we turn to the time evolution of the entanglement entropy. In the last section we draw our conclusions and in two appendices we report some additional technical details.
2 The model, the quench protocol, and the observables
We consider here the non-equilibrium dynamics of the transverse field Ising chain with Hamiltonian
where , is the Pauli matrix at site of the chain, is the transverse field, and we impose periodic boundary conditions . The transverse field Ising chain is a crucial paradigm of quantum phase transitions  because at zero temperature and in the thermodynamic limit it exhibits ferromagnetic () and paramagnetic () phases, separated by a quantum critical point at .
2.1 The diagonalisation of the Hamiltonian
The diagonalisation of Hamiltonian (1) is presented in several textbooks (see e.g. ), but we repeat it here to make the paper self consistent. The first step in diagonalising the Hamiltonian (1) is to introduce a set of spinless fermion annihilation and creation operators through the non-local Jordan-Wigner transformation
where the operators , satisfy anti-commutation relation
Then the Fourier modes are defined as
Since the transformation is unitary, also the operators satisfy anti-commutation relations. In terms of operators, the Hamiltonian is
where the sum over the modes runs over integers or half-integers depending on the parity of the fermion number, see e.g. the Appendix of Ref.  for a detailed discussion. This is however superfluous for our goals since we will only consider the thermodynamic limit in which the momentum becomes continuous.
It is now necessary to apply one more unitary transformation to cast the Hamiltonian in a diagonal form: the transformation is a Bogoliubov rotation which takes the fermionic operators , to the new ones ,
where the coefficients and are chosen in such a way as to make diagonal
and the angle is defined by the relation
Again, the unitarity of the transformation ensures the validity of the usual anti-commutation relation for the operators. Moreover, the characteristic of (6) of mixing only the modes and of the operators allows us to re-write the Bogoliubov transformation in a more compact way as a rotation in a Hilbert space
where and are the two-component vectors
and the matrix is a special case of
The Hamiltonian can be written in terms of the Bogoliubov quasi-particles as
where the one-particle dispersion relation is
In the thermodynamic limit, the momentum becomes a continuos variable living in the interval . In the following we will often use instead of , being sure that it will always be clear whether we refer to continuous or discrete momenta.
Other useful fermionic operators
For the calculation of correlation functions it is convenient to introduce some other sets of fermionic operators. First, we introduce the Majorana fermions
which satisfy the algebra
While the spin operator is local in terms of these Majorana fermions, , the operator has the non-local representation
which, as we shall see, is particularly useful in the calculation of real space correlation functions and the entanglement entropy.
The operators and can be collected together in two different ways. On the one hand, one can introduce a single set of operators at the price of doubling their number per site in the following way
and they satisfy the algebra
On the other hand, and can be collected in a two-component vector operator
We will denote the Fourier transform of with and the Fourier transform of with , and, with a slight abuse of notation, we will denote the Fourier transform of the vector (19) simply as :
2.2 The exact spectrum
The ground state of is the vacuum of Bogoliubov operators, i.e. it is annihilated by all . Its energy is . However, the exact diagonalisation of the model gives not only the ground state properties but all the eigenstates and their energies. In the basis of free fermions, the excited states are classified according to the occupation numbers of the single-particle basis. An eigenstate can be then written as
where is a characteristic function of the state representing the set of occupied momenta, i.e. if the momentum is occupied and if not. All observables can be written in terms of the characteristic function as for example the energy in Eq. (21). While in finite system the characteristic function can assume only the values and , in the thermodynamic limit it becomes an arbitrary function of with the restriction to be in the interval . This function is a coarse-grained version of , see e.g.  for specific examples. The ground state corresponds to for every .
2.3 The quench protocol
The time dependence of the system after a quench of the transverse field starting from the ground state has been studied extensively in a series of works, as reported in the introduction. Here we are interested in the case when the initial state is an excited state of the pre-quench Hamiltonian, i.e. for the system is in an excited state of the Hamiltonian with field equal to At time the value of the field is suddenly quenched to and all the following time evolution is governed by this new Hamiltonian. We will denote by and the fermionic mode operators that diagonalise the Hamiltonian with and , respectively. Similarly, primed symbols will be used to denote all pre-quenched operators and variables while not-primed ones for post-quench operators and variables. The initial state is then one of those in Eq. (21) for the pre-quench Hamiltonian, i.e.
As already discussed this state is fully specified by the characteristic function (in finite systems). The time dependent state is given by
with being the post-quench Hamiltonian with transverse field .
We point out that even if the states are a basis for the many-body Hilbert space, they do not represent the most generic excited state because the spectrum of the Ising chain is highly degenerate and linear combination of degenerate states are still eigenstates, but they cannot be written as (21). One property of the states is that they do not break the symmetry of the Hamiltonian even in the ferromagnetic phase, unless i.e. in the ground state.
The relation between the pre- and post-quench Bogoliubov operators is given by the combination of the two Bogoliubov rotations in Eq. (6) with angles and (i.e. post- and pre-quench Bogoliubov angles, respectively). Thus the overall rotation is
where we defined . Being the main quench variable entering in all the following calculations and results, it is worth writing explicitly its form in terms of and
2.4 Observables, Relaxation and Generalised Gibbs ensemble
The most important observable in the study of quantum quenches is the reduced density matrix of a block built with contiguous spins. Indeed from all the correlation functions local within can be obtained. Since the system is in a pure state at any time, the density matrix of the entire system is
The reduced density matrix of a subsystem is defined as
where is the complement of . The importance of stems from the fact that it is the quantity which generically displays a stationary behaviour described by some statistical ensemble, while the full density matrix always corresponds to a pure state with zero entropy.
exists. This is described by a given statistical (mixed state) ensemble with full density matrix if its reduced density matrix restricted to equals , i.e. if for and for any finite subsystem
In particular, this implies that arbitrary local multi-point correlation functions within subsystem can be evaluated as averages within the . By no means this implies that equals the full density matrix of the system which is clearly impossible being the former a mixed state and the latter a pure one.
When a system thermalises, is the standard Gibbs distribution and this is expected to be the case when the model is non-integrable. However, for an integrable model, the proper statistical ensemble describing the system for long time is a generalised Gibbs ensemble (GGE) rather than a thermal one. The density matrix of the GGE is defined as 
where is set of commuting integrals of motion, i.e. , and is a normalisation constant . It is important for a proper definition of GGE to specify which conserved charges enter in the GGE density matrix above. It has been understood recently [22, 26] that only local integrals of motion should be used in Eq. (30) if we are interested in the expectation values of local observables such as the reduced density matrix.
While in general it is a formidable task to calculate a reduced density matrix even for an integrable system, in the case of a model that can be mapped to free fermions it is a rather straightforward application of the Wick theorem to write it in terms of only two-point correlators of fermions. To this aim, let us first introduce the correlation matrix of Majorana fermions with the definition
For consecutive fermions ( Majoranas), using explicitly the periodicity of the chain, this matrix has the block structure
where the ’s are matrices with entries equal to the correlations of Majorana fermions, which explicitly are
where the correlations can clearly be taken to start from an arbitrary site . Because of its periodic structure, the matrix turns out to be a block Toeplitz matrix, i.e. its constituent blocks depend only on the difference between row and column indices of
We can now use Wick’s theorem to construct all correlation functions in the Ising chain. As shown in Refs. [76, 77], the matrix determines entirely the reduced density matrix of the block of contiguous fermions in the chain (and hence spins, because contiguous spins are mapped to contiguous fermions, see Eq. (16)), with a final result that can be written in the compact way [76, 78]
Given we can calculate any local correlation function with support in . Very importantly, because of this direct relation between reduced density matrix and correlation matrix, it is sufficient to prove that the two ensembles or states have the same correlation matrix in order to prove that they are equal. This is an immense simplification because while has elements, has only elements.
From the reduced density matrix, another fundamental observable is easily constructed, namely the entanglement entropy which is the von Neumann entropy of ,
Using again Wick’s theorem , can be related to the eigenvalues of the matrix . Indeed, denoting the eigenvalues of as , (being an antisymmetric matrix, its eigenvalues are purely imaginary complex conjugate pairs), the entanglement entropy is 
Apart from the reduced density matrix, we will also consider the transverse magnetisation
and the two-point function of the order parameter at a distance
where is given by
Notice that although the matrices and in Eq. (33) look very similar and they have the same block-diagonal elements, the off-diagonal ones are different since the second operator is shifted by . We note that in the literature the two matrices are often denoted by the same symbol and it is very easy to mix them up.
3 The infinite time limit and the generalised Gibbs ensemble
In this section we consider the infinite time limit of the reduced density matrix of a subsystem A composed of contiguous spins. To analyse , we first consider the time evolution of the long time limit of its building blocks, that, according to Eq. (34), are the two-point real-space correlation functions of fermions.
3.1 Time evolution of the fermionic two-point function
The Bogoliubov rotations diagonalising the pre- and post-quench Hamiltonians only couple modes with opposite momenta, cf. Eq. (9). It is then convenient to cast the two-point correlation functions of pre-quench Bogoliubov modes in the matrix
in which is the initial state specified by the function as in Eq. (22). Combining the two Bogoliubov rotations for pre- and post-quench Hamiltonian as in Eq. (24), we write the expectation value of post-quench Bogoliubov operators in the initial state as
which is the initial condition for the fermionic two-point functions. The time evolution can be worked out in the Heisenberg picture where the post-quench operators evolve according to the Hamiltonian (12), so , where is the restriction of the time evolution operator to the subset of the Hilbert space with momenta and , i.e.
It is now possible to evaluate the expectation value of expressions bilinear in the fermions and at any time. In order to do so, it is enough to invert Eq. (4) and express the operators in terms of which evolve according to Eq. (46). From , we can write
and similarly for . Hence from Eq. (20) we obtain for the Fourier transform of Majorana operators
The vector operator has the form
where we defined
which stand for the even and odd part of respectively. It is straightforward to check that when and are set to zero this expression reduces to the one obtained in the case when the initial state is the ground state of the initial Hamiltonian [17, 26]. Note that if then Eq. (50) is constant in time, which is a manifestation of the fact that if for any the state is not only an eigenstate of the pre-quench Hamiltonian but also of the post-quench one. We will define the states with as parity invariant states (PIS) in which all the positive and negative momentum modes are populated with the same weights. Note that all the PIS have zero momentum, but the condition for PIS is more restrictive than that. We will refer to the states with for some as non parity invariant states (NPIS).
Equation (50) is the final expression for the two-point function of fermions in momentum space from which, by Fourier transform, the one for real-space fermions given by Eq. (33) can be straightforwardly obtained. As a usual feature of free systems, each momentum mode oscillates in time with typical frequency proportional to . However, when taking the Fourier transform, in the thermodynamic limit the various modes interfere in a destructive way and their long-time expectation is the time-average of the expression above, i.e.
Thus, in order to show that the reduced density matrix attains a stationary behaviour described by GGE, it is sufficient to show that the GGE prediction for equals Eq. (53). By no means this implies that has a long-time limit, on the contrary, it oscillates forever as a consequence of the fact that the state is pure for any time and the Hamiltonian governing the evolution is diagonal in the modes.
3.2 GGE expectation value of the fermionic two-point function
The GGE density matrix for the full system is given by Eq. (30) constructed with the local integrals of motion. However, it has been shown that for the transverse field Ising chain the post-quench occupation number operators
although non-local quantities, can be written as linear combinations of the local integrals of motion  (see also the next subsection). Thus the GGE density matrix constructed with local integrals of motion and the one constructed with are equivalent. We will then consider
where we use the same symbols for the Lagrange multipliers and in Eq. (30) since we will never use the two concomitantly. The are fixed by matching the expectation values of the occupation numbers with their values in the initial state, i.e. imposing
The left-hand-side of this equation can be read out from Eq. (45)
while the right-hand-side is
Equating the two expressions we get the equation determining :
The components of can be readily calculated in the GGE, for example
Performing similar calculations for the other three elements of the matrix we finally get
which coincides with Eq. (53). This proves that the GGE two-point functions of fermions at arbitrary distance are equal to the long-time limit of the same two-point function after a quench from an excited state of the initial Hamiltonian. Since the reduced density matrix can be constructed solely from the fermionic two-point functions as in Eq. (34), this also proves that any local multipoint correlation function of spins and fermions will be described by the GGE for long times.
3.3 Local Conservation laws in the TFIC
It is instructive to have a look at the behaviour of the local conservation laws in the chain in order to check whether they bring any further understanding in the non-equilibrium quench dynamics. In the thermodynamic limit there is an infinite number of local conserved charges which can be written in terms of the post-quench occupation numbers as [79, 28]
where the apex refers to their parity properties: are even and are odd under spatial reflections. The expectation values of these conserved charges in the initial state (and in the subsequent time evolution) is obtained by inserting the expectation value of of Eq. (57) into Eq. (62):
Hence and are both, in general, non-vanishing for an initial excited state. This is different from what happens if the initial state is the ground state of the pre-quench Hamiltonian, when all the parity-odd charges vanish. However, is zero every time that , i.e. for PIS with . This will represent an important class of states for which the calculation of some local observables in the following will be much easier. It is surely interesting to understand whether the increased number of non-zero conservation laws alters somehow the time-dependence of the asymptotic behaviour of local correlation functions.
4 Transverse magnetisation
The first observable we consider is the transverse magnetisation which is particularly easy to calculate because it has a local expression in terms of fermions, cf. Eq. (38).
For an excited state with characteristic function , can be easily found expressing the in terms of the and substituting the matrix elements of Eq. (50). After simple algebra we obtain in the thermodynamic limit
which again reduces to the ground-state evolution in the case , i.e. . Quite remarkably, we have that only the symmetric part of the characteristic function contributes to the time evolution of the transverse magnetisation and so the odd conserved charges in Eq. (63) have no influence at all on this observable.
Furthermore, the result can be divided into a stationary and a time-dependent part. As for the ground-state case, it is particularly interesting to understand the approach to the stationary value which can be evaluated by a stationary phase approximation. The stationary points are the zeros of in the interval , which are . However, when the characteristic function is not an analytic function in , possible new extrema of the integration domain need to be taken into account in the calculation. Indeed, these extrema need not coincide with , as it happens for example in the case in which the initial state excitations are non-zero only in a particular subinterval of .
In the general case under consideration, we have integrals of the form
where and do not need to coincide with and . Notice also that in Eq. (65) the stationary points of are always zeros of , requiring to go to the second order in stationary phase approach. There are three different classes of extremal points which should be summed up in order to have the complete large time behaviour. Denoting with each of the points, the three classes are
The extremal point is internal to the domain of integration, i.e. . In this case must be stationary to be extremal, i.e. , and then we have
where . Because of the leading behaviour is instead of .
The extremal point is at the boundary of the integration domain, i.e. or and it is also stationary, i.e. , in which case we have
Finally, the extremal point can be at the boundary of the integration domain, i.e. or , but it is not stationary i.e. . For and not stationary we have
and the same without the minus sign for .
It is clear that cases (1) and (2) give a different contribution compared to case (3). Indeed the first two cases generically give a large time behaviour of the type (whenever , else the contribution would be faster like ), like in the quench from the ground state. The third case instead will generically produce a slower decay of the type (unless , in which case we will have ).
For a generic initial state with characteristic function , the general strategy would be the following: (i) divide the integral in Eq. (65) in pieces in which is continuous, (ii) treat each of the integrals as in Eq. (66), and finally (iii) sum up all the contributions with the slowest power-law behaviour. In order to show the differences between the various states, we compare the evolution from the ground state with the ones from three representative excited states. We choose the following three states
where all are clearly defined in the interval .
For the initial ground state we regain the large time behaviour 
For the initial state characterised by we obtain
The leading large time behaviour going like comes from the extrema (non-stationary points) while the from the stationary point. Conversely, for the initial states characterised by and we get a power law behaviour due to the stationary points falling into cases and above, similarly to the ground state but with different calculable coefficients that we do not report here, but are easily obtained from Eqs. (67) and (68).
In order to show the reliability and the range of validity of the stationary phase approximations, in Fig. 1 we report the time dependent part of the transverse magnetisation (i.e. we subtract the stationary behaviour). We compare the exact results from the numerical determination of the integral in Eq. (65) with the stationary phase approximation up to order . It is evident that even for not so-large time the oscillating power-law decay of the stationary phase correctly describes the data. We mention that in the case of Eq. (73) it is important to keep the term to describe the data for not too large times.
To conclude this section we would like to emphasise the main difference we have found in the large time behaviour of the transverse magnetisation starting from the ground-state or from excited states of the pre-quench Hamiltonian. For quenches starting from the ground-state we always have a power-law tail of the form . While several excited states have the same power-law behaviour, this is not true in general. The state with above presents a much slower relaxation going like . We stress that this is not at all an academic state, quite the reverse, it is the most physical among the ones presented above because it has all the modes larger than a given one (, but this is not essential) occupied. Furthermore by choosing very particular momentum occupation functions , it is not difficult to cook up quite untypical power-law behaviours: for example, considering , since it vanishes in all the stationary points of the phase , we obtain a power-law decay like .
5 Equal-time two point longitudinal correlation function
In this section we investigate the longitudinal spin-spin correlation function between two spins at a distance at the same time , i.e.
The most interesting regime of this two-point function is the so-called space-time scaling limit [22, 25] defined as the limit with their ratio kept fixed. In general, the space-time scaling limit does not have to commute with the limit or either if taken before or after the thermodynamic limit.
The two-point function is the Pfaffian of the matrix in Eq. (42) the elements of which are the already calculated two-point fermion functions in Eq. (43). The fermionic correlators in Eq. (43) can be identified looking at Eq. (50), obtaining that the two-by-two constituent blocks have the form
Thus in the thermodynamic limit we have