A finite quantum system evolving unitarily equilibrates in a probabilistic fashion. In the general many-body setting the time-fluctuations of an observable are typically exponentially small in the system size. We consider here quasi-free Fermi systems where the Hamiltonian and observables are quadratic in the Fermi operators. We first prove a novel bound on the temporal fluctuations and then map the equilibration dynamics to a generalized classical XY model in the infinite temperature limit. Using this insight we conjecture that, in most cases, a central limit theorem can be formulated leading to what we call Gaussian equilibration: observables display a Gaussian distribution with relative error where is the dimension of the single particle space. The conjecture, corroborated by numerical evidence, is proven analytically under mild assumptions for the magnetization in the quantum XY model and for a class of observables in a tight-binding model. We also show that the variance is discontinuous at the transition between a quasi-free model and a non-integrable one.
Out of equilibrium and equilibration dynamics of closed quantum systems have been recently at the center of a renewed and intense interest (Kinoshita et al., 2006; Calabrese and Cardy, 2006; Hofferberth et al., 2007; Rigol et al., 2008; Linden et al., 2009; Campos Venuti and Zanardi, 2010a; Campos Venuti et al., 2011). In particular the issue whether quantum integrability plays a key role in equilibration and, if so, which one, has been investigated by several authors (Rigol et al., 2007; Kollath et al., 2007; Barthel and Schollwöck, 2008). In this paper we will address this problem along the lines of the probabilistic approach to quantum equilibration advocated in (Campos Venuti and Zanardi, 2010a, b; Campos Venuti et al., 2011). Here the central object is the (infinite time) full time statistics of the expectation value of a quantum observable. We will focus on a particular, yet very important class of quantum integrable systems: quasi-free fermionic systems, i.e. systems where both the Hamiltonian and the observable are quadratic in the canonical Fermi operators. Using fairly general central-limit type arguments as well as explicit analytic examples we will argue that a sharp distinctive feature of these systems, as opposed to the general interacting ones, is an exponential enhancement of the amplitude of the temporal fluctuations of a quadratic observable around its mean value. This appears to be a precise and quantitative way to make sense of the common folklore that integrability leads to a poorer (or no) equilibration.
The system is initialized in a generic state with particles and both the evolution Hamiltonian and the observable are quadratic in the fermionic operators 111Throughout the paper we use interchangeably the terms quasi-free or quadratic, for observables quadratic in Fermi operators. . The Hamiltonian is (notation , number of sites). The general quadratic observable has the form . We will assume that 222To achieve this, for example, in the translation invariant case, it suffices to have sufficiently fast decaying. as this guarantees that the expectation values of scales at most extensively with the system size 333By diagonalizing and exploiting unitary invariance of the operator norm, one finds . Here the ’s are the eigenvalues of and the ’s the fermionic operators associated to the corresponding eigenvectors. . The main object of investigation is . Exploiting the quadratic nature of the problem and introducing the covariance matrix () one can show that the expectation value reduces to a trace in the one-particle space:
Eq. (1) is perfectly analogous to its many-body version with playing the role of the initial state . There is however one importance difference: while one has , i.e. is extensive (we defined the filling factor).
For unitary evolution in finite systems, the density matrix does not converge neither in the strong nor in the weak topology (Campos Venuti and Zanardi, 2010a). Equilibration must be formulated in probabilistic terms. Given the observation time window , the observable has probability of being in the interval , where the probability density is given by , and the time average operation is . For simplicity we will always take the limit when taking time average and write simply in place of . Roughly speaking an observable equilibrates if its probability density is highly peaked around its mean (Campos Venuti and Zanardi, 2010a). The role of equilibrium state is played by the time-averaged density matrix
The question we are going to address here is: what is the size of the fluctuations of for observables in this quasi-free setting?
Before tackling this question and concentrate on fluctuations, let us point out a few remarks on the nature of the equilibrium state itself. If the spectrum is non-degenerate, the average, dephased, state has the form where are many-body eigenstate of corresponding to energy and . The powers of the Hamiltonian , are linearly independent if there are different eigenvalues and if the spectrum is non-degenerate coincides with the Hilbert space dimension 444Note that one can have a non-degenerate spectrum even in free systems. It suffices, e.g., to have a rational independent one-particle spectrum.. In this latter case the average state can always be written in the form , the so called GGE (generalized Gibbs ensemble) (Rigol et al., 2007). The coefficients depend on the initial state and on the eigenvectors . The condition to write the coefficients in terms of the is precisely that of the invertibility of the Vandermonde matrix (), i.e., once again, non-degeneracy of the spectrum, since . The relation expressing the ’s in terms of the ’s is: . Since , the inverse relation is, in vector notation . The inverse of the Vandermonde matrix can be found explicitly by expanding the identity , multiplying by , summing over and exponentiating see 555Explicitly where , and .
Let us now go back to the context of Eq. (1). Since the time averaged covariance satisfies and so defines a Gaussian state with covariance . Moreover, since , for all , , is a Gaussian state with particles. Now, for what concerns quadratic observables of the kind , their time average expectation value is the same as that obtained with : . In other words the states and are the same when restricted to quadratic observables. A generic Gaussian state can be written in the form where (Peschel, 2003) and is a normalization constant. From this it immediately follows that can be written as where ’s are eigenmodes of and coincides with equation (8) of (Rigol et al., 2007). When the system equilibrates i.e., this remark shows the validity of the GGE for any initial state and all quadratic observables. Indications that converges to in some sense as the size increases were indeed already present in the literature. For instance in (Lanford and Robinson, 1972) it was shown that, in the thermodynamic limit, weakly as , while in (Cramer et al., 2008; Cramer and Eisert, 2010) a particular form of strong convergence was derived, when considering subsystems (though for a Bosonic system).
Ii A bound on the variance
Assuming the non-resonant condition on the energies ( implies and or and ), Reimann has shown (Reimann, 2008) that the temporal fluctuations satisfy where is the maximum minus the minimum eigenvalue of . Now in general, for most initial states , the purity is exponentially small in the system size implying exponentially small fluctuations. Here is one argument. First note that with being the Loschmidt echo: which admits the following cumulant expansion (Campos Venuti and Zanardi, 2010a). Here are the cumulants of computed with . The point is that if is sufficiently clustering, but not an eigenstate of , all the cumulants are extensive in the system size and non-zero. Then where does not depend on the size and, for sufficiently large one has 666When exists, for large one can use the more precise estimate .
In the quasi-free setting the non-resonant condition does not hold. Let us then seek for the analogous of the bound of Reimann in our quasi-free case. Let the one-particle Hamiltonian have the following diagonal form . The time averaged covariance matrix is then . We also define . Assuming the non-resonance condition for the one-particle spectrum, one gets . Now is a non-negative operator and so induces a (possibly degenerate) scalar product which satisfies Cauchy-Schwarz inequality: . This leads us to
Now, since , , we finally obtain While equation (2) is the quasi-free analog of the Reimann’s bound it implies some important differences with respect to the general (non-free) case. Consider the situation where the observable is extensive. In the non-free case the diameter of is extensive, i.e. in spatial dimensions. Moreover, the minimum value of the purity is and so is exponentially small in the system size. In the quasi-free setting, instead, the minimum value of in Eq. (2) is , since defines a density matrix for which the minimum purity is the inverse of the matrix dimension . All in all, recalling that for extensive observables, the minimum of the bound to the variance is , , (for systems of linear size ), whereas in the quasi-free setting one has . This seems to hint at the fact that fluctuations in the quasi-free setting are proportional to the system size and are hence much larger than in the non-free case where they are exponentially small in the volume. Of course Eq. (2) is just an upper bound, and nothing prevents, in principle, from having a much smaller variance. For example, whenever the initial state or the observable commute with the Hamiltonian, is constant and its fluctuations are zero. We will always avoid such pathological situations. In the following, instead, we will argue that the extensive behavior of the fluctuations, in quasi-free systems, is in fact quite general, leading to a scaling of the relative error for a generic observable. Indeed such a scaling has been observed to hold for a quadratic Hamiltonian even for more general observables (see supplementary material of (Cassidy et al., 2011)) and even in presence of disorder except for quenches into a localized phase (Gramsch and Rigol, 2012).
Iii Mapping to a classical XY model
Let us write again the generic expectation value (1) in the basis which diagonalizes : with . To obtain information on the probability density we consider the generating function . Now we observe that if the (one-particle) energies are rationally independent (RI), as a consequence of a theorem on the averages, the infinite time average of is the same as the uniform average over the torus . In this case the generating function is exactly given by the partition function of the generalized, classical XY model with energy and inverse temperature . The matrix defines the lattice of the interactions while the phases give the offset from which the angles are measured. Note that the behavior of the density is dictated by in a neighborhood of which corresponds to infinite temperature of the classical XY model.
It is not difficult to engineer a situation which exactly reproduces the standard XY model in -dimension. For example, it suffices to consider the Hamiltonian , with RI ( is a point of a -dimensional lattice), choose the observable ( indicates nearest neighbor) and initial state . In this case the partition function (free-energy ) of the classical -dimensional XY model is precisely the characteristic function of the observable : (). In fact in this case while so that which defines the nearest neighbor hyper-cubic graph.
We would like to stress here that the one-particle space has a natural underlying geometric structure. For instance, the labels represent points in momentum (real) space in a superfluid (localized) phase and the distance is well defined. Now, when the matrix elements decay sufficiently fast as the corresponding XY model is well defined in the thermodynamic limit, i.e. the intensive free energy has a limit as . This happens for instance in case decays exponentially in or if one has with . When this is the case one has where is the free energy per site. Moreover, under these conditions, one expects to be analytic in the high temperature, , limit, implying that all the cumulants of are extensive. From this we immediately draw the central limit theorem (CLT): as the variable tends in distribution to a Gaussian with zero mean and finite variance given by . We call this situation Gaussian equilibration. It is important to stress that one cannot have Gaussian equilibration in the non-free setting otherwise all the cumulants would scale in the same way. Instead in the non-free setting one expects extensive average but variance exponentially small in the volume.
We will now further corroborate these arguments with two examples where indeed Gaussian equilibration can be proven or shown.
Iv Quench on the quantum XY model
The Hamiltonian is given in terms of Pauli spin operators (we use periodic boundary conditions)
In the quench scenario, the initial state is the ground state of the Hamiltonian with parameters . The parameters are then suddenly changed and is evolved with the Hamiltonian corresponding to .
The model in Eq. (3) has been long used as a testbed for the study of quantum phase transitions in many body systems, and more recently in the realm of out-of equilibrium unitary dynamics. See (Lieb et al., 1961; Barouch et al., 1970) and the more recent monumental (Calabrese et al., 2012) for more details and references. A Jordan-Wigner transformation brings Eq. (3) to a quadratic form in Fermi operators. Since in terms of Fermi operators is , the transverse, total magnetization is a quadratic observable. Its expectation value in the quench setting is given by (Barouch et al., 1970; Campos Venuti and Zanardi, 2010a): where , and are the one-particle energies. The quasi-momenta are quantized according to , . At this point it seems quite natural to expect that the energies are rationally independent. Indeed one can show that the numbers , are rationally independent for prime (Campos Venuti et al., 2011). Given the form of the dispersion we may expect that the requirement that is prime may be removed.
Assuming rational independence of the one-particle energies , the corresponding classical XY model has energy with . Each classical spin ’’ interacts with an external field along a fixed axis with strength . The partition function factorizes, each integral over gives a Bessel function and we obtain . Clearly is a sum of independent random variables, each with zero mean and variance .
Now, under the –quite reasonable– assumption of rational independence of the one-particle energies, one can prove Gaussian equilibration for the observable . More precisely one can show that for any value of parameters , the variable as tends in distribution to a Gaussian with zero mean and variance .
We can prove this result by showing that the Lyapunov condition is satisfied so that the central limit theorem follows from Lindeberg’s theorem (see e.g. (Billingsley, 2012)). Following the notation of (Billingsley, 2012) we have . Then, with . The Lyapunov’s condition with amounts to the vanishing of the following quantity as :
Indeed the RHS of Eq. (4) goes to zero as since for all and for almost any for .
For small quench close to a critical point, the function ( is the quenched variable, ), becomes highly peaked (for instance, the peak is around close to the Ising critical point at and diverges as ). For finite and sufficiently small quench, few terms dominate and one can obtain (for finite ) a non-Gaussian distribution. Indeed, as discussed in detail in (Campos Venuti and Zanardi, 2010b), this is the case in general: for small quenches close to a quantum critical point, observables become a sum of few independent random variables and the distribution acquires a universal double-peaked form.
In figure 1 we show a plot of the distribution of the transverse magnetization per site for the equilibration dynamics undergoing a quantum quench. The distribution agrees very well with a Gaussian with variance obtained considering as a sum of independent variables (Campos Venuti and Zanardi, 2010a). This in turns shows that the assumptions of rational independence seems to be justified or at least that the number of relations among the frequencies is sufficiently small as not to break the CLT. Note that is not a prime in fig. 1. By only checking that the one particle spectrum is non-degenerate, it is easy to prove that the variance is (Campos Venuti and Zanardi, 2010a). Since is a bounded function this implies immediately that , in accordance with the Gaussian equilibration prediction.
V Tight binding model
We consider here the 1- tight binding model with twisted boundary conditions as proposed in (Tasaki, 2010). As quadratic observable we take with diagonal one-particle matrix: . The system is initialized setting all the particles say to the left of the chain, i.e. the initial covariance matrix is with ones and zeros on the diagonal. The observable is extensive for , and the thermodynamic limit is given by , constant and . The time evolved observable reads where the function is given by . The matrix depends only on the difference : . The eigen-energies are given by and the quasimomenta can be considered quantized as , . For the energies are degenerate as but for most of the , does not belong to the Brillouin zone and the energies are non-degenerate. In this case the average is . Still the cannot be RI as , however it is likely that there are few relations among the energies. In fact Tasaki showed (Tasaki, 2010) that for odd and for most of the , the one-particle energies satisfy the non-resonant condition. This implies that, for most , the variance is given by
The double sum above can be evaluated going back to real space with result (assuming otherwise swap with ). This proves analytically that, besides the average, also the second cumulant –the variance– is extensive and in particular one has .
On the other hand, using methods borrowed from statistical mechanics, we can get an approximate analytic expression for the whole cumulant generating function. In this approximation all the cumulants will turn out to be extensive implying Gaussian equilibration for (see figure 2). We first assume rational independence of the energies so that the problem is mapped to an classical XY model. A single (or a few) relation among the energies is a sort of boundary condition for the classical model and is not expected to change the leading, bulk, term. The energy of the classical model for the translation invariant case is ( and is a periodic extension of , i.e. ). Now we note that is highly peaked around , so we approximate the energy keeping only the nearest neighbor term, i.e. . This is precisely a one-dimensional (classical) XY model with periodic boundary condition (and off-set ) and can be solved via transfer matrix method (Mattis, 1984). The partition function becomes where the transfer matrix operator is . The transfer matrix is non-hermitian because of . Using the identity where are Bessel functions (satisfying and ) one sees that plane waves are eigenfunction of with eigenvalue . The largest eigenvalue in modulus, with , gives immediately the free energy in the large size limit: . As expected, in this approximation, the cumulant generating function is extensive and analytic in . One then has again the CLT in its standard form: the variable , tends in distribution to a Gaussian in the thermodynamic limit with variance . The Gaussian prediction is clearly confirmed by a numerical experiment see fig. 2.
Vi Fluctuations and integrability
From the previous discussion it appears that the variance , for a quadratic observable , can be used as an effective tool for the characterization of integrable-non-integrable transition, at least in the case where integrable models are identified with quadratic systems. In particular we expect an increase of when the evolution Hamiltonian crosses a quasi-free point. To check this conjecture let us add a next-nearest neighbor interaction term to Eq. (3) (for simplicity we set ) which is known to break integrability (see e.g. (Beccaria et al., 2006) and references therein). Hence we consider the model
and perform numerically simulated quenches from to . At the model can be mapped to quasi-free fermions while for any non zero the Hamiltonian is not integrable. As shown in fig. 3, the variance of is discontinuous at the integrable point . As expected is larger than . We would like to stress here that the appearance of a discontinuity in an infinity time average, is perfectly legitimate even for finite size systems. The origin of the discontinuity of the variance stems from a massive violation of the non-resonant condition at . Note that the average (not shown) is smooth at , indicating that the degeneracies of the energies are constant around . We also performed numerical simulations keeping the observation time window finite. The effect of a finite is to make the variance smooth at , approaching the discontinuous value with corrections of order .
As noted previously, it is easy to prove that at the integrable point the variance scale as . On the contrary, checking the exponentially small scaling expected at non-integrable points is very difficult numerically as the computation of the variance requires full exact diagonalization which limits the analysis to very short sizes.
Finally, a compelling question is whether similar results generalize to more complex integrable models such as those integrable by Bethe Ansatz. To investigate this scenario we performed preliminary numerical simulations with the following Hamiltonian
with periodic boundary condition. At the Hamiltonian is integrable by Bethe Ansatz. In our simulations we quenched from to , and looked at the statistics of as this is a quadratic observable in the fermionic setting. Our numerics shows that the fluctuations are smooth as crosses the integrable point, which indicates that the non-resonant condition is satisfied (when restricted to the relevant subspace) also at the integrable point. This indeed has to be expected from the form of the many-body energies arising from the Bethe Ansatz 777V. Gritsev, private communication.. However it is still possible that a lack of rational independence of the many-body energies, for Bethe Ansatz integrable models, will show up in a non-analytic behavior of higher cumulants of as a function of the non-integrability parameter .
In figure 4 we plot the distribution function for for various quench experiment. At the integrable point the distribution looks qualitatively different from the distribution obtained at . However further investigations are needed to ascertain whether the dependence of on is analytic at the integrable point . A lack of analyticity would be a signature of a lack of rational independence of the many-body energies.
In this paper we addressed the question of equilibration in quasi-free Fermi systems. The initial state is a general state while the evolution Hamiltonian as well as the observable are quadratic in the Fermi operators. While for general non-free systems the variance is typically exponentially small in the volume, we find that in the quasi-free setting (for extensive observables) both the mean and the variance are proportional to the volume. This hints at the possibility of an underlying central limit theorem, a circumstance that we termed Gaussian equilibration. In this case the properly rescaled observable becomes Gaussian in the large limit, and in general the relative error satisfies . We proved Gaussian equilibration for the magnetization in the quantum XY model assuming rational independence of the one body energies and gave evidence for a class of observables/initial states evolving with a tight-binding model. In all cases Gaussian equilibration was confirmed by numerical simulations. The key insight is a mapping of the equilibration dynamics to a generalized classical XY model at infinite temperature. As a by product we obtain a quantum setting (initial state, Hamiltonian, and observable ), such that the equilibration dynamics of gives the solution of a classical XY model in -dimension, and vice-versa.
A consequence of the above scenario –confirmed by numerical simulations– is that the variance turns out to be discontinuous at the quasi-free point of an otherwise non-integrable Hamiltonian. This shows that the enhancement of the temporal fluctuations of (for a quadratic observable ) may provide a universal and experimentally testable signature of integrability in the context of out of equilibrium dynamics, at least when integrable systems are identified with quadratic models.
Ultimately the origin of the exponential decrease of the signal to noise ratio illustrated in this paper lies in the exponential reduction of the effective phase space entailed by the quasi-free nature of our setup i.e., the many body space gets effectively replaced by the single particle one. A natural question for future research is whether a similar mechanism is at place for a more general class of integrable systems such as Bethe-Anstanz integrable models. Preliminary numerical results indicate that, in this case the variance is smooth and a possible non-analytic behavior at the integrable point must be sought in higher order cumulants.
LCV wishes to thank KITP for the kind hospitality. This research was partially supported by the ARO MURI under grant W911NF-11-1-0268 and by the National Science Foundation under Grant No. NSF PHY11-25915. PZ also acknowledges partial support by NSF grants No. PHY-969969 and No. PHY-803304.
- Kinoshita et al. (2006) T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 900 (2006).
- Calabrese and Cardy (2006) P. Calabrese and J. Cardy, Phys. Rev. Lett. 96, 136801 (2006).
- Hofferberth et al. (2007) S. Hofferberth, I. Lesanovsky, B. Fischer, T. Schumm, and J. Schmiedmayer, Nature 449, 324 (2007).
- Rigol et al. (2008) M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
- Linden et al. (2009) N. Linden, S. Popescu, A. J. Short, and A. Winter, Phys. Rev. E 79, 061103 (2009).
- Campos Venuti and Zanardi (2010a) L. Campos Venuti and P. Zanardi, Phys. Rev. A 81, 022113 (2010a).
- Campos Venuti et al. (2011) L. Campos Venuti, N. T. Jacobson, S. Santra, and P. Zanardi, Phys. Rev. Lett. 107, 010403 (2011).
- Rigol et al. (2007) M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Phys. Rev. Lett. 98, 050405 (2007).
- Kollath et al. (2007) C. Kollath, A. M. Läuchli, and E. Altman, Phys. Rev. Lett. 98, 180601 (2007).
- Barthel and Schollwöck (2008) T. Barthel and U. Schollwöck, Phys. Rev. Lett. 100, 100601 (2008).
- Campos Venuti and Zanardi (2010b) L. Campos Venuti and P. Zanardi, Phys. Rev. A 81, 032113 (2010b).
- Peschel (2003) I. Peschel, Journal of Physics A: Mathematical and General 36, L205 (2003).
- Lanford and Robinson (1972) O. E. Lanford and D. W. Robinson, Communications in Mathematical Physics (1965-1997) 24, 193 (1972).
- Cramer et al. (2008) M. Cramer, C. M. Dawson, J. Eisert, and T. J. Osborne, Phys. Rev. Lett. 100, 030602 (2008).
- Cramer and Eisert (2010) M. Cramer and J. Eisert, New Journal of Physics 12, 055020 (2010).
- Reimann (2008) P. Reimann, Phys. Rev. Lett. 101, 190403 (2008).
- Cassidy et al. (2011) A. C. Cassidy, C. W. Clark, and M. Rigol, Phys. Rev. Lett. 106, 140405 (2011).
- Gramsch and Rigol (2012) C. Gramsch and M. Rigol, arXiv:1206.3570 (2012).
- Lieb et al. (1961) E. Lieb, T. Schultz, and D. Mattis, Annals of Physics 16, 407 (1961).
- Barouch et al. (1970) E. Barouch, B. M. McCoy, and M. Dresden, Phys. Rev. A 2, 1075 (1970).
- Calabrese et al. (2012) P. Calabrese, F. H. L. Essler, and M. Fagotti, arXiv:1204.3911 (2012).
- Billingsley (2012) P. Billingsley, Probability and Measure (John Wiley & Sons, 2012).
- Tasaki (2010) H. Tasaki, arXiv:1003.5424 (2010).
- Mattis (1984) D. C. Mattis, Physics Letters A 104, 357 (1984).
- Beccaria et al. (2006) M. Beccaria, M. Campostrini, and A. Feo, Phys. Rev. B 73, 052402 (2006).