Density-matrix quantum Monte Carlo method

We present a quantum Monte Carlo method capable of sampling the full density matrix of a many-particle system at finite temperature. This allows arbitrary reduced density matrix elements and expectation values of complicated nonlocal observables to be evaluated easily. The method resembles full configuration interaction quantum Monte Carlo but works in the space of many-particle operators instead of the space of many-particle wave functions. One simulation provides the density matrix at all temperatures simultaneously, from $T=\infty$ to $T=0$, allowing the temperature dependence of expectation values to be studied. The direct sampling of the density matrix also allows the calculation of some previously inaccessible entanglement measures. We explain the theory underlying the method, describe the algorithm, and introduce an importance-sampling procedure to improve the stochastic efficiency. To demonstrate the potential of our approach, the energy and staggered magnetization of the isotropic antiferromagnetic Heisenberg model on small lattices, the concurrence of one-dimensional spin rings, and the Renyi $S_2$ entanglement entropy of various sublattices of the $6\times 6$ Heisenberg model are calculated. The nature of the sign problem in the method is also investigated.

Figure 1

The inverse-temperature dependence of the fractions of psips on the first four excitation levels during a DMQMC simulation of a 4$\times$4 square antiferromagnetic Heisenberg lattice [40]. A single $\beta$ loop was carried out using an initial population of ${10}^5$ psips on the diagonal elements. The results shown in (a) were obtained without importance sampling. The results shown in (b) were obtained using importance sampling to ensure that every excitation level had a roughly equal number of psips at zero temperature.

Figure 2

The exact and DMQMC energies for the 4$\times$4 antiferromagnetic Heisenberg lattice (Hilbert space dimension $D\approx 1.29\times {10}^4$) with periodic boundary conditions [40]. The $\beta$ step size, $\Delta \beta$, obeys $JN\Delta \beta =0.1$, where $N=16$, and results are shown every 30 iterations. Each simulation consisted of 1000 $\beta$ loops. For the simulation reported in (a), 100 psips were introduced at the start of each $\beta$ loop; for (b), ${10}^5$ psips were introduced. The error bars in (b) are smaller than the size of the markers.

Figure 3

The square of the staggered magnetization for an $8\times 8$ antiferromagnetic Heisenberg lattice ($D\approx 1.83\times {10}^{18}$) using $1.4\times {10}^7$ psips and 143 $\beta$ loops [40]. The importance-sampling procedure described in the main text was applied. The DMQMC value of $\langle {\widehat{M}}^2\rangle$ is plotted every 50th iteration. Error bars, where not visible, are smaller than the size of the marker. The ground-state value of $\langle {\widehat{M}}^2\rangle$ obtained from an SSE simulation [41] is plotted for comparison.

Figure 4

Renyi-2 entropy for square sublattices of a $6\times 6$ antiferromagnetic Heisenberg lattice. Error bars, where not shown, are smaller than the line thickness. The inset zooms in on the $4\times 4$ sublattice, which has the largest error bars. $S_2$ values were taken from every 25th iteration and averaged over $16–180$ loops over the temperature range.

Figure 5

Renyi-2 entropy for strip sublattices of a $6\times 6$ lattice. Error bars are shown in the inset (color scheme as in the main plot) for large $\beta$, where they are largest, and are smaller than the linewidth where not visible. $S_2$ values were taken from every 25th iteration and averaged over $32–240$ loops over the temperature range. As expected, $S_2$ for the $2\times 6$ and $4\times 6$ sublattices tend to the same value in the zero-temperature limit.

Figure 6

Population dynamics and energy estimate for the antiferromagnetic Heisenberg model on a 4$\times$4 triangular lattice with periodic boundary conditions. A fixed shift was used and a single $\beta$ loop performed. (a) shows the emergence of a plateau in the psip population and (b) shows the exact energy and the DMQMC estimate of the energy as functions of $\beta$. Once the plateau has been exited, the DMQMC energy is in good agreement with the exact results. The severity of the sign problem increases with inverse temperature.