Propagation of quantum correlations after a quench in the Mott-insulator regime of the Bose-Hubbard model

Abstract

We study a quantum quench in the Bose-Hubbard model where the tunneling rate J is suddenly switched from zero to a finite value in the Mott regime. In order to solve the many-body quantum dynamics far from equilibrium, we consider the reduced density matrices for a finite number of lattice sites and split them up into on-site density operators, i.e., the mean field, plus two-point and three-point correlations etc. Neglecting three-point and higher correlations, we are able to numerically simulate the time-evolution of the on-site density matrices and the two-point quantum correlations (e.g., their effective light-cone structure) for a comparably large number $\mathcal{O}({10}^3)$ of lattice sites.

1 Introduction

The exponential growth of the underlying Hilbert space is one of the main reasons for our limited ability to simulate quantum many-body systems. In contrast to classical many-particle systems, it is not sufficient to treat the state of each particle or lattice site separately - one has to consider the quantum correlations (entanglement) as well. If these quantum correlations obey certain properties, for example if they are sufficiently short-ranged, there are quite efficient methods for approximating the quantum state of the system, such as matrix-product states [1–4]. However, in typical non-equilibrium situations [5–15] these quantum correlations tend to spread [16–26] and thus are not short-ranged after some time - which is the reason why these methods typically break down in such cases eventually. In the following, we present an alternative method which fully incorporates quantum correlations between pairs of lattice sites at arbitrary distances. The price we have to pay is the neglect of genuine three-point (and higher) correlations. However, comparison with exact numerical diagonalization (for small systems) reveals that this is a reasonable approximation. Even though the method we shall present can be applied to quite general lattice systems, we are going to focus on the Bose-Hubbard model, because the non-equilibrium many-body quantum dynamics in this system can be tested in high-precision experiments using ultra-cold atoms in optical lattices [27–29].

2 Bose-Hubbard model

We consider a system of bosons with local interactions in a lattice described by the Bose-Hubbard Hamiltonian

$$\hat{H}=-\frac{J}{Z}\sum _{{\mu }_1{\mu }_2}{T}_{{\mu }_1{\mu }_2}{\hat{b}}_{{\mu }_1}^{†}{\hat{b}}_{{\mu }_2}+\frac{U}{2}\sum _{\mu }{\hat{n}}_{\mu }({\hat{n}}_{\mu }-1),$$

(1)

where J is the tunneling rate for the nearest neighbors and $U>0$ is the interaction parameter. The creation and annihilation operators ${\hat{b}}_{{\mu }_1}^{†}$ and ${\hat{b}}_{{\mu }_2}$ at the lattice sites ${\mu }_1$ and ${\mu }_2$, respectively, satisfy the standard bosonic commutation relations, and ${\hat{n}}_{\mu }={\hat{b}}_{\mu }^{†}{\hat{b}}_{\mu }$ is the particle-number operator. The lattice structure is encoded in the adjacency matrix $T_{{\mu }_1{\mu }_2}$ which equals unity if the sites labeled by ${\mu }_1$ and ${\mu }_2$ are tunneling neighbors (i.e., if a particle can hop from ${\mu }_1$ to ${\mu }_2$) and zero otherwise. The number of tunneling neighbors at a given site ${\mu }_1$ yields the coordination number $Z={\sum }_{{\mu }_2}{T}_{{\mu }_1{\mu }_2}$. (We assume a translationally invariant lattice with periodic boundary conditions.)

In order to study the time evolution of the system after a sudden change of the parameter J, we consider the density matrix describing the quantum state of the complete lattice which obeys the von Neumann equation ($ħ=1$)

$$i{\partial }_t\hat{\rho }=[\hat{H},\hat{\rho }].$$

(2)

Then we introduce a set of ℓ-point reduced density matrices for the sites ${\mu }_1{\mu }_2\cdots {\mu }_{\ell }$ by tracing over all other sites:

(3)

These reduced density matrices can be split up into their correlated parts ${\hat{\rho }}_{{\mu }_1{\mu }_2\cdots {\mu }_{\ell }}^{\mathrm{corr}}$ plus all possible products of reduced density matrices for single lattice sites and the correlated parts of the reduced density matrices for smaller number of sites - such that the total number of sites for the terms within each product equals to ℓ. For instance, the correlated parts of the two-point and three-point reduced density matrices are in analogy with the cumulant expansion given by

$$\begin{array}{r}{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}={\hat{\rho }}_{{\mu }_1{\mu }_2}-{\hat{\rho }}_{{\mu }_1}{\hat{\rho }}_{{\mu }_2},\\ {\hat{\rho }}_{{\mu }_1{\mu }_2{\mu }_3}^{\mathrm{corr}}={\hat{\rho }}_{{\mu }_1{\mu }_2{\mu }_3}-{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}{\hat{\rho }}_{{\mu }_3}-{\hat{\rho }}_{{\mu }_1{\mu }_3}^{\mathrm{corr}}{\hat{\rho }}_{{\mu }_2}-{\hat{\rho }}_{{\mu }_2{\mu }_3}^{\mathrm{corr}}{\hat{\rho }}_{{\mu }_1}-{\hat{\rho }}_{{\mu }_1}{\hat{\rho }}_{{\mu }_2}{\hat{\rho }}_{{\mu }_3}.\end{array}$$

Inserting this split into the von Neumann equation (2) yields the equations of motion for the reduced density matrices ${\rho }_{{\mu }_1}$ for single lattice sites (i.e., the on-site matrices) and the correlated parts ${\hat{\rho }}_{{\mu }_1{\mu }_2\cdots {\mu }_{\ell }}^{\mathrm{corr}}$, which were already derived in Ref. [30]. The equation for the on-site matrix ${\hat{\rho }}_{{\mu }_1}$ reads

$$i{\partial }_t{\hat{\rho }}_{{\mu }_1}=[{\hat{H}}_{{\mu }_1},{\hat{\rho }}_{{\mu }_1}]+\frac{1}{Z}\sum _{{\mu }_2\ne {\mu }_1}{tr}_{{\mu }_2}\{{\hat{H}}_{{\mu }_1{\mu }_2}+{\hat{H}}_{{\mu }_2{\mu }_1},{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}+{\hat{\rho }}_{{\mu }_1}{\hat{\rho }}_{{\mu }_2}\},$$

(4)

where ${\hat{H}}_{\mu }=(U/2){\hat{n}}_{\mu }({\hat{n}}_{\mu }-1)$ and ${\hat{H}}_{{\mu }_1{\mu }_2}=-J{T}_{{\mu }_1{\mu }_2}{\hat{b}}_{{\mu }_1}^{†}{\hat{b}}_{{\mu }_2}$ are the local and non-local parts of the Hamiltonian.

The analogous equation for the two-point correlators ${\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}$ contains the on-site matrices ${\hat{\rho }}_{{\mu }_1}$, but also the three-point correlations ${\hat{\rho }}_{{\mu }_1{\mu }_2{\mu }_3}^{\mathrm{corr}}$. In the following, we assume that we can neglect these three-point correlations - which is the main assumption of our method. The quality of this approximation will be discussed in the next section. With this assumption, we obtain the approximate equation for ${\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}$:

$$\begin{array}{rl}i{\partial }_t{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}=& [{\hat{H}}_{{\mu }_1},{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}]+\frac{1}{Z}[{\hat{H}}_{{\mu }_1{\mu }_2},{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}+{\hat{\rho }}_{{\mu }_1}{\hat{\rho }}_{{\mu }_2}]\\ -\frac{1}{Z}{\hat{\rho }}_{{\mu }_1}{tr}_{{\mu }_1}\left\{[{\hat{H}}_{{\mu }_1{\mu }_2}+{\hat{H}}_{{\mu }_2{\mu }_1},{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}+{\hat{\rho }}_{{\mu }_1}{\hat{\rho }}_{{\mu }_2}]\right\}\\ +\frac{1}{Z}\sum _{{\mu }_3\ne {\mu }_1,{\mu }_2}{tr}_{{\mu }_3}\left\{[{\hat{H}}_{{\mu }_1{\mu }_3}+{\hat{H}}_{{\mu }_3{\mu }_1},{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}{\hat{\rho }}_{{\mu }_3}+{\hat{\rho }}_{{\mu }_2{\mu }_3}^{\mathrm{corr}}{\hat{\rho }}_{{\mu }_1}]\right\}\\ +({\mu }_1\leftrightarrow {\mu }_2).\end{array}$$

(5)

As a consequence, the two equations (4) and (5) form an approximate closed set of coupled non-linear equations, which we solve numerically.

In the present work, we are interested in the dynamics of the system within the Mott-insulator phase. Due to the fact that the $U(1)$-symmetry is not broken, ${\hat{\rho }}_{\mu }$ is diagonal in the basis of the local Fock states $|n{〉}_{\mu }$ and the equation of motion (4) simplifies in this basis to

$$i{\partial }_t〈{\hat{O}}_{{\mu }_1}^{n_1{n}_1}〉=-\frac{J}{Z}\sum _{{\mu }_2}{T}_{{\mu }_1{\mu }_2}(\sqrt{n_1}{\psi }_{{\mu }_1{\mu }_2}^{n_1,n_1-1}-\sqrt{n_1+1}{\psi }_{{\mu }_1{\mu }_2}^{n_1+1,n_1}-\mathrm{c}.\mathrm{c}.),$$

(6)

where ${\hat{O}}_{\mu }^{n_1{n}_2}=|n_1{〉}_{\mu }〈n_2|$ and we use the notations

$$\begin{array}{r}〈{\hat{O}}_{\mu }^{n_1{n}_2}〉=tr\left\{{\hat{\rho }}_{\mu }{\hat{O}}_{\mu }^{n_1{n}_2}\right\},\\ {〈{\hat{O}}_{{\mu }_1}^{n_1{n}_2}{\hat{O}}_{{\mu }_2}^{n_3{n}_4}〉}^{\mathrm{corr}}=tr\left\{{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}{\hat{O}}_{{\mu }_1}^{n_1{n}_2}{\hat{O}}_{{\mu }_2}^{n_3{n}_4}\right\},\\ {\psi }_{{\mu }_1{\mu }_2}^{n_1{n}_2}={〈{\hat{O}}_{{\mu }_1}^{n_1{n}_2}{\hat{b}}_{{\mu }_2}〉}^{\mathrm{corr}}=\sum _{n_3=0}^{\mathrm{\infty }}\sqrt{n_3+1}{〈{\hat{O}}_{{\mu }_1}^{n_1{n}_2}{\hat{O}}_{{\mu }_2}^{n_3,n_3+1}〉}^{\mathrm{corr}}.\end{array}$$

(7)

Note that in the regime we are interested in (vanishing order parameter: $〈{\hat{b}}_{\mu }〉=0$), $〈{\hat{O}}_{\mu }^{n_1{n}_2}〉$ vanishes for $n_1\ne n_2$. The diagonal terms $〈{\hat{O}}_{\mu }^{nn}〉$ yield the probabilities $p_{\mu }(n)$ to have n particles at the lattice site μ.

As an additional simplification, we neglect the term $[{\hat{H}}_{{\mu }_1{\mu }_2},{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}]/Z$ in Eq. (5) because it is of order $1/Z^2$ (see below). As a result, the correlations ${〈{\hat{O}}_{{\mu }_1}^{n_1,n_2}{\hat{O}}_{{\mu }_2}^{n_3,n_4}〉}^{\mathrm{corr}}$ vanish unless $n_1=n_2\pm 1$ and $n_3=n_4\mp 1$ and thus Eq. (5) simplifies to

$$\begin{array}{r}i{\partial }_t{〈{\hat{O}}_{{\mu }_1}^{n_1+1,n_1}{\hat{O}}_{{\mu }_2}^{n_2,n_2+1}〉}^{\mathrm{corr}}\\ =U(n_2-n_1){〈{\hat{O}}_{{\mu }_1}^{n_1+1,n_1}{\hat{O}}_{{\mu }_2}^{n_2,n_2+1}〉}^{\mathrm{corr}}\\ -\frac{J}{Z}{T}_{{\mu }_1{\mu }_2}\sqrt{n_1+1}\sqrt{n_2+1}(〈{\hat{O}}_{{\mu }_1}^{n_1+1,n_1+1}〉〈{\hat{O}}_{{\mu }_2}^{n_2,n_2}〉-〈{\hat{O}}_{{\mu }_1}^{n_1,n_1}〉〈{\hat{O}}_{{\mu }_2}^{n_2+1,n_2+1}〉)\\ -\frac{J}{Z}\sum _{{\mu }_3\ne {\mu }_1,{\mu }_2}{T}_{{\mu }_3{\mu }_1}\sqrt{n_1+1}(〈{\hat{O}}_{{\mu }_1}^{n_1+1,n_1+1}〉-〈{\hat{O}}_{{\mu }_1}^{n_1,n_1}〉){\left({\psi }_{{\mu }_2{\mu }_3}^{n_2+1,n_2}\right)}^{\ast }\\ -\frac{J}{Z}\sum _{{\mu }_3\ne {\mu }_1,{\mu }_2}{T}_{{\mu }_3{\mu }_2}\sqrt{n_2+1}(〈{\hat{O}}_{{\mu }_2}^{n_2,n_2}〉-〈{\hat{O}}_{{\mu }_2}^{n_2+1,n_2+1}〉){\psi }_{{\mu }_1{\mu }_3}^{n_1+1,n_1}.\end{array}$$

(8)

We have checked numerically that neglecting or including the term $[{\hat{H}}_{{\mu }_1{\mu }_2},{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}]/Z$ does not produce any visible difference in the results presented in this work. For other correlation functions such as $〈{\hat{n}}_{{\mu }_1}{\hat{n}}_{{\mu }_2}〉$, however, the term $[{\hat{H}}_{{\mu }_1{\mu }_2},{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}]/Z$ becomes important.

In the following we will study the dynamics of the system with one boson per site which is initially in the ground state of the Bose-Hubbard Hamiltonian at $J=0$. In this case, $〈{\hat{O}}_{\mu }^{nn}〉={\delta }_{n,1}$ and all correlations ${〈{\hat{O}}_{{\mu }_1}^{n_1+1,n_1}{\hat{O}}_{{\mu }_2}^{n_2,n_2+1}〉}^{\mathrm{corr}}$ vanish. The system of non-linear equations (6) and (8) can be solved using different strategies. One possibility is to make a perturbative expansion with respect to the inverse coordination number $1/Z$ assuming that the reduced density matrices scale as [30]

$${\hat{\rho }}_{\mu }=\mathcal{O}(1),{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}=\mathcal{O}\left(\frac{1}{Z}\right),{\hat{\rho }}_{{\mu }_1{\mu }_2{\mu }_3}^{\mathrm{corr}}=\mathcal{O}\left(\frac{1}{Z^2}\right),\text{etc.}$$

(9)

This method allows to get approximate analytical results, see, e.g., [30–34]. However, exact solutions of the truncated set of non-linear equations (6) and (8) can only be obtained numerically. Before we start to discuss the results obtained using the two strategies, we would like to provide a justification of the neglect of three-point correlations which will be done in the next section.

3 Two-point versus three-point correlations in 1D and 2D

In this section, we present exact numerical results for two-point and three-point correlation functions in a one-dimensional chain consisting of 11 lattice sites and in a two-dimensional square lattices of $3\times 3$ lattice sites. They are obtained by full diagonalization of the Bose-Hubbard Hamiltonian with periodic boundary conditions without any truncation of the Hilbert space. This allows us to calculate exactly the complete time evolution of any quantity as well as their mean values averaged over an infinite time.

In Figure 1, we display the time evolution of the nearest-neighbor two-point correlations $〈{\hat{b}}_1^{†}{\hat{b}}_2〉$ as well as two exemplary three-point correlations $〈{\hat{b}}_1{({\hat{b}}_2^{†})}^2{\hat{b}}_3〉$ and ${〈{\hat{n}}_1{\hat{b}}_2^{†}{\hat{b}}_3〉}^{\mathrm{corr}}\equiv 〈{\hat{n}}_1{\hat{b}}_2^{†}{\hat{b}}_3〉-〈{\hat{n}}_1〉〈{\hat{b}}_2^{†}{\hat{b}}_3〉$. By comparing these three-point correlators to others such as $〈{({\hat{b}}_2^{†})}^2{\hat{b}}_3{\hat{b}}_4〉$, we found that the former $〈{\hat{b}}_1{({\hat{b}}_2^{†})}^2{\hat{b}}_3〉$ is larger than the latter $〈{({\hat{b}}_2^{†})}^2{\hat{b}}_3{\hat{b}}_4〉$ - i.e., the quantities plotted in Figure 1 yield an upper bound. In one dimension, the two-point correlation $〈{\hat{b}}_1^{†}{\hat{b}}_2〉$ oscillates around an average value of about 0.2 while the largest three-point correlator $〈{\hat{b}}_1{({\hat{b}}_2^{†})}^2{\hat{b}}_3〉$ always lies far below $〈{\hat{b}}_1^{†}{\hat{b}}_2〉$ and oscillates around an average value of roughly 0.03. Thus, already in one dimension, the neglect of the three-point correlations seems to be an approximation which is not too bad. As expected from the $1/Z$-arguments above, the two-point correlation $〈{\hat{b}}_1^{†}{\hat{b}}_2〉$ is weaker in two dimensions with an average value of approximately 0.12 while the three-point correlator $〈{\hat{b}}_1{({\hat{b}}_2^{†})}^2{\hat{b}}_3〉$ is even more suppressed and oscillates around an average value of roughly 0.01.

Figure 1

Two-point and three-point correlations. Correlation functions $〈{\hat{b}}_1^{†}{\hat{b}}_2〉$ (red), $〈{\hat{b}}_1{({\hat{b}}_2^{†})}^2{\hat{b}}_3〉$ (blue), ${〈{\hat{n}}_1{\hat{b}}_2^{†}{\hat{b}}_3〉}^{\mathrm{corr}}$ (green) after quench $J/U=0\to 0.1$ calculated by exact diagonalization for one-dimensional lattice of 11 sites (1D) and for two-dimensional lattice of $3\times 3$ sites (2D). Straight horizontal lines show the values averaged over an infinite evolution time.

In summary, although the three-point correlations are not zero, they are indeed smaller than the two-point functions $〈{\hat{b}}_1^{†}{\hat{b}}_2〉$, which justifies our approximation up to a certain accuracy - even in one dimension. Let us roughly estimate the impact of the neglected three-point correlations on a local observable such as $〈{\hat{n}}_{\mu }^2〉$. The time-derivative ${\partial }_t〈{\hat{n}}_{\mu }^2〉$ yields terms containing two-point correlations such as $J〈{\hat{n}}_{\mu }{\hat{b}}_{\mu }^{†}{\hat{b}}_{\nu }〉$ which are fully included in our method. However, the time-derivative of that term gives contributions containing three-point correlations, for example $J^2〈{\hat{b}}_{\lambda }{({\hat{b}}_{\mu }^{†})}^2{\hat{b}}_{\nu }〉$, which are neglected within our method. Thus, in order to achieve an accuracy ε (say, one percent) for the local observable $〈{\hat{n}}_{\mu }^2〉$, we could integrate the evolution equation for a time t until ${(Jt)}^2〈{\hat{b}}_1{({\hat{b}}_2^{†})}^2{\hat{b}}_3〉=\mathcal{O}(\epsilon )$ before the accumulated error induced by the neglect of the three-point correlations may become too large. With $J/U=0.1$, this gives a time of order $10/U$ in two dimensions and a somewhat shorter time in one dimension.

4 First order in $1/Z$

In this section we discuss the approximate analytical solutions of Eqs. (6) and (8) obtained in Ref. [30] in first order of $1/Z$. The probabilities of the occupation numbers larger than two vanish and the probabilities to have zero and two particles are equal to each other. Their time dependence is given by

$$p_{\mu }(0)=p_{\mu }(2)=\frac{4J^2}{N}\sum _{\mathbf{k}}{T}_{\mathbf{k}}^2\frac{1-cos({\omega }_{\mathbf{k}}t)}{{\omega }_{\mathbf{k}}^2},$$

(10)

where ${\omega }_{\mathbf{k}}$ denotes the dispersion relation of the particle-hole excitations

$${\omega }_{\mathbf{k}}=\sqrt{U^2-6JU{T}_{\mathbf{k}}+J^2{T}_{\mathbf{k}}^2},$$

(11)

and $T_{\mathbf{k}}$ is the Fourier transform of the adjacency matrix

$$T_{\mathbf{k}}=\frac{1}{NZ}\sum _{{\mu }_1{\mu }_2}{T}_{{\mu }_1{\mu }_2}exp[i\mathbf{k}\cdot ({\mathbf{x}}_{{\mu }_1}-{\mathbf{x}}_{{\mu }_2})].$$

(12)

For hypercubic lattices in D dimensions (with $Z=2D$), we have

$$T_{\mathbf{k}}=\frac{1}{D}\sum _{d=1}^{D}cos{k}_d,k_d\in \frac{2\pi }{L_d}\mathbb{Z},$$

(13)

where $L_d$ is the size of the lattice along the direction d. The summation over k in Eq. (10) is restricted to the first Brillouin zone.

The time evolution of the one-body density matrix is described by the equation

$$〈{\hat{b}}_{{\mu }_1}^{†}{\hat{b}}_{{\mu }_2}〉=\frac{4JU}{N}\sum _{\mathbf{k}}{T}_{\mathbf{k}}\frac{1-cos({\omega }_{\mathbf{k}}t)}{{\omega }_{\mathbf{k}}^2}exp[i\mathbf{k}\cdot ({\mathbf{x}}_{{\mu }_1}-{\mathbf{x}}_{{\mu }_2})].$$

(14)

For large distances ${\mathbf{x}}_{{\mu }_1}-{\mathbf{x}}_{{\mu }_2}$, we may approximate the k-summation/integration via the saddle-point method. Then, for a given direction in k-space, the maximum group velocity ${\mathbf{v}}_{max}=max{\mathrm{\nabla }}_{\mathbf{k}}{\omega }_{\mathbf{k}}$ determines the maximum propagation speed of correlations, i.e., the effective light cone. In a hypercubic lattice in D dimensions with small J, for example, it is given by $v_{max}\approx 3J/D$ along the lattice axes and by $v_{max}\approx 3J/\sqrt{D}$ along the diagonal (where all the components of ${\mathbf{v}}_{max}$ are equal to each other). A similar result has been obtained in Ref. [20] for the one-dimensional Bose-Hubbard model. For an experimental realization, see, e.g., Ref. [28].

The above approximate analytical solution (10) for the probability of zero occupation number $p_{\mu }(0)$ is compared in Figure 2 with the results obtained by exact diagonalization for one-dimensional and two-dimensional lattices - which was already presented in Ref. [30]. Although Eq. (10) predicts the correct behavior for short times $t=\mathcal{O}(1/U)$, the discrepancy with exact numerical data on a longer time scale is quite noticeable, especially in one dimension. The same feature was also observed for the two-point correlation function $〈{\hat{b}}_{{\mu }_1}^{†}{\hat{b}}_{{\mu }_2}〉$. As we will see in the next section, a significant improvement can be achieved if the calculations are done in a non-perturbative manner.

Figure 2

Approximate analytical solutions. Probability to have no particles on a lattice site in one-dimensional lattice of 11 sites (1D) and for two-dimensional lattice of $3\times 3$ sites (2D). Red lines - exact diagonalization, black lines - Eq. (10). Dashed horizontal lines show the values averaged over an infinite evolution time.

5 Numerical solution of coupled equations

Now we turn to the numerical solution of the coupled set of equations (6) and (8). In order to check the quality of the obtained results, we do calculations first for small systems, where we can compare with exact diagonalization. The data presented in Figure 3 indicate that there is certainly an improvement compared to the perturbative solutions discussed in Section 4. The major difference between the two methods is the following: The perturbative approach discussed in Section 4 is based on a linearization around the Gutzwiller solution ${\hat{\rho }}_{\mu }^0=|1{〉}_{\mu }〈1|$ which facilitates (approximate) analytical solutions - whereas coupled equations (6), (8) fully include the evolution of ${\hat{\rho }}_{\mu }(t)$ and its back-reaction onto the dynamics of the correlations.

Figure 3

Numerical solutions for small lattices. Probability to have no particles on a lattice site in one-dimensional lattice of 11 sites (1D) and for two-dimensional lattice of $3\times 3$ sites (2D). Red lines - exact diagonalization (the same as in Figure 2), black lines - numerical solution of Eqs. (6), (8).

As can be observed in Figure 3, we obtain a quantitative agreement not just for short times $t=\mathcal{O}(1/U)$, but also for intermediate times $t=\mathcal{O}(10/U)$ in one dimension and $t=\mathcal{O}(30/U)$ in two dimensions. These time scales are consistent with the estimate discussed at the end of Section 3 because the data in Figure 3 are directly related to the quantity $〈{\hat{n}}_{\mu }^2〉\approx p_{\mu }(1)+4p_{\mu }(2)\approx 1+2p_{\mu }(0)$ considered there due to $p_{\mu }(0)\approx p_{\mu }(2)\gg p_{\mu }(3)$. Moreover, even for longer times $t=\mathcal{O}(100/U)$, we find that our method reproduces qualitative features reasonably well, especially in two dimensions. There seems to be a shift of the time coordinate of a few percent (for certain modes), which might be explainable by an effective renormalization of J due to the neglected higher-order contributions as discussed in Ref. [30].

Having found that the neglect of three-point correlations yields quantitative agreement for short and intermediate times and reproduces qualitative features for longer time scales, we now apply the same method to larger lattices for long times (which are hardly reachable by other methods). From a minimal standpoint, this is a study of what happens if these three-point correlations are neglected. However, in view of the agreement observed above and since the three-point correlators in large lattices are presumably comparable to those in Figure 1, we expect that this procedure again yields quantitative agreement for short and intermediate times and reproduces qualitative features for longer time scales.

5.1 Dynamics after quench in a large one-dimensional chain

The time dependence of the probabilities $p_{\mu }(n)$ as well as the one-body density matrix $〈{\hat{b}}_{{\mu }_1}^{†}{\hat{b}}_{{\mu }_2}〉$ for a chain of 50 sites is shown in Figure 4. The probabilities $p_{\mu }(0)$ and $p_{\mu }(2)$ almost coincide and higher occupations $p_{\mu }(n)$ with $n\ge 3$ are negligible.

Figure 4

Numerical solutions for large lattice in 1D. Calculations for 50 sites. Upper panel: Probabilities to have zero (red) and two (green) particles on a lattice site. Lower panel: $〈{\hat{b}}_{\mu }^{†}{\hat{b}}_{\mu +s}〉$, $s=1,\dots ,49$. Red line $s=v_{max}t$, $v_{max}=3J$.

After several oscillations at the initial stage, $p_{\mu }(n)$ stabilizes - which is often referred to as prethermalization [35–37]. However, at a later time $t>150/U$, we see a revival of oscillations (followed again by stabilization). Comparison with the dynamics of the correlations displayed in the lower panel of Figure 4 reveals that this revival time roughly coincides with the time the correlations need to move through the whole chain. (Note that this picture is symmetric with respect to $s=25$ due to the imposed periodic boundary conditions.)

This coincidence can be explained via the following intuitive quasi-particle picture: The quantum quench generates correlated particle-hole pairs (analogous to cosmological particle creation, for example [38, 39]) which move away in opposite directions due to quasi-momentum conservation. Initially, the wave-functions of these particles and holes have a large overlap and thus their phase coherence (i.e., their correlation) affects the local probabilities $p_{\mu }(0)$ and $p_{\mu }(2)$ leading to an oscillation (destructive versus constructive interference). When they move away from each other, this overlap decreases and thus the local probabilities $p_{\mu }(0)$ and $p_{\mu }(2)$ settle down to a quasi-stationary value (loss of phase coherence). However, when these correlated particles and holes meet again (after one round trip), their wave-functions overlap once more and thus the oscillation of $p_{\mu }(n)$ induced by (destructive or constructive) interference is repeated. Therefore, these revivals can be regarded as a finite-size effect induced by the periodic boundary conditions: enlarging the system size shifts these revivals to later times.

Consistent with this quasi-particle picture, one can clearly see a propagating front of correlations moving with a constant velocity $v_{\max }$, which corresponds to the time-dependent distance between the correlated particles and holes. This velocity $v_{max}$ agrees quite well with the analytical result obtained in the first order of the $1/Z$ expansion, see the previous section.

5.2 Dynamics after quench in a large two-dimensional square lattice

Numerical results for a two-dimensional lattice of $30\times 30$ sites are presented in Figures 5 and 6. The spatial dependences of the function $〈{\hat{b}}_{{\mu }_1}^{†}{\hat{b}}_{{\mu }_2}〉$ at different moments of time are shown in Figure 5, where each panel extends over the whole lattice and the reference site labeled by ${\mu }_1$ is in the middle of the panels. The four-fold symmetry of these images reflects the lattice geometry. The propagation is anisotropic with the velocities being maximal along the lattice diagonals and minimal along the axes. Similar results were recently obtained by the variational Monte Carlo method [24].

Figure 5

Anisotropic propagation of correlations in 2D. Calculations for $30\times 30$ sites. Coordinate dependences of the correlation function $〈{\hat{b}}_{{\mu }_1}^{†}{\hat{b}}_{{\mu }_2}〉$ at different times after quench: $tU=10$ (a), 20 (b), 30 (c), 40 (d), 50 (e), 60 (f), 70 (g), 80 (h), 90 (i).

Figure 6

Numerical solutions for large lattice in 2D. Calculations for $30\times 30$ sites. Upper panel: Probabilities to have zero (red) and two (green) particles on a lattice site. Middle panel: $〈{\hat{b}}_{\mu }^{†}{\hat{b}}_{\mu +{\mathbf{e}}_{1}s}〉$, $s=1,\dots ,29$. Red line $s=v_{max}t$, $v_{max}=3J$. Lower panel: $〈{\hat{b}}_{\mu }^{†}{\hat{b}}_{\mu +({\mathbf{e}}_1+{\mathbf{e}}_2)s}〉$, $s=1,\dots ,29$. Red line $s=\sqrt{2}{v}_{max}t$, $v_{max}=3J/\sqrt{2}$. ${\mathbf{e}}_1$ and ${\mathbf{e}}_2$ are unit vectors along the lattice axes.

Figure 6 shows the time evolution of the probabilities $p_{\mu }(n)$ as well as of the correlation function $〈{\hat{b}}_{{\mu }_1}^{†}{\hat{b}}_{{\mu }_2}〉$ along a lattice axis (middle panel) and along a diagonal (lower panel). Note that the distances in the latter plot are rescaled by a factor of in comparison to the former one.

Consistent with the $1/Z$ expansion, the overall amplitude of $p_{\mu }(n)$ and $〈{\hat{b}}_{{\mu }_1}^{†}{\hat{b}}_{{\mu }_2}〉$ is reduced by roughly a factor of $1/2$. The propagation velocity of the wave front is also in a good agreement with the analytical predictions of the $1/Z$ expansion.

As in the one-dimensional setup, we see again a revival of oscillations of $p_{\mu }(n)$ after the propagation of the correlations through the whole lattice. However, this revival effect is not so pronounced as in one dimension. We attribute this reduction to the enlarged phase space of the quasi-particles in two dimensions, which results in a weaker phase coherence. Intuitively speaking, the reunions of particle-hole pairs with different momenta do not occur at the same time.

6 Conclusions

We have developed a method which allows us to simulate the dynamics of quantum correlations in lattices of comparably large sizes for relatively long time scales. The formalism is based on the equations of motion (4) and (5) for the reduced density matrices for one and two lattice sites etc. This infinite set of equations is truncated such that the on-site density matrices ${\rho }_{{\mu }_1}$ and two-point correlations ${\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}$ are taken into account exactly whereas three-point and higher correlations are neglected.

This approximation is motivated by a perturbative expansion in powers of the inverse coordination number $1/Z$. As a result, our method is very suited for high-dimensional lattices - which are quite hard to treat within most other approaches (such as matrix product states). However, exact diagonalization for small lattices in one and two dimensions (see Figure 1) shows that the underlying approximation is not too bad in this regime as well: The three-point correlations such as $〈{\hat{b}}_1{({\hat{b}}_2^{†})}^2{\hat{b}}_3〉$ are much smaller than the two-point function $〈{\hat{b}}_1^{†}{\hat{b}}_2〉$, see Figure 1, while other three-point correlations such as ${〈{\hat{n}}_1{\hat{b}}_2^{†}{\hat{b}}_3〉}^{\mathrm{corr}}\equiv 〈{\hat{n}}_1{\hat{b}}_2^{†}{\hat{b}}_3〉-〈{\hat{n}}_1〉〈{\hat{b}}_2^{†}{\hat{b}}_3〉$ are even more suppressed. This observation is partly explained by the fact that we are still comparably deep in the Mott phase. Using strong-coupling perturbation theory type arguments, the two-point function $〈{\hat{b}}_1^{†}{\hat{b}}_2〉$ scales with the first order in $J/U$ while the three-point correlations are second-order effects. Note, however, that we do not use any expansion in $J/U$ - our only approximation is the neglect of the three-point and higher correlations.

Comparison to exact diagonalization for small lattices (in one and two dimensions, see Figure 3) shows that our truncation scheme yields quantitative agreement for short $t=\mathcal{O}(1/U)$ and intermediate $t=\mathcal{O}(10/U)$ times and reproduces qualitative features for longer time scales $t=\mathcal{O}(100/U)$. Since the three-point correlations in large lattices are presumambly comparable to those in small lattices (depicted in Figure 1), we expect that this agreement applies to larger lattices as well.

As an application, we study the dynamics of interacting bosons in a lattice as described by the Bose-Hubbard Hamiltonian (1) after quench within the Mott-insulator regime. We find that the time evolution of the local particle-number distribution $p_{\mu }(n)$ is directly related to the propagation of the correlations $〈{\hat{b}}_1^{†}{\hat{b}}_2〉$. In particular, we find the revival of oscillations of the local quantities $p_{\mu }(n)$ after the correlations traverse the whole lattice. In two dimensions, the propagation of correlations is anisotropic with the velocity being minimal along the lattices axes and maximal along the diagonals.

Of course, we do not claim that our method is generally superior to other approaches such as matrix-product states (MPS) or the density-matrix renormalization group (DMRG). As always, different approaches have their own advantages and drawbacks. Apart from the aforementioned applicability to higher-dimensional lattices and the ability to treat long-range correlations, the advantages of our method are the following: Most importantly, the computational complexity of our method scales polynomially (instead of exponentially) with system size - even in the presence of long-range correlations. Note that this is different for matrix-product states, for example, where the required internal matrix dimension scales exponentially with the entanglement entropy which is proportional to the system size in these non-equilibrium situations [4]. Furthermore, our method can be applied to general lattice structures with periodic or open boundary conditions and does not contain any free fit parameter. In contrast to weak-coupling methods (see, e.g., [10]), it is particularly suited for the regime of strong interactions U.

As an outlook, our method can be extended easily to inhomogeneous lattices with manageable computational overhead - which facilitates taking into account the trap potential as well as disorder potentials. Even though the initial state used in this work had no correlations, it is trivial to incorporate initial correlations (such as a thermal state with finite J, for example). In order to improve the accuracy, it is also possible to shift the truncation by taking into account three-point correlations (see the Note added below), either fully or only within a finite range - or to approximate them by suitable functions of the on-site density matrices ${\rho }_{{\mu }_1}$ and two-point correlations ${\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}$. Finally, our method can be applied to other systems such as spin lattices or the Fermi-Hubbard model, for instance. In this context, it would be interesting to study the relation between our method and time-dependent (i.e., non-equilibrium) dynamical mean-field theory (t-DMFT). This approach is based on the approximation of the self-energy by a local quantity which becomes exact in the limit $Z\to \mathrm{\infty }$. Thus, it displays some similarities to the $1/Z$-expansion in Eq. (9) but with the important difference that the hopping term in Eq. (1) scales with $1/Z$ instead of the $1/\sqrt{Z}$-scaling used in (fermionic) t-DMFT. As a result, the limit $Z\to \mathrm{\infty }$ in t-DMFT becomes more complicated than in our method where all the correlations decay as a power of $1/Z$, see Eq. (9). Nevertheless, it would be interesting to compare our method to t-DMFT, which should be the subject of further studies.

7 Note added

After submitting our manuscript, we completed the first run of simulations which fully include all three-point correlations (in addition to the on-site density matrices and two-point correlations) but neglect four-point and higher correlators. Even though these results are still preliminary, they indicate that the accuracy and reliability for longer time scales can be improved significantly, see the plots for small lattices presented in Figure 7 and compare with those in Figure 3. The same calculations for larger lattices lead to results which are very similar to those presented in Figures 4–6. More details will be published elsewhere.

Figure 7

Numerical solutions for small lattices: three-point correlations. Probability to have no particles on a lattice site in one-dimensional lattice of 11 sites (1D) and for two-dimensional lattice of $3\times 3$ sites (2D). Red lines - exact diagonalization (the same as in Figure 2), black lines - numerical solution of the equations of motion for the reduced density matrices including three-point correlations.

Appendix

Since we approximate the full von Neumann equation (2) within our truncation scheme (by omitting three-point and higher correlations), it is a good idea to ask the question of whether this scheme breaks some of the important conservation laws of the underlying system - which could then lead to unphysical results.

First we check that our truncated set of equations (6) and (8) conserves the trace of the density operator $tr\{\hat{\rho }\}=1$ (conservation of probability). For the reduced density matrices, this implies $tr\{{\hat{\rho }}_{\mu }\}=1$ and $tr\{{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}\}=0$. By taking the trace of equations (4) and (5), we see that these traces are conserved exactly within our truncated scheme. However, one should be aware that the positivity of ${\hat{\rho }}_{\mu }$ is not guaranteed. Even though this was not a problem for our simulations, we encountered numerical instabilities related to this issue in some other extremal cases.

By multiplying Eq. (4) with ${\hat{n}}_{{\mu }_1}$, taking the trace, and summing over ${\mu }_1$, we find that the total number of particles is also conserved

$${\partial }_t〈\hat{N}〉=\sum _{\mu }{\partial }_t〈{\hat{n}}_{\mu }〉=0.$$

(15)

Similarly, one can show that the variance $〈{\hat{N}}^2〉$ of the total particle number is also conserved exactly by equations (4) and (5).

Finally, equations (6) and (8) exactly conserve the total energy which is given by

$$E=〈\hat{H}〉=-\frac{J}{Z}\sum _{{\mu }_1{\mu }_2}{T}_{{\mu }_1{\mu }_2}〈{\hat{b}}_{{\mu }_1}^{†}{\hat{b}}_{{\mu }_2}〉+\frac{U}{2}\sum _{\mu }\sum _{n=0}^{\mathrm{\infty }}n(n-1)p_{\mu }(n),$$

(16)

where

$$〈{\hat{b}}_{{\mu }_1}^{†}{\hat{b}}_{{\mu }_2}〉\equiv tr\left\{{\hat{\rho }}_{{\mu }_1{\mu }_2}^{\mathrm{corr}}{\hat{b}}_{{\mu }_1}^{†}{\hat{b}}_{{\mu }_2}\right\}=\sum _{n=0}^{\mathrm{\infty }}\sqrt{n+1}{\psi }_{{\mu }_1{\mu }_2}^{n+1,n}$$

(17)

are the entries of the one-body density matrix.