Equations of Motion for the Out-of-Equilibrium Dynamics of Isolated Quantum Systems from the Projection Operator Technique

We present a rigorous framework to obtain evolution equations for the momentum distribution and higher order correlation functions in weakly interacting systems based on the Projection Operator Technique. These equations can be numerically solved in an efficient way. We compare the solution of the equations with known results for 1D models and find an excellent agreement.


Introduction
The out-of-equilibrium dynamics of closed quantum systems is nowadays a very active field of research (see [1] for a review). The central question behind these efforts is to describe the dynamics of an isolated many-body system towards an eventual stationary state after it has been driven out of equilibrium. Of particular interest are the properties of the equilibrium state itself. If simple correlation functions in this equilibrium state can be reproduced by means of a canonical ensemble the system is said to thermalize, which, interestingly, has been shown not to occur in certain situations [2,3]. These issues, that lie at the heart of statistical mechanics, have recently received a very strong impulse from the experiments with cold atomic gases. These systems are almost perfectly isolated from the environment and can remain quantum coherent for a long time compared with typical experiments duration allowing to study in situ the dynamics of large aggregates of interacting quantum particles [2].
From the theoretical perspective, the non-equilibrium dynamics of large, interacting manybody systems settles a formidable challenge. Among the several theoretical approaches at disposal, different kinds of evolution equations for few-body correlations are being widely employed in the contemporary literature (see, for example [4,5,6,7]). In contrast to numerical methods, which are confined either to small systems or short times, as for instance exact diagonalization or time-dependent density matrix renormalization group (t-DMRG), the solution of evolution equations allows to investigate the dynamics of large systems for relatively long times and usually demands a rather modest numerical effort. However, the deduction of these kind of equations starting from first principles, i.e., the many-body Schrödinger equation or its equivalents, often involves a series of uncontrolled approximations [7,8] that, although routinely utilized, cannot be justified in all cases.
On the other hand, the most studied type of evolution equations are kinetic equations, which are markovian, i.e. memoryless, equations for one-particle distribution functions, a typical example being the Boltzmann equation. Even though the markovian limit is well defined in some situations [9,10] it may not be always justified. In particular, in situations in which the system is reluctant to loose memory of the initial conditions (such as integrable models) markovian kinetic equations may be loosing important features of the relaxation dynamics.
In this work we present a method for constructing evolution equations for few body correlation functions in weakly interacting systems that involves no approximations beyond a series expansion in the small parameter quantifying the strength of the interaction. The resulting equations of motion are inherently non-markovian. The approach is based on the Projection Operator Technique (POT) [9,10,11]. We benchmark the method by comparing the solution of the equations with known results from one dimensional (1D) systems.

Model and statement of the problem
We consider a system of interacting spinless fermions with Hamiltonian H = H 0 + αH 1 with where c † (k) and c(k) are fermionic creation and annihilation operators satisfying canonical anticommutation relations, V k 1 ,k 2 k 3 ,k 4 is the momentum-space matrix element of the interaction, (k) is the dispersion relation, n(k) = c † (k)c(k) is the number operator and α is the strength of the interaction. Our results can be easily extended to the bosonic case. The hermiticity of the Hamiltonian and the symmetry in the sum indices impose k 2 ,k 1 , whereV denotes the complex conjugate. Furthermore, we will be interested in the special case of a translationally invariant Hamiltonian in which the particles interact via a pair potential v(x − y). In such case wherev(k) = dr v(r)e ik·r is the Fourier transform of the potential and we have written the antisymmetrized version in order to respect the symmetry conditions of the potential. We are interested in the evolution of the system starting from an arbitrary initial condition given by a density matrix ρ(0) which we leave unspecified until the next section.

Evolution equation for the momentum distribution
Our starting point is the Liouville equation in the interaction representation (h=1): whereÕ(t) = e iH 0 t Oe −iH 0 t is the interaction representation of the operator O and we have introduced the Liouville superoperator Our task is to find approximate solutions to the microscopic dynamics described by the Liouville equation. The POT defines a program for achieving this. We need to first identify the "slow" or "macroscopic" variables in our system and then project the dynamics into the subspace of these slow variables. In a weakly interacting system the occupation number operators emerge as natural slow variables since [H, n(k)] = O(α). To perform the projection we first introduce the "relevant" density matrix where the time-dependent partition function is given by Note thatσ(t) = σ(t). The Lagrange multipliers λ(k, t) enforce the relation: In other words, σ(t) is the density matrix that maximizes the Von Neumann entropy subject to the constraints (5). We note that although the dynamics ofρ(t) is generated by the Liouville equation, the relevant density matrix σ(t) is time-dependent independently of the representation. The projection of the dynamics consists in finding an equation of motion for σ(t). To this end we introduce a projection super-operator P (t) that projects the relevant density matrix P (t)ρ(t) =σ(t). We refer the interested reader to the literature to find an explicit expression for P (t) [9,10,11,12,13]. It is also useful to define the complementary projector Q(t) = 1 − P (t).
Following the usual steps [9,10,11,12], from the Liouville equation we obtain an equation for the slow dynamics is an anti-chronologically ordered exponential. The first term in Eq. (6) is a mean-field term that vanishes due to momentum conservation, the second one is a memory term that can be completely expressed in terms of the past history of the n(k) t 's and the third term is a microscopic noise term that cannot be expressed solely as a function of the slow variables. This last term vanishes if we choose an initial condition that has the same form as the relevant density matrices, i.e., if ρ(0) = σ(0), and we shall in the following circumscribe to this case.
Eq. (6) is equivalent to the Liouville dynamics and, in general, as difficult to solve as the original problem. It sets, however, a good starting point for approximations. To render Eq. (6) tractable we perform a perturbative expansion in the interaction strength using that G(t, s) = I + O(α). Taking the trace n(k) t = Tr[n(k)σ(t)] we finally obtain A great simplification arises since, given the Gaussian structure of σ(t), we can use the Wick pairing rule to evaluate the trace in (7). After a straightforward (but potentially tedious) calculation we obtain the explicit equation of motion where ∆e k,k 2 k 3 ,k 4 = (k) + (k 2 ) − (k 3 ) − (k 4 ) and, in order to ease the notation, we have defined f (k, t) ≡ n(k) t andf (k, t) ≡ 1 − n(k) t . This equation, in slightly different versions, has appeared many times in the literature. In Ref. [10] it was derived using the same tools that we present here but it was used only as an intermediate step to derive the Boltzmann equation whereas in Refs. [6,14] it was heuristically derived and used to study the dynamics of infinite dimensional models and to derive a quantum version of the Boltzmann equation, respectively. We want to stress that (8) is valid also for lattice systems in which only quasi-momentum is conserved.
Although the approach is clearly perturbative, the results go beyond conventional lowest order perturbation theory because the perturbation expansion is performed inside the integrodifferential equations and, therefore, the coupling is involved in a highly non-linear way in the final expressions. To asses the accuracy of the approximation is thus not a straightforward task. One possible alternative is to calculate higher order corrections to Eq. (8), which can be done systematically within the context of the projection operator technique. Then, the error could be calculated as the relative difference between the lowest order and the next-to-leading order solutions [11]. We relegate such error analysis to future work. An alternative approach is to compare results with other techniques. For instance, performing a suitable short-times approximation on (8) (details can be found in [6]) it is possible to rederive results first obtained by Moeckel and Kehrein using the flow-equation approach in lowest order perturbation theory [15]. Moreover, in [6] it was found that the solution of Eq. (8) compares well on the accessible timescales with dynamical mean-field theory results for the infinite dimensional Hubbard model. Below we provide an example for a 1D model in which the solution of the equation (8) captures highly non trivial features of the short-to-intermediate times relaxation. Regarding the validity of the approach for longer timescales than those accessible for current numerical techniques, we should mention that under certain assumptions it can be shown [14] that Eq. (8) reduces, for long enough times, to the quantum Boltzmann equation, which is expected to successfully describe the thermalization of isolated weakly interacting quantum systems after the system has effectively lost the information about the details of the initial conditions [8].
With respect to implementation details, Eq. (8) can be solved using standard techniques for systems of Volterra integral equations [19]. A straightforward algorithm for the solution using, for instance, the trapezoidal rule to perform the time integral, implies a calculation time that scales as L 3D × N 2 , where N is the number of times steps and D the space dimension. We have found an algorithm whose execution time scales as L 3D × N allowing us to reach large sizes and times. We finally note that the evolution equations are very suitable for parallel computing.

Results for a 1D model
Now we present results obtained using (8) in a 1D model of spinless fermions with nearestneighbor hopping and interactions where n j = c † j c j and sums run for j = 0, L − 1 (L is the length of the system) and we impose periodic boundary conditions. This model can be mapped on to the XXZ model by means of a Jordan-Wigner transformation and, in equilibrium, exhibits a Luttinger liquid phase for α < 2J. We study a system at half-filling (k F = π 2 ) and take as initial state the ground state of the Hamiltonian with α = 0, i.e., we consider the "quantum quench" scenario. The same situation was studied in Ref. [16] using t-DMRG where, remarkably, it was found that the short time dynamics were accurately described by the quench dynamics of the Luttinger model (LM) [17,18], raising the question of to what extent the equilibrium low energy fixed point of a given model can describe the out of equilibrium dynamics. In Fig. 1 we show results for two different observables, the jump Z(t) of the momentum distribution at k F and the kinetic energy e kin (t), and compare them with the LM predictions. For Z(t) we find that, after an initial Gaussian decay at very short times, the intermediate time decay is well described by a power law whose exponent is given precisely by the LM predictions. At longer times, beyond an α-dependent time scale t * , deviations occur. These deviations will be studied in detail in a future publication. For e kin (t) we find that the LM prediction works well even beyond the time scale t * associated with Z(t).

Evolution equation for higher order correlations
To climb another step in the hierarchy of correlations we must switch to the Heisenberg representation: Note the difference with the Liouvillian defined in (3). The trace operation defines a dual projection operator over the observables of the Hilbert space: where O is an observable, µ a density matrix and P(t) is the observable space projection operator. In the Heisenberg representation the strategy is to separate the slow and fast components of the evolution operator e iLt using the projector P(t) and its complement Q(t). This allows to find an operator Langevin-like equation for the slow variables [9], where δn(k, t) = e iLt n(k) − n(k) t . The memory function Φ k,k (s, t) can be expressed in terms of the n(k) t 's. In particular, it is possible to show that [9] Φ k,k (t, u) = δ where In the last equation a perturbative expansion analogous to that performed in the interaction representation has been used. The microscopic noise term η k (t) cannot be written entirely in terms of the n(k) t 's. Again, we refer the interested reader to the literature to find an explicit expression for η k (t) [9] since for our immediate purposes we shall not need it. Using the Wick rule and taking the functional derivative Eq. (13) we find an explicit evolution equation for the time-correlation function of the slow variables F k,k (t) ≡ Tr [ρ(0)δn(k, t)δn(k , 0)] where and A k (t, s) = k B k ,k (t, s). To the best of our knowledge, it is the first time that Eq. (15) is derived. In order to solve Eq. (15) we need first to calculate the solution to Eq. (8) to determine the coefficients A k (t, s), B k,q (t, s) and C k,q (t, s). Once this is done, Eq. (15) is as amenable to numerical solution as Eq. (8). Solutions of this equation for specific 1D models shall be presented elsewhere. We remark here, however, that the Heisenberg representation POT could be useful to obtain evolution equations for other correlation functions using the same scheme.

Conclusions
We have presented a rigorous framework to obtain evolution equations for various observables and correlations functions in weakly interacting systems. These equations can be numerically solved in an efficient way. We have put to test the accuracy of the approach comparing with known results for 1D models. Moreover, the POT formalism could allow to obtain higher order corrections to these equations. The inherent limitation of the approach is that it is valid only in the perturbative regime.