Time and spatially resolved quench of the fermionic Hubbard model showing restricted equilibration

We investigate the quench of half-filled one-dimensional and two-dimensional fermionic Hubbard models to models without Coulomb interaction. Since the time propagation is Gaussian, we can use a variety of time-dependent quantum Monte Carlo methods to tackle this problem without generating a dynamical sign problem. Using a continuous time quantum Monte Carlo method, we achieve a system size of 128 sites in one dimension, and using a Blankenbecler-Scalapino-Sugar type algorithm, we were able to simulate $20\times 20$ square lattices. Applying these methods to study the dynamics after the quench, we observe that the final state of the system can be reasonably well described by a thermal single-particle density matrix that takes the initial single-particle conservation laws into account. The characteristic decay toward this limit is found to be oscillatory with an additional power law decay that depends on the dimensionality. This numerically exact result is shown to compare favorably to mean-field approximations as well as to perturbation theory. Furthermore, we observe the information propagation in the one-dimensional case in the charge-charge and spin-spin correlations and find that it is linear with a velocity of roughly $v\approx 4$ in units of the hopping amplitude.

Figure 1

Momentum-resolved $\eta$-pairing correlation function $\langle {\eta }_q^{†}{\eta }_q\rangle \left(t\right)$ for an initial Hubbard interaction of $U=2$ on a $L=128$ site chain. One can clearly see the horizontal black line in the middle of the figure, which is present for all times. This corresponds to the fact that this special value of $q=\pi$ is a conserved quantity.

Figure 2

Time-dependent decay of the magnetization of the mean-field solution with $T=0$ and $U=2$. (a) The behavior of Eq. (16) in dimensions from 1 to 4. The equations for the magnetization were solved self-consistently and afterward the time propagation was calculated. The colors are consistent through all four plots. Now in (b) every data set is multiplied by the expected $t^{\frac{D}{2}}$ behavior. We see that all data sets approach sines with constant amplitude, thereby confirming our conjecture Eq. (17) without introducing the numerical artifacts of the log-log-plot (d). (c) The Fourier transform of the data sets in (b). Finally, (d) is a log-log plot of the absolute value of the data sets in (a). The kinks in the data for 1D and 3D are due to the fact that we have data with a very fine time resolution. The oscillations in these graphs would extend all the way to zero in the plot, thereby effectively coloring the plot in black. Nevertheless, the power-law decay of the envelopes is clearly visible.

Figure 3

Spin-spin correlation function $S(\vec{q}={(\pi ,\pi )}^T,t)$ of a two-dimensional Monte Carlo simulation of a $20\times 20$ lattice at $U=8$ and $\beta =2.5$. We expect this to be in the high-temperature regime where Eq. (33) is valid. The dashed black line is $\propto J_0^4\left(4t\right)$, the result from perturbation theory. The offset is taken from the large time behavior of the QMC data and the amplitude was taken from the value at $t=0$. The inset shows a magnified view of the region below the inset from $t=1$ to $t=4.5$. We see that the approximation works almost flawlessly in that regime. The deviation for $t\to 5$ has its root in boundary effects that set in for approaching that point in time.

Figure 4

Full contour $C$ that enables us to cover imaginary-time evolution and real-time evolution on a common footing. $t^{-}$ is a time on the forward branch, $t^{+}$ is a time on the backward branch, and $\tau$ is a time on the imaginary branch. This contour is parametrized by the contour time $s$, that runs from 0 to $t_{\text{exp}}$ on the forward contour, from $t_{\text{exp}}$ to $2t_{\text{exp}}$ on the backward contour, and from $2t_{\text{exp}}$ to $2t_{\text{exp}}+\beta$ on the imaginary branch. Note that this mapping of imaginary and real time into the complex plane leads imaginary time (the one familiar from thermodynamic perturbation theory) to end up on the imaginary axis; therefore, it's a purely imaginary number.

Figure 5

Double occupancy $\langle n_{\uparrow }{n}_{\downarrow }\rangle$ of the one-dimensional Monte Carlo simulations for different values of the Hubbard interaction $U$ of a 64-site system at $\beta =10$ seems to decay to $0.25$, the value of a free model.

Figure 6

Spin-spin correlations in 1D and 2D. (a) The spin-spin correlation functions $S(q=\pi ,t)$ decay toward the values given by the effective model (dashed line). (b) In the log-log plot we see the power-law-like behavior of $|S(q,t)-S_{\text{eff}}\left(q\right)|$ for $q=\pi$. The decay is roughly $t^{-1}$. Here $N=64$ and $\beta =10$. (c) Also in 2D we see the decay of $S(\vec{q}={(\pi ,\pi )}^T,t)$ toward the effective model. Since these are simulations of a $20\times 20$ system at $\beta =10$, the maximum time is limited to $t\approx 4.5$ as determined from finite-size scaling. (d) In 2D, the envelope of the oscillations of $|S(\vec{q},t)-S_{\text{eff}}\left(\vec{q}\right)|$ for $\vec{q}={(\pi ,\pi )}^T$ follows a long-time decay law $\propto t^{-2}$. For $U=8$ and $\beta =2.5$, this is supported by the dashed line which is $\propto J_0^4\left(4t\right)$

Figure 7

Snapshots of the spatially resolved charge-charge correlation: $|\langle n_0\left(t\right)n_i\left(t\right)\rangle -\langle n_0\left(t\right)\rangle \langle n_i\left(t\right)\rangle |$ for different times. For $t=0$ we see the characteristic exponential decay of an insulator. Between $t=1.6$ and $t=3.2$ we see that a characteristic front forms that is propagating through the lattice. The area behind this front seems to be metallic as evidenced by the lack of an exponential decay. This is a lattice of 128 sites at $\beta =10$ with an initial $U=1$.

Figure 8

For a 64-site lattice, we see that the time in which the correlation functions $C_{\mathrm{corr}}(r,t)=\langle n_0\left(t\right)n_r\left(t\right)\rangle -\langle n_0\left(t\right)\rangle \langle n_r\left(t\right)\rangle$ equilibrate depends on the initial conditions, the chosen $U$. However, for $U=1$ it seems reasonable to think of the short-ranged correlation functions as equilibrated to the effective values $C_{\text{eff}}\left(r\right)$. Note that the $y$ axis has a logarithmic scale.

Figure 9

We see the causality cone in the charge-charge correlation function $|\langle n_0\left(t\right)n_i\left(t\right)\rangle -\langle n_0\left(t\right)\rangle \langle n_i\left(t\right)\rangle |$. Figure 11 corresponds to a vertical cut along the $t$ axis.

Figure 10

Spin-spin correlation functions $\langle S_i^z\left(t\right)S_j^z\left(t\right)\rangle$ also show the horizon, but also more noise. This plot as well as Fig. 9 is for a 128-site lattice at $U=1$ and $\beta =10$.

Figure 11

Exponential suppression of $|\langle n_0\left(t\right)n_i\left(t\right)\rangle -\langle n_0\left(t\right)\rangle \langle n_i\left(t\right)\rangle |$ outside of the causality cone. The blue lines are exponentials meant as a guide to the eye. The indices $n$ correspond to different distances in the measurement of the correlation function.

Figure 12

Comparison of numerically gained mean-field data of a finite chain of 4096 sites with the asymptotic expansion given in Eq. (A8). Here $U=2$, which gives a self-consistently determined $\Delta \approx 0.34$.