Abstract
We consider the nonequilibrium dynamics after a sudden quench of the magnetic field in the transverse field Ising chain starting from excited states of the prequench 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 twopoint 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
[2.3em]
Leda Bucciantini, Márton Kormos, Pasquale Calabrese
Dipartimento di Fisica dell’Università di Pisa and INFN, 56127 Pisa, Italy
MTABME “Momentum” Statistical Field Theory Research Group,
1111 Budapest, Budafoki út 8, Hungary
1 Introduction
Recent experiments in the field of ultracold atoms [2, 3, 4, 5, 6] triggered an intense theoretical activity aimed to understand the unitary nonequilibrium evolution of isolated manybody 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 [7]. The most remarkable results that emerged from these theoretical and experimental investigations are probably the following two: (i) there is a lightcone 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 nonequilibrium dynamics. For nonintegrable 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 [12]. For integrable models, the existence of nontrivial local conservation laws strongly constrains the dynamics and the stationary values are expected to be described by a generalised Gibbs ensemble (GGE) [13]. 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 freelike models have been studied as well. [32, 33, 34, 35, 36, 37, 38, 39, 40]. For integrable interacting models, i.e. with a nontrivial scattering matrix between quasiparticle excitations, the exact solution of the quench dynamics is still an outstanding problem [41], 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 [57]. Clearly, for nonintegrable 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 groundstate of a given prequench 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 [72] (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 socalled tensor network techniques used to effectively simulate manybody quantum systems (see e.g. Refs. [73] as reviews). Somewhat more general conclusions can be drawn by considering a quantum quench in which the initial state is not a groundstate, but an excited state of the prequench Hamiltonian, which generically follows a volume law for the entanglement entropy [74] 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 nonlocal 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 twopoint 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 nonequilibrium dynamics of the transverse field Ising chain with Hamiltonian
(1) 
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 [75] 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. [75]), 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 nonlocal JordanWigner transformation
(2) 
where the operators , satisfy anticommutation relation
(3) 
Then the Fourier modes are defined as
(4) 
Since the transformation is unitary, also the operators satisfy anticommutation relations. In terms of operators, the Hamiltonian is
(5) 
where the sum over the modes runs over integers or halfintegers depending on the parity of the fermion number, see e.g. the Appendix of Ref. [25] 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 ,
(6) 
where the coefficients and are chosen in such a way as to make diagonal
(7) 
and the angle is defined by the relation
(8) 
Again, the unitarity of the transformation ensures the validity of the usual anticommutation relation for the operators. Moreover, the characteristic of (6) of mixing only the modes and of the operators allows us to rewrite the Bogoliubov transformation in a more compact way as a rotation in a Hilbert space
(9) 
where and are the twocomponent vectors
(10) 
and the matrix is a special case of
(11) 
The Hamiltonian can be written in terms of the Bogoliubov quasiparticles as
(12) 
where the oneparticle dispersion relation is
(13) 
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
(14) 
which satisfy the algebra
(15) 
While the spin operator is local in terms of these Majorana fermions, , the operator has the nonlocal representation
(16) 
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
(17) 
and they satisfy the algebra
(18) 
On the other hand, and can be collected in a twocomponent vector operator
(19) 
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 :
(20) 
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 singleparticle basis. An eigenstate can be then written as
(21) 
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 coarsegrained version of , see e.g. [74] 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 prequench 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 prequenched operators and variables while notprimed ones for postquench operators and variables. The initial state is then one of those in Eq. (21) for the prequench Hamiltonian, i.e.
(22) 
As already discussed this state is fully specified by the characteristic function (in finite systems). The time dependent state is given by
(23) 
with being the postquench Hamiltonian with transverse field .
We point out that even if the states are a basis for the manybody 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 postquench Bogoliubov operators is given by the combination of the two Bogoliubov rotations in Eq. (6) with angles and (i.e. post and prequench Bogoliubov angles, respectively). Thus the overall rotation is
(24) 
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
(25) 
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
(26) 
The reduced density matrix of a subsystem is defined as
(27) 
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.
More precisely, following Refs. [33, 26] it is usually said that a system reaches a stationary state if a long time limit of the reduced density matrix exists, i.e. if the limit
(28) 
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
(29) 
In particular, this implies that arbitrary local multipoint 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 nonintegrable. 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 [13]
(30) 
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 twopoint correlators of fermions. To this aim, let us first introduce the correlation matrix of Majorana fermions with the definition
(31) 
For consecutive fermions ( Majoranas), using explicitly the periodicity of the chain, this matrix has the block structure
(32) 
where the ’s are matrices with entries equal to the correlations of Majorana fermions, which explicitly are
(33) 
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]
(34) 
where [78]
(35) 
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 ,
(36) 
Using again Wick’s theorem [76], 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 [76]
(37) 
Apart from the reduced density matrix, we will also consider the transverse magnetisation
(38) 
and the twopoint function of the order parameter at a distance
(39) 
While is local within fermions, is not. However, in computing only the string of Majorana fermions between sites and matters and thus it takes the form [14, 75]
(40) 
By means of Wick’s theorem [14, 75], can be written as the Pfaffian of a skew symmetric matrix
(41) 
where is given by
(42) 
where
(43) 
Notice that although the matrices and in Eq. (33) look very similar and they have the same blockdiagonal elements, the offdiagonal 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 twopoint realspace correlation functions of fermions.
3.1 Time evolution of the fermionic twopoint function
The Bogoliubov rotations diagonalising the pre and postquench Hamiltonians only couple modes with opposite momenta, cf. Eq. (9). It is then convenient to cast the twopoint correlation functions of prequench Bogoliubov modes in the matrix
(44) 
in which is the initial state specified by the function as in Eq. (22). Combining the two Bogoliubov rotations for pre and postquench Hamiltonian as in Eq. (24), we write the expectation value of postquench Bogoliubov operators in the initial state as
(45) 
which is the initial condition for the fermionic twopoint functions. The time evolution can be worked out in the Heisenberg picture where the postquench 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.
(46) 
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
(47) 
and similarly for . Hence from Eq. (20) we obtain for the Fourier transform of Majorana operators
(48) 
The vector operator has the form
(49) 
Thus
(50) 
where we defined
(51)  
(52) 
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 prequench Hamiltonian but also of the postquench 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 twopoint function of fermions in momentum space from which, by Fourier transform, the one for realspace 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 longtime expectation is the timeaverage of the expression above, i.e.
(53) 
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 longtime 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 twopoint 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 postquench occupation number operators
(54) 
although nonlocal quantities, can be written as linear combinations of the local integrals of motion [28] (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
(55) 
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
(56) 
The lefthandside of this equation can be read out from Eq. (45)
(57) 
while the righthandside is
(58) 
Equating the two expressions we get the equation determining :
(59) 
The components of can be readily calculated in the GGE, for example
(60)  
Performing similar calculations for the other three elements of the matrix we finally get
(61) 
which coincides with Eq. (53). This proves that the GGE twopoint functions of fermions at arbitrary distance are equal to the longtime limit of the same twopoint function after a quench from an excited state of the initial Hamiltonian. Since the reduced density matrix can be constructed solely from the fermionic twopoint 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 nonequilibrium quench dynamics. In the thermodynamic limit there is an infinite number of local conserved charges which can be written in terms of the postquench occupation numbers as [79, 28]
(62)  
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):
(63)  
Hence and are both, in general, nonvanishing for an initial excited state. This is different from what happens if the initial state is the ground state of the prequench Hamiltonian, when all the parityodd 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 nonzero conservation laws alters somehow the timedependence 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).
In the case when the initial state is the prequench ground state the transverse magnetisation is in the thermodynamic limit [14, 26]
(64) 
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
(65) 
which again reduces to the groundstate 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 timedependent part. As for the groundstate 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 nonzero only in a particular subinterval of .
In the general case under consideration, we have integrals of the form
(66) 
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
(67) 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
(68) 
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
(69) 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 powerlaw 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
(70) 
where all are clearly defined in the interval .
For the initial ground state we regain the large time behaviour [26]
(71) 
where
(72) 
For the initial state characterised by we obtain
(73) 
where
(74) 
The leading large time behaviour going like comes from the extrema (nonstationary 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 solarge time the oscillating powerlaw 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 groundstate or from excited states of the prequench Hamiltonian. For quenches starting from the groundstate we always have a powerlaw tail of the form . While several excited states have the same powerlaw 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 powerlaw behaviours: for example, considering , since it vanishes in all the stationary points of the phase , we obtain a powerlaw decay like .
5 Equaltime two point longitudinal correlation function
In this section we investigate the longitudinal spinspin correlation function between two spins at a distance at the same time , i.e.
(75) 
The most interesting regime of this twopoint function is the socalled spacetime scaling limit [22, 25] defined as the limit with their ratio kept fixed. In general, the spacetime scaling limit does not have to commute with the limit or either if taken before or after the thermodynamic limit.
The twopoint function is the Pfaffian of the matrix in Eq. (42) the elements of which are the already calculated twopoint fermion functions in Eq. (43). The fermionic correlators in Eq. (43) can be identified looking at Eq. (50), obtaining that the twobytwo constituent blocks have the form
(76) 
with elements
(77)  
(78)  
(79) 
Thus in the thermodynamic limit we have