Symmetric minimally entangled typical thermal states for canonical and grand-canonical ensembles

TL;DR

This paper introduces symmetric minimally entangled typical thermal states (METTS) for efficient simulation of finite-temperature quantum systems, leveraging symmetries and novel collapse bases to improve convergence and reduce computational costs.

Contribution

It presents new symmetric Fourier and Haar-random collapse bases for METTS, enabling efficient simulation of canonical and grand-canonical ensembles with reduced autocorrelations.

Findings

Symmetric bases improve autocorrelation decay.

Efficient simulation of grand-canonical ensembles.

Enhanced convergence speeds in spin chain models.

Abstract

Based on the density matrix renormalization group (DMRG), strongly correlated quantum many-body systems at finite temperatures can be simulated by sampling over a certain class of pure matrix product states (MPS) called minimally entangled typical thermal states (METTS). When a system features symmetries, these can be utilized to substantially reduce MPS computation costs. It is conceptually straightforward to simulate canonical ensembles using symmetric METTS. In practice, it is important to alternate between different symmetric collapse bases to decrease autocorrelations in the Markov chain of METTS. To this purpose, we introduce symmetric Fourier and Haar-random block bases that are efficiently mixing. We also show how grand-canonical ensembles can be simulated efficiently with symmetric METTS. We demonstrate these approaches for spin-1/2 XXZ chains and discuss how the choice of the collapse bases influences autocorrelations as well as the distribution of measurement values and, hence, convergence speeds.

Equations

$${\hat{\varrho }}_{\beta ,Q}^{\mathrm{c}}:=e^{-\beta {\hat{H}}_Q}\text{on HQ}$$

$${\hat{\varrho }}_{\beta ,\alpha }^{\mathrm{gc}}=e^{-\beta (\hat{H}+\alpha \hat{Q})}.$$

$$|\psi ⟩=\sum _{\mathbf{\sigma }}{A}_1^{{\sigma }_1}{A}_2^{{\sigma }_2}\cdots A_L^{{\sigma }_L}|\mathbf{\sigma }⟩,$$

$$[A_i^{{\sigma }_i}{]}_{a_i,b_i}\ne 0\text{only if}q(a_i)+q({\sigma }_i)=q(b_i).$$

$${\hat{\varrho }}_{\beta }=\sum _{\mathbf{n}}{P}_{\mathbf{n}}|{\varphi }_{\mathbf{n}}⟩⟨{\varphi }_{\mathbf{n}}|,$$

$$|{\varphi }_{\mathbf{n}}⟩:=\frac{1}{\sqrt{P_{\mathbf{n}}}}{e}^{-\beta \hat{K}/2}|\mathbf{n}⟩,P_{\mathbf{n}}:=⟨\mathbf{n}|e^{-\beta \hat{K}}|\mathbf{n}⟩.$$

$${\varphi }_{\mathbf{n}}\to {\varphi }_{{\mathbf{n}}^{\prime }}\to {\varphi }_{{\mathbf{n}}^{\prime \prime }}\to \dots$$

$$P_{\mathbf{n}}{p}_{\mathbf{n}\to {\mathbf{n}}^{\prime }}=|⟨{\mathbf{n}}^{\prime }|e^{-\beta \hat{K}/2}|\mathbf{n}⟩{|}^2=P_{{\mathbf{n}}^{\prime }}{p}_{{\mathbf{n}}^{\prime }\to \mathbf{n}}.$$

$$\{|\mathbf{n}⟩|\hat{Q}|\mathbf{n}⟩=Q|\mathbf{n}⟩\}$$

$$P_{\mathbf{n}}{\sum }_{\tilde{\mathbf{n}}}{p}_{\mathbf{n}\to \tilde{\mathbf{n}}}{p}_{\tilde{\mathbf{n}}\to {\mathbf{n}}^{\prime }}=P_{{\mathbf{n}}^{\prime }}{\sum }_{\tilde{\mathbf{n}}}{p}_{{\mathbf{n}}^{\prime }\to \tilde{\mathbf{n}}}{p}_{\tilde{\mathbf{n}}\to \mathbf{n}},$$

$$|\tilde{k}⟩:=\frac{1}{\sqrt{\mathcal{D}}}\sum _{x=1}^{\mathcal{D}}{e}^{2\pi ikx/\mathcal{D}}|x⟩\text{for k=1,…,D,}$$

$$\text{SWAP}=|\uparrow \uparrow ⟩⟨\uparrow \uparrow |+|\uparrow \downarrow ⟩⟨\downarrow \uparrow |+|\downarrow \uparrow ⟩⟨\uparrow \downarrow |+|\downarrow \downarrow ⟩⟨\downarrow \downarrow |.$$

$$\{|{\uparrow}{\uparrow}\rangle,\frac{1}{\sqrt{2}}\big{(}|{\uparrow}{\downarrow}\rangle+|{\downarrow}{\uparrow}\rangle\big{)},\frac{1}{\sqrt{2}}\big{(}|{\uparrow}{\downarrow}\rangle-|{\downarrow}{\uparrow}\rangle\big{)},|{\downarrow}{\downarrow}\rangle\},$$

$$\hat{H}=\sum _i\left({\hat{S}}_i^x{\hat{S}}_{i+1}^x+{\hat{S}}_i^y{\hat{S}}_{i+1}^y+\Delta {\hat{S}}_i^z{\hat{S}}_{i+1}^z\right),$$

$$|{\mathbf{n}}^{\prime }⟩=\sum _j{\alpha }_j|{\mathbf{n}}_j^{\prime }⟩,$$

$$P_{{{\boldsymbol{n}}}^{\prime}_{j}}=\Big{[}\prod_{s}r^{(j)}_{s}\Big{]}^{2}=\exp\Big{[}2\sum_{s}\log r^{(j)}_{s}\Big{]}.$$

$$|{\varphi }_{{\mathbf{n}}^{\prime }}⟩=\sum _j{c}_j|{\varphi }_{{\mathbf{n}}_j^{\prime }}⟩,\text{where}{c}_j\equiv \frac{{\alpha }_j\sqrt{P_{{\mathbf{n}}_j^{\prime }}}}{\sqrt{{\sum }_k|{\alpha }_k{|}^2{P}_{{\mathbf{n}}_k^{\prime }}}},$$

$$⟨{\varphi }_{{\mathbf{n}}^{\prime }}|\hat{O}|{\varphi }_{{\mathbf{n}}^{\prime }}⟩=\sum _j|c_j{|}^2⟨{\varphi }_{{\mathbf{n}}_j^{\prime }}|\hat{O}|{\varphi }_{{\mathbf{n}}_j^{\prime }}⟩.$$

Full text

I Introduction

The density matrix renormalization group (DMRG) is a powerful numerical technique for the simulation of one-dimensional (1D) strongly correlated quantum systems White1992-69 ; White1993-48 ; Schollwoeck2005 . Its concise formulation in terms of matrix product states (MPS) Fannes1992-144 ; Rommer1997-55 ; Schollwoeck2011-326 provides a framework for the efficient computation of ground states and the real-time evolution of pure states. Three quite different approaches were developed that extend it to the simulation of finite-temperature properties. The historically first was the quantum transfer-matrix renormalization group Nishino1995-64 ; Bursill1996-8 ; Shibata1997-66 ; Wang1997-56 which has some technical complications such as the non-Hermiticity of the quantum transfer matrix. The second one relies on a purification of the mixed state Uhlmann1976 ; Uhlmann1986 ; Nielsen2000 that can be encoded in matrix product form Verstraete2004-93 ; Barthel2016-94 . These matrix product purifications (MPP) were successfully applied to calculate for example finite-temperature correlation and response functions of quantum spin chains Feiguin2005-72 ; Barthel2005 ; Barthel2009-79b . As, e.g., described in Ref. Barthel2013-15 , the employed MPPs are in one-to-one relation with matrix product density operators (MPDO) Zwolak2004-93 . However, MPP computations remain challenging, because the simulation on an enlarged Hilbert space can lead to a considerable growth of entanglement, making them costly at low temperatures. More recently, an alternative approach was developed that avoids the purification and hence the enlarged Hilbert space. Instead, one samples over a cleverly chosen set of pure quantum states, called minimally entangled typical thermal states (METTS) White2009-102 ; Stoudenmire2010-12 . The efficiency of the METTS algorithm is limited by the statistical error induced by the sampling Binder2015-92 .

In this paper, we describe and demonstrate how symmetries can be utilized to improve the efficiency of the METTS algorithm. Recently, it was shown in Ref. Bruognolo2015-92 how grand-canonical METTS simulations of response functions can be made substantially more efficient by switching to symmetric states just before the real-time evolution. However, the actual METTS sampling was unmodified, i.e., symmetries were not employed in the imaginary-time evolution and transitions.

Here, we discuss symmetric METTS algorithms for both canonical and grand-canonical ensembles. If the system and its environment exchange energy and there is a conserved quantity $\hat{Q}$, the equilibrium state of the system is given by the (here, unnormalized) canonical ensemble

[TABLE]

with $\hat{H}_{Q}$ being the projection of the Hamiltonian onto the quantum number $Q$ subspace $\mathcal{H}_{Q}$ of the full Hilbert space $\mathcal{H}=\bigoplus_{Q}\mathcal{H}_{Q}$. Similarly, if system and environment also exchange the quantity associated with $\hat{Q}$, the equilibrium state is the grand-canonical ensemble

[TABLE]

Here, the Lagrange multiplier $\alpha$ fixes the expectation value of $\hat{Q}$. In more complex cases with multiple conserved quantities $\hat{Q}^{(j)}$, one can also consider ensembles like $\exp[{-\beta(\hat{H}_{Q^{(1)}}+\alpha_{2}\hat{Q}^{(2)})}]$.

The transitions in the Markov chain of METTS samples are determined by projective measurements with respect to a collapse basis that can be freely chosen. This choice strongly affects the statistical properties of the resulting samples Stoudenmire2010-12 ; Binder2015-92 . In this work, we introduce novel collapse bases for symmetric METTS simulations. To be able to conserve global quantum numbers $Q$ and to increase efficiency, we go beyond bases of single-site product states and carry out the projective measurements on blocks of several sites. We discuss the influence of the collapse basis choice on the convergence of the algorithm and introduce Fourier and Haar-random block bases which are, as we call it, efficiently mixing.

The structure of this work is as follows. In Sec. II we review how symmetries can be utilized in the matrix product state representation to achieve a significant speedup in the simulations. We briefly summarize the original METTS algorithm without the use of symmetries in Sec. III. In Sec. IV, we discuss how to use symmetries in the simulation of canonical ensembles. We go on to introduce maximally and efficiently mixing (symmetric) collapse bases in Sec. V and summarize the factors influencing convergence speeds in Sec. VI. Section VII applies the techniques to spin-$1/2$ XXZ chains in the canonical ensemble. A symmetric METTS algorithm for the simulation of grand-canonical ensembles is introduced and demonstrated in Sec. VIII. We summarize and conclude in Sec. IX.

II Symmetries in matrix product states

Let us consider a lattice system with $L$ sites and orthonormal on-site basis states $\{|\sigma_{i}\rangle\,|\,\sigma_{i}=1,\dotsc,d\}$. Matrix product states (MPS) for the system have the form

[TABLE]

with ${{\boldsymbol{\sigma}}}:=(\sigma_{1},\dotsc,\sigma_{L})$ and $D_{i-1}\times D_{i}$ matrices $A^{\sigma_{i}}_{i}$. The $D_{i}$ are also called bond dimensions and we require $D_{0}=D_{L}=1$ such that the matrix product yields a scalar coefficient.

We are concerned with states $|\psi\rangle$ that are eigenstates of a conserved quantity $\hat{Q}$. For simplicity, we restrict our considerations to a single Abelian symmetry such as conservation of total particle number or magnetization. However, everything generalizes in a very similar manner to the cases of multiple conserved quantities and non-Abelian symmetries McCulloch2002-57 ; Weichselbaum2012-327 . For the latter, one exploits that dependencies inside each multiplet are given by Clebsch-Gordan coefficients as exemplified by the Wigner-Eckart theorem.

For an Abelian symmetry, the conserved quantity may have the form $\hat{Q}=\sum_{i}\hat{q}_{i}$ with $\hat{q}_{i}|\sigma_{i}\rangle=q(\sigma_{i})|\sigma_{i}\rangle$ 111The total quantum number can arise less trivially than by summing single-site quantum numbers, e.g., by addition modulo $N$ for a cyclic group of order $N$.. We can construct an MPS (3) with quantum number $Q$ by imposing selection rules on the tensor elements $[A^{\sigma_{i}}_{i}]_{a_{i},b_{i}}$. Specifically, one assigns quantum numbers $q(a_{i})$ and $q(b_{i})$ to the matrix indices and imposes selection rules like

[TABLE]

Setting $q(a_{1})=0$ and $q(b_{L})=Q$, this ensures that $\hat{Q}|\psi\rangle=Q|\psi\rangle$.

Explicitly enforcing these conditions by decomposing the tensors into symmetry blocks leads to a significant speedup and improved accuracy of the MPS algorithms. The numerically most costly operations are typically singular value decompositions (SVD) of the tensors in the MPS. These then reduce to cheaper SVDs of the symmetry blocks.

III Minimally entangled typical thermal states

The strategy employed in the minimally entangled typical thermal states (METTS) algorithm White2009-102 is to decompose the thermal density matrix ${\hat{\varrho}}_{\beta}:=\exp(-\beta\hat{K})$, with $\hat{K}=\hat{H}_{Q}$ or $\hat{K}=\hat{H}-\alpha\hat{Q}$, into a sum of projectors

[TABLE]

with METTS

[TABLE]

Here, $\mathcal{B}:=\{|{{\boldsymbol{n}}}\rangle\}$ is an appropriate orthonormal basis for the full system with $\langle{{\boldsymbol{n}}}|{{\boldsymbol{n}}}^{\prime}\rangle=\delta_{{{\boldsymbol{n}}},{{\boldsymbol{n}}}^{\prime}}$. As discussed later, we choose block product states. Block product states may be entangled within blocks of a certain number of lattice sites but not across the block boundaries. Correspondingly, the MPS bond dimensions $D_{i}$ are 1 at the block boundaries.

As described in Refs. White2009-102 ; Stoudenmire2010-12 ; Binder2015-92 , one can efficiently generate a Markov chain

[TABLE]

of METTS in MPS form according to their (unnormalized) probabilities $P_{{\boldsymbol{n}}}$ by repeated imaginary-time evolution steps $|{{\boldsymbol{n}}}\rangle\to|\phi_{{\boldsymbol{n}}}\rangle$ and projective measurements $|\phi_{{\boldsymbol{n}}}\rangle\to|{{\boldsymbol{n}}}^{\prime}\rangle$. The evolution step can be executed with time-dependent DMRG (tDRMG) White2004 ; Daley2004 . The transition probabilities $p_{{{\boldsymbol{n}}}\to{{\boldsymbol{n}}}^{\prime}}=|\langle{{\boldsymbol{n}}}^{\prime}|\phi_{{\boldsymbol{n}}}\rangle|^{2}$ obey detailed balance

[TABLE]

Thermal expectation values $\langle\hat{O}\rangle_{\beta}=\operatorname{Tr}({\hat{\varrho}}_{\beta}\hat{O})/\operatorname{Tr}({\hat{\varrho}}_{\beta})$ can then be computed by averaging $\langle\phi_{{{\boldsymbol{n}}}^{\nu}}|\hat{O}|\phi_{{{\boldsymbol{n}}}^{\nu}}\rangle$ over the Markov chain. If the states $|{{\boldsymbol{n}}}\rangle$ are (block) product states, the projective measurements can be executed in a sweep through the lattice by doing local projective measurements White2009-102 ; Stoudenmire2010-12 ; Binder2015-92 .

In the following, we discuss how symmetries can be utilized in METTS simulations for canonical and grand-canonical ensembles, i.e., how the conservation of $\hat{Q}=\sum_{i}\hat{q}_{i}$ eigenvalues can be used to substantially reduce computation costs. Please note that Ref. Bruognolo2015-92 shows how to produce symmetric METTS for the grand-canonical ensemble for the evaluation of time-dependent quantities (study of quenches or response functions). To this purpose, non-symmetric METTS, which are not $\hat{Q}$ eigenstates, have been generated. Symmetric METTS are then obtained from these in subsequent symmetric collapses. While this does not provide any computational advantage for the imaginary-time evolution and the evaluation of static quantities, it can make subsequent real-time evolutions of the METTS, in which the symmetries are exploited, much more efficient Bruognolo2015-92 . What is described in the following offers an efficient way to already utilize symmetries during the imaginary-time evolution.

IV Symmetries for the canonical ensemble

Conceptually, it is straightforward to simulate canonical ensembles (1) using symmetric METTS. One simply needs to restrict the initial state and the collapse basis $\{|{{\boldsymbol{n}}}\rangle\}$ to the correct symmetry sector. In particular, one should work with an orthonormal basis

[TABLE]

of $\mathcal{H}_{Q}$. If these states are (block) product states, they can be easily encoded as symmetric MPS (3) with small bond dimensions, where matrix elements obey the constraint (4). As the Hamiltonian $\hat{H}_{Q}$ commutes with $\hat{Q}$, the symmetry constraints on the MPS also hold during the imaginary-time evolution $|{{\boldsymbol{n}}}\rangle\to|\phi_{{{\boldsymbol{n}}}}\rangle$. Now, the projective measurements $|\phi_{{\boldsymbol{n}}}\rangle\to|{{\boldsymbol{n}}}^{\prime}\rangle$ need to be done such that also $|{{\boldsymbol{n}}}^{\prime}\rangle$ has quantum number $Q$. If we use a symmetric collapse basis (i.e., every basis state is a $\hat{Q}$ eigenstate), we always stay in the same symmetry sector as transition probabilities $p_{{{\boldsymbol{n}}}\to{{\boldsymbol{n}}}^{\prime}}=|\langle{{\boldsymbol{n}}}^{\prime}|\phi_{{\boldsymbol{n}}}\rangle|^{2}$, in the projective measurements, vanish for states $|{{\boldsymbol{n}}}^{\prime}\rangle$ with quantum number $Q^{\prime}\neq Q$. Inside the symmetry sector with quantum number $Q$, detailed balance is fulfilled as in Eq. (5).

In practice, using the same symmetric collapse basis (6) for every transition in the Markov chain can however be very inefficient, because it often leads to strong autocorrelations between subsequent METTS samples. This is obvious for infinite temperature, where we would be stuck in the initial state of the Markov chain. To arrive at an efficient algorithm, one needs to alternate between different symmetric bases. To this purpose, we introduce novel symmetric collapse bases and discuss their properties in the following (Sec. V).

After completion of this work, we noticed Ref. Lacki2015-91 . To our knowledge it is the only previous work trying to simulate canonical ensembles with METTS. In particular, a canonical ensemble for the Bose-Hubbard model in the gapped Mott regime was simulated, using only the $\{\hat{n}_{i}\}$ eigenbasis, i.e., Fock states. Because of the strong autocorrelations, only every 200th METTS sample was included in the final ensemble.

V Efficient collapse bases

It is possible and often advantageous to switch between different collapse bases in order to decrease autocorrelation times in the Markov chain. For example, one can do projective measurements using a basis $\{|{{\boldsymbol{n}}}\rangle\}$ for all odd iteration steps and a second basis $\{|\tilde{{{\boldsymbol{n}}}}\rangle\}$ for all even iteration steps. Detailed balance is still fulfilled in every second iteration step, as

[TABLE]

where $|{{\boldsymbol{n}}}\rangle$ and $|{{\boldsymbol{n}}}^{\prime}\rangle$ are from basis 1 and $|\tilde{{{\boldsymbol{n}}}}\rangle$ from basis 2.

A simple example is to collapse alternatingly to $\{\hat{S}^{z}_{i}\}$ and $\{\hat{S}^{x}_{i}\}$ eigenstates, respectively, for a spin-$1/2$ system as described in Ref. Stoudenmire2010-12 . For a general system with a $d$-dimensional local state space (e.g., the Bose-Hubbard model with a maximum of $n_{\max}=d-1$ particles per site), one can generate Haar-random collapse bases 222With Haar random, we refer to bases being drawn from a probability distribution that is uniform with respect to the Haar measure, i.e., invariant under unitary transformations. Of course, one can also use other probability distributions. for each iteration step and lattice site Binder2015-92 . Note that in both cases, these measurements break the symmetry associated with the conservation of total magnetization $\hat{S}^{z}_{\mathrm{tot}}=\sum_{i}\hat{S}^{z}_{i}$ or total particle number $\hat{N}_{\mathrm{tot}}=\sum_{i}\hat{n}_{i}$, respectively, because the basis states are not symmetry eigenstates. Such METTS computations hence simulate the grand-canonical ensemble and symmetries can in general not be utilized.

V.1 Maximally mixing bases

If we use a single collapse basis $\{|{{\boldsymbol{n}}}\rangle\}$ such as $\{\hat{S}^{z}_{i}\}$ eigenstates for a spin system or $\{\hat{n}_{i}\}$ eigenstates for a system of bosons or fermions, there is no dynamics at all at infinite temperature ($\beta=0$). Starting from an arbitrary initial state $|{{\boldsymbol{n}}}\rangle$, transitions to all other basis states are impossible such that the METTS simulation is stuck in the state $|{{\boldsymbol{n}}}\rangle$.

Having the METTS dynamics at high temperatures in mind, we can minimize autocorrelation times by switching between collapse bases $\{|{{\boldsymbol{n}}}\rangle\}$ and $\{|\tilde{{{\boldsymbol{n}}}}\rangle\}$ for which the distribution of overlaps $|\langle\tilde{{{\boldsymbol{n}}}}|{{\boldsymbol{n}}}\rangle|$ is as flat as possible. In other words, we want that all overlaps $|\langle\tilde{{{\boldsymbol{n}}}}|{{\boldsymbol{n}}}\rangle|$ are as small as possible. This guarantees that, at least at high temperatures, transitions to many states are possible and of similar probability. For $\beta=0$ and a Hilbert space of dimension $\mathcal{D}$ 333$\mathcal{D}=d^{N}$ for the Hilbert space of an $N$-site system with single-site Hilbert space dimension $d$., an optimal combination of bases yields transition probabilities $p_{{{\boldsymbol{n}}}\to\tilde{{{\boldsymbol{n}}}}}=|\langle\tilde{{{\boldsymbol{n}}}}|{{\boldsymbol{n}}}\rangle|^{2}=1/\mathcal{D}$ $\forall_{{{\boldsymbol{n}}},\tilde{{{\boldsymbol{n}}}}}$ 444This is the case where the expansion coefficients of any state $|{{\boldsymbol{n}}}\rangle$ with respect to the second basis $\{|\tilde{{{\boldsymbol{n}}}}\rangle\}$ all have absolute value $1/\sqrt{\mathcal{D}}$. and also $p^{(2)}_{{{\boldsymbol{n}}}\to{{\boldsymbol{n}}}^{\prime}}:=\sum_{\tilde{{{\boldsymbol{n}}}}}p_{{{\boldsymbol{n}}}\to\tilde{{{\boldsymbol{n}}}}}p_{\tilde{{{\boldsymbol{n}}}}\to{{\boldsymbol{n}}}^{\prime}}=1/\mathcal{D}$ $\forall_{{{\boldsymbol{n}}},{{\boldsymbol{n}}}^{\prime}}$. More generally, we call a sequence of $K$ bases maximally mixing, if $p^{(K)}_{{{\boldsymbol{n}}}\to{{\boldsymbol{n}}}^{\prime}}=1/\mathcal{D}$ $\forall_{{{\boldsymbol{n}}},{{\boldsymbol{n}}}^{\prime}}$ at infinite temperature, where $p^{(K)}$ are the transition probabilities after $K$ steps, i.e., when having cycled once through all $K$ collapse bases.

One example of maximally mixing bases are the bases $\{|{{\boldsymbol{\sigma}}}\rangle\}$ and $\{|\tilde{{{\boldsymbol{\sigma}}}}\rangle\}$ of $\{\hat{S}^{z}_{i}\}$ and $\{\hat{S}^{x}_{i}\}$ eigenstates for a spin-$1/2$ system Stoudenmire2010-12 . Here, $|\sigma_{i}\rangle\in\{|{\uparrow}\rangle,|{\downarrow}\rangle\}$ and $|\tilde{\sigma}_{i}\rangle=\frac{1}{\sqrt{2}}(|{\uparrow}\rangle\pm|{\downarrow}\rangle)$. Unfortunately, this choice of bases is not applicable if we want to exploit the conservation of the total magnetization.

We may call a sequence of $K$ bases efficiently mixing if, at infinite temperature, there are many nonzero $K$-step transition probabilities $p^{(K)}_{{{\boldsymbol{n}}}\to{{\boldsymbol{n}}}^{\prime}}$ of comparable (small) amplitude for every ${{\boldsymbol{n}}}$, and if the transitions grant ergodicity.

In the following, we give specific examples for maximally mixing or at least efficiently mixing collapse bases that are also applicable for symmetric METTS, i.e., when global quantum numbers are conserved.

V.2 Symmetric collapse bases with efficient mixing

Inspired by the discrete Fourier transform, one can construct maximally mixing bases for any $\mathcal{D}$-dimensional Hilbert space. With a first orthonormal basis $\{|x\rangle\,|\,x=1,\dotsc,\mathcal{D}\}$, we can choose a second basis as

[TABLE]

which we call the Fourier basis with respect to $\{|x\rangle\}$. As $|\langle\tilde{k}|x\rangle|^{2}=1/\mathcal{D}$ $\forall{x,k}$ we have indeed $p^{(1)}_{x\to\tilde{k}}=1/\mathcal{D}$ and also $p^{(2)}_{x\to x^{\prime}}=1/\mathcal{D}$ $\forall{x,x^{\prime}}$ at infinite temperature. Note that we are free to reorder the states $|x\rangle$. Hence, permutations can be used to construct different versions of the Fourier basis.

In principle, we could use this approach to construct global symmetric collapse bases that are maximally mixing. Let us discuss this using the example of a spin system. Given an orthonormal basis $\{|{{\boldsymbol{\sigma}}}\rangle\}$ of $\{\hat{S}^{z}_{i}\}$ eigenstates for the magnetization $M$ subspace $\mathcal{H}_{M}$, we can identify the states $|{{\boldsymbol{\sigma}}}\rangle$ with $|x\rangle$, where $x=1,\dotsc,\dim\mathcal{H}_{M}$, and obtain their Fourier basis $\{|\tilde{k}\rangle\}$ according to Eq. (8). These two bases of $\mathcal{H}_{M}$ are maximally mixing. We call the METTS scheme in which one alternates between collapses in this Fourier basis and the $\{\hat{S}^{z}_{i}\}$ eigenbasis “SF-Sz” (Symmetric Fourier w.r.t. $\hat{S}^{z}_{i}$). It is illustrated in Fig. 1.

However, such a global Fourier transform is technically infeasible because the basis states $|\tilde{k}\rangle$ are in general highly entangled. On average, their entanglement entropy is extensive, causing exponentially growing computation costs. To avoid this problem while retaining the efficient mixing property, we divide the lattice into blocks of $b$ lattice sites. For each of these blocks, we have an $\{\hat{S}^{z}_{i}\}$ eigenbasis and can construct a Fourier basis for each symmetry sector of that block. This is scheme “SF$b$-Sz” in Fig. 1. For a spin-$1/2$ system with magnetization conservation, the $2^{b}$-dimensional Hilbert space of a $b$-site block is decomposed into $b+1$ symmetry sectors with magnetizations $-b/2,-b/2+1,\dots,+b/2$ and corresponding subspace dimensions $\tbinom{b}{0},\tbinom{b}{1},\dots,\tbinom{b}{b}$. For each of these subspaces, we construct collapse bases according to Eq. (8) with $|x\rangle$ referring in this case to the $\{\hat{S}^{z}_{i}\}$ eigenstates of some fixed magnetization on the block (in some ordering). These bases are maximally mixing within each symmetry sector, but there are no transitions between different symmetry sectors of a block. Finally, in order to achieve ergodicity also at infinite temperature and to enhance the dynamics of fluctuations in general, we can shift blocks in every second block collapse by $b/2$ sites as illustrated in Fig. 1.

It is instructive to shortly discuss the nature of the resulting METTS dynamics for an example. The simplest case is a spin-$1/2$ system in the symmetry sector with a single up-spin, $M=-L/2+1$. First, consider infinite temperature, $\beta=0$. With global maximally mixing bases such as the symmetric Fourier $\hat{S}_{i}^{z}$ bases (SF-Sz), every second collapse, the up-spin jumps with equal probability to any new site. In the corresponding block scheme SF$b$-Sz, the up-spin does a random walk (diffusive) whose average step size is proportional to the block size $b$. In contrast, for zero temperature ($\beta\to\infty$), every METTS is equal to the ground state in $\mathcal{H}_{M}$, for which, for Hamiltonians of interest, the up-spin is delocalized. So, at low temperatures, the change in position of the up-spin after each Sz collapse is not caused by the mixing property of the bases, but mainly by the delocalization due to the imaginary-time evolution.

Several variations of the above construction of efficiently mixing symmetric bases are conceivable. One is to use Haar-random collapse bases instead of the Fourier bases. So, instead of applying Eq. (8) to obtain the second basis, one can draw a Haar-random basis for each symmetry sector of each $b$-site block. While such collapse bases, for the same block size $b$, perform very similarly at high temperatures, their efficiencies at finite temperatures may be quite different and will in general also depend on the system parameters. For the spin systems, we call the corresponding symmetric METTS schemes “SR$b$-Sz”.

For both the SF$b$-Sz and SR$b$-Sz bases, where blocks are shifted in every second collapse, we noticed that in practice, the intermediate Sz collapses for every second sample are not really necessary. They actually result in somewhat slower convergence in our exemplary benchmark simulations. This can again be understood by considering the overlaps of the different basis states. We denote the schemes where Sz collapses are omitted by SF$b$ and SR$b$.

The symmetric Fourier and random block bases naturally have non-symmetric counterparts that can be applied in simulations of grand-canonical ensembles without symmetries. One simply omits the partitioning of the Hilbert space into symmetry sectors and uses Eq. (8) or the Haar-random choice in the full $b$-site Hilbert space. Such non-symmetric Fourier block bases are always maximally mixing.

Let us briefly mention another variation that we pursued for the construction of efficiently mixing symmetric collapse bases. It is based on the swap operator that acts on a pair of lattice sites and swaps their quantum states. For two spins-$1/2$, it is

[TABLE]

Let $\text{SWAP}_{i,j}$ act on sites $i$ and $j$ of the lattice. As $[\text{SWAP}_{i,j},\hat{S}^{z}_{\mathrm{tot}}]=0$, we can use a symmetric eigenbasis of the swap operator,

[TABLE]

as a collapse basis on pairs of sites. Ideally, we would randomly select arbitrary pairs of lattice sites $(i,j)$ and choose the swap eigenbasis (9) on each of these pairs, forming the global collapse basis as their tensor product (“SWAP” in Fig. 1). To avoid extensive entanglement, we again restrict the random pairing of sites to blocks of $b$ lattice sites (“SWAP$b$” in Fig. 1). As before, we also shift the blocks by $b/2$ lattice sites in every second collapse. For block size $b=2$, the symmetric Fourier basis coincides with the swap basis.

VI Factors influencing convergence speeds

Clearly, the model and the system parameters strongly affect the convergence of observables in the METTS sampling algorithm. For a given system, we can significantly influence the convergence properties by the choice of the collapse bases. This influence is mediated by two key mechanisms. First, the choice of the collapse bases $\{|{{\boldsymbol{n}}}^{(\kappa)}\rangle\}$ ($\kappa=1,\dotsc,K$) determines the METTS ensembles $\{|\phi^{(\kappa)}_{{\boldsymbol{n}}}\rangle\}$ from which we sample. If we could draw independent samples without any autocorrelations, the statistical error of an observable $\langle\hat{O}\rangle$ would solely depend on the distribution of its measurement values $\{\langle\phi^{(\kappa)}_{{\boldsymbol{n}}}|\hat{O}|\phi^{(\kappa)}_{{\boldsymbol{n}}}\rangle\}$ and their corresponding probabilities $P^{(\kappa)}_{{\boldsymbol{n}}}$ in the METTS ensembles. Second, the collapse bases influence the strength of autocorrelations between METTS samples.

To see the interplay of collapse bases and observables, consider, e.g., the operator $\hat{S}^{+}_{0}\hat{S}^{-}_{3}$ in a spin-$1/2$ system with spin-flip symmetry. At high temperatures ($\beta\to 0$), every METTS sample equals its corresponding basis state up to corrections of order $\beta$, i.e., $|\phi_{{\boldsymbol{n}}}\rangle=|{{\boldsymbol{n}}}\rangle+\mathcal{O}(\beta)$. Now, if we collapse into states $|{{\boldsymbol{n}}}\rangle$ that are eigenstates of $\hat{S}^{z}_{0}$ and $\hat{S}^{z}_{3}$, we have $\langle{{\boldsymbol{n}}}|\hat{S}^{+}_{0}\hat{S}^{-}_{3}|{{\boldsymbol{n}}}\rangle=0\ \forall{{\boldsymbol{n}}}$ and, hence, $\langle\phi_{{\boldsymbol{n}}}|\hat{S}^{+}_{0}\hat{S}^{-}_{3}|\phi_{{\boldsymbol{n}}}\rangle=\mathcal{O}(\beta)$. The distribution of measurement values in this METTS ensemble is peaked around the expectation value $\langle\hat{S}^{+}_{0}\hat{S}^{-}_{3}\rangle=\mathcal{O}(\beta)$, leading to small statistical errors. If, however, we choose basis states $|{\boldsymbol{m}}\rangle$ that are eigenstates of $\hat{S}^{x}_{0}$ and $\hat{S}^{x}_{3}$, we have $\langle{\boldsymbol{m}}|\hat{S}^{+}_{0}\hat{S}^{-}_{3}|{\boldsymbol{m}}\rangle=\pm\frac{1}{4}\ \forall{\boldsymbol{m}}$, and $\langle\phi_{{\boldsymbol{m}}}|\hat{S}^{+}_{0}\hat{S}^{-}_{3}|\phi_{{\boldsymbol{m}}}\rangle=\pm\frac{1}{4}+\mathcal{O}(\beta)$ with probabilities $\frac{1}{2}$, leading to large statistical errors. Note that if our observable of interest happened to be, say, $\hat{S}^{z}_{0}$, the effect would be exactly reversed. Hence, there is no generally optimal choice of collapse bases as the answer also depends on the observables of interest. At low temperatures, the METTS become similar to the ground state of the system and, provided it is non-degenerate, the distribution of measurement values is strongly peaked around the ground state expectation value for any collapse basis.

For the effect of autocorrelations, consider first the naïve way of simulating canonical ensembles by using the same symmetry eigenbasis $\{|{{\boldsymbol{n}}}\rangle\}$ for every collapse (e.g., “Sz” in Fig. 1). At infinite temperature, as discussed previously, the Markov chain cannot leave its initial state because $p_{{{\boldsymbol{n}}}\to{{\boldsymbol{n}}}^{\prime}}=|\langle{{\boldsymbol{n}}}^{\prime}|\phi_{{\boldsymbol{n}}}\rangle|^{2}=\delta_{{{\boldsymbol{n}}},{{\boldsymbol{n}}}^{\prime}}$. For small values of $\beta$, the probability that two subsequent samples are equal is still high because $|\langle{{\boldsymbol{n}}}^{\prime}|\phi_{{\boldsymbol{n}}}\rangle|^{2}=\delta_{{{\boldsymbol{n}}},{{\boldsymbol{n}}}^{\prime}}+\mathcal{O}(\beta)$. Thus, the collapse may occasionally induce a few transitions, but autocorrelations between subsequent samples remain high. One obtains slow diffusive dynamics in the Markov chain (left panel in Fig. 2). This can be resolved by using a sequence of efficiently mixing bases. For a spin-$1/2$ system with magnetization conservation, we can, e.g., alternate between Sz collapses and the corresponding symmetric Fourier or Haar-random block bases. The larger we choose the blocks, the faster the Markov chain explores the state space due to reduced autocorrelations (center and right panel in Fig. 2).

The influence of the system parameters becomes more important at lower temperatures. As the temperature is decreased, the METTS $|\phi_{{\boldsymbol{n}}}\rangle$ get closer and closer to the ground state. Typically, deviations from the ground state are localized if the system is gapped, and delocalized if the system is critical. For our block collapse bases, localized deviations then lead again to diffusive dynamics in the METTS and autocorrelations may considerably depend on specifics of the chosen bases. Delocalized deviations generally lead to short autocorrelation times.

VII Results for XXZ chains in the canonical ensemble

In the following, we demonstrate the influence of the collapse bases on convergence speeds for spin-$1/2$ XXZ chains Bethe1931 ; Cloizeaux1966 ; Mikeska2004 with Hamiltonian

[TABLE]

with varying values of the anisotropy parameter $\Delta$. We focus on symmetric METTS computations for the canonical ensemble (CE) with zero magnetization and also compare the convergence speeds to those of non-symmetric simulations of the grand-canonical ensemble (GCE) with $\langle\hat{S}^{z}_{\mathrm{tot}}\rangle=0$. The convergence behavior for different temperatures, observables, and collapse bases is shown in Fig. 3. It demonstrates how our novel collapse bases can significantly improve the convergence in practice.

The statistical error of the METTS sampling is quantified as follows. First, we apply the purification approach Verstraete2004-93 , which has been adapted to also describe canonical ensembles and utilize symmetries Nocera2016-93 ; Barthel2016-94 , to compute reference values of observables. For these simulations, we use the same system parameters, ensemble, and system size as for the corresponding METTS computation, and set a very low DMRG truncation threshold such that we can regard the obtained reference values as quasi-exact. For the error of $N$ METTS (ErrorN in Figs. 3, 4, and 7), we generate several sets of $N$ subsequent samples, and take the root mean square of the average (absolute) deviations of the observable from the quasi-exact reference data in each set 555In practice, when reference values are not available, one can apply standard techniques for the estimation of statistical errors in Monte Carlo simulations, e.g., jackknife resampling or bootstrapping, to probe the efficiency of different collapse bases.. For the METTS simulations, the DMRG truncation weight $\epsilon$ and the Trotter time step $\Delta\tau$ were chosen such that the corresponding errors are negligible compared to the statistical error. See, e.g., Ref. Binder2015-92 for details or Ref. Schollwoeck2005 for a general review.

In Fig. 3, the first two columns show results for the isotropic antiferromagnet ($\Delta=1$ in Eq. (10)) at inverse temperatures $\beta=2$ and $8$, and the third column shows results for $\Delta=3$ in the gapped Néel phase at inverse temperature $\beta=8$. The rows correspond to different observables: the energy per site $\langle\hat{H}/L\rangle$ and the correlators $\langle\hat{S}^{+}_{0}\hat{S}^{-}_{1}\rangle$ and $\langle\hat{S}^{+}_{0}\hat{S}^{-}_{3}\rangle$. Site $i=0$ is located at the center of the lattice, i.e., correlators are evaluated near the center to minimize finite-size effects 666Specifically, site 0 is the left of the two central sites for even system sizes $L$. We tested the influence of the exact position of site 0 on the correlators. It is negligible except for very high temperatures. Also, one could exploit spatial self-averaging Binder2015-92 , but we did not in this study.. At high temperatures ($\beta=2$), using only the Sz collapse leads to strong autocorrelations that slow down the convergence for all three observables shown here. This problem can be alleviated by choosing the efficiently mixing symmetric Fourier or symmetric Haar-random collapse bases. In many cases, the symmetric Fourier bases work slightly better than the symmetric random bases. Errors reduce with increasing block sizes $b$. Larger block sizes typically also increase the computation costs per sample. For the correlator $\langle\hat{S}^{+}_{0}\hat{S}^{-}_{3}\rangle$, the simple Sz collapse leads, for the chosen temperatures, to relatively small statistical errors. This can be attributed to the fact that the distribution of measurement values is peaked around the (small) expectation values $\langle\hat{S}^{+}_{0}\hat{S}^{-}_{3}\rangle$, while it is broader for our novel block bases SF$b$ and SR$b$. Still, the latter reduce autocorrelations between samples: The curves corresponding to our new block bases roughly follow the $1/\sqrt{N}$ convergence, while the Sz-curves start rather flat, which is due to the autocorrelations. In the combination of both effects, the novel block bases typically outperform the simple Sz collapse, sometimes reducing errors by an order of magnitude.

For the critical system at low temperatures (center column in Fig. 3), the METTS errors are relatively independent of the chosen collapse bases. Our interpretation is that, in this case, the METTS are similar to the ground state with deviations from it being delocalized such that the particular choice of the collapse bases in the blocks is not as decisive. For the gapped phase at low temperatures (right column in Fig. 3), however, the METTS errors vary significantly for the different bases. Our interpretation is that, in this case, deviations from the ground state are much more localized such that the distributions of measurement values and autocorrelations depend considerably on the basis choice. For both values of $\Delta$, the Fourier bases (SF$b$) provide the best results.

To compare convergence speeds, Fig. 3 also shows METTS errors for a simulation of the GCE using the non-symmetric METTS algorithm, where one alternates between the $\{\hat{S}_{i}^{z}\}$ and $\{\hat{S}_{i}^{x}\}$ collapse bases. The comparability is of course somewhat limited as the CE and GCE are not equivalent for our finite systems. Overall, the statistical errors in the simulation of the CE with our novel symmetric collapse bases are comparable to and sometimes considerably smaller than those of the GCE simulation using the Sx-Sz bases. Hence, beyond cases where one specifically wants or needs to study the CE, simulating the CE with symmetries can be an efficient variant of the METTS algorithm when one is interested in the thermodynamic limit. This is possible because, as we learn from statistical mechanics, the different equilibrium ensembles such as the CE (1) with quantum number $Q$ and the GCE (2), with $\alpha$ tuned such that $\langle\hat{Q}\rangle^{\mathrm{gc}}_{\alpha}=Q$ are equivalent in the thermodynamic limit.

The symmetry-breaking versions of the Fourier and Haar-random bases can also be applied in simulations of the GCE without symmetries. In many cases, this improves the convergence of the sampling as illustrated for the antiferromagnetic spin-$1/2$ Heisenberg chain in Fig. 4. Here, we use Haar-random symmetry-breaking collapse bases (R$b$) on blocks of $b=1,2,4,$ and $8$ sites and compare the statistical errors of the thermal energy per site. The behavior observed here is typical: The statistical errors reduce as the block size is increased. Note, however, that larger block sizes also tend to increase the computation cost per sample.

VIII Symmetric simulation of grand-canonical ensembles

Except for Ref. Lacki2015-91 , METTS have so far been used exclusively to simulate grand-canonical ensembles (2) White2009-102 ; Stoudenmire2010-12 ; Alvarez_2013 ; Bonnes_2014 ; Binder2015-92 ; Bruognolo2015-92 , alternating between different collapse bases that break the symmetry of the system to reduce autocorrelations. In this case, basis states are not eigenstates of the conserved quantities ($\hat{Q}$) and it is hence in general not possible to exploit symmetries in the MPS representation (3) of the METTS.

In the following, we discuss a modification of the algorithm to simulate grand-canonical ensembles under utilization of symmetries.

VIII.1 The symmetric METTS algorithm

As illustrated in Fig. 5, we employ small collections of symmetric METTS to exploit symmetries in the DMRG time evolution. Starting from an initial symmetric basis state $|{{\boldsymbol{n}}}\rangle$ with some quantum number $Q$, we use imaginary-time evolution to obtain the corresponding symmetric METTS sample $|\phi_{{\boldsymbol{n}}}\rangle$, where symmetries can be utilized in the tDMRG algorithm as described in Sec. II. After the evaluation of observables, we apply a collapse that breaks the symmetry in such a way that the resulting basis state $|{{\boldsymbol{n}}}^{\prime}\rangle$ only has contributions from a small number of different symmetry sectors. We write the new basis state as a superposition of (normalized) symmetry eigenstates,

[TABLE]

where $|{{\boldsymbol{n}}}^{\prime}_{j}\rangle$ is the component with quantum number $Q_{j}$. Subsequently, the $|{{\boldsymbol{n}}}^{\prime}_{j}\rangle$ are evolved in imaginary time separately, again exploiting symmetries.

The DMRG time evolution White2004 ; Daley2004 usually entails a Trotter decomposition. After each time step $s\Delta\tau\to(s+1)\Delta\tau$, the MPS should be renormalized to avoid numerical problems due to the exponential norm decay. So, after each time step, we multiply the evolved state $j$ by a factor $1/r^{(j)}_{s}$. In order to keep the different renormalizations in the separate time evolutions of the states $|{{\boldsymbol{n}}}^{\prime}_{j}\rangle$ consistent, we have to determine their relative weights $P_{{{\boldsymbol{n}}}^{\prime}_{j}}=\langle{{\boldsymbol{n}}}^{\prime}_{j}|e^{-\beta\hat{H}}|{{\boldsymbol{n}}}^{\prime}_{j}\rangle$. These can be obtained by multiplying the renormalization factors $r^{(j)}_{s}$,

[TABLE]

As indicated, we do not actually multiply the $r^{(j)}_{s}$ but rather accumulate the sum of their logarithms because the $P_{{{\boldsymbol{n}}}^{\prime}_{j}}$ decay exponentially in the system size and inverse temperature $\beta$. One should take care that the factors $r^{(j)}_{s}$ do not comprise the norm change due to the DMRG truncations of bond dimensions but only the norm change due to the application of evolution operators $e^{-\Delta\tau\hat{H}}$.

With the resulting normalized states $|\phi_{{{\boldsymbol{n}}}^{\prime}_{j}}\rangle$, the normalized METTS sample is then given by

[TABLE]

and $|\phi_{{{\boldsymbol{n}}}^{\prime}_{j}}\rangle=P_{{{\boldsymbol{n}}}^{\prime}_{j}}^{-1/2}\,e^{-\beta\hat{H}/2}|{{\boldsymbol{n}}}^{\prime}_{j}\rangle$. However, note that we can evaluate any observable $\hat{O}$ without explicitly encoding the superposition $|\phi_{{{\boldsymbol{n}}}^{\prime}}\rangle$ as an MPS according to $\langle\phi_{{{\boldsymbol{n}}}^{\prime}}|\hat{O}|\phi_{{{\boldsymbol{n}}}^{\prime}}\rangle=\sum_{kj}c_{k}^{*}c_{j}\langle\phi_{{{\boldsymbol{n}}}^{\prime}_{k}}|\hat{O}|\phi_{{{\boldsymbol{n}}}^{\prime}_{j}}\rangle$. For symmetry-conserving observables $[\hat{O},\hat{Q}]=0$ we have $\langle\phi_{{{\boldsymbol{n}}}^{\prime}_{k}}|\hat{O}|\phi_{{{\boldsymbol{n}}}^{\prime}_{j}}\rangle=0$ for $k\neq j$ as $Q_{k}\neq Q_{j}$ $\forall\ k\neq j$, and the expression reduces to the simple sum

[TABLE]

Subsequently, the small collection of symmetric METTS is collectively collapsed to a new single symmetric basis state $|{{\boldsymbol{n}}}^{\prime\prime}\rangle$, using, e.g., any of the symmetric collapse bases described in Sec. V. The transition probabilities for this projective measurement can be obtained using $\hat{O}=|{{\boldsymbol{n}}}^{\prime\prime}\rangle\langle{{\boldsymbol{n}}}^{\prime\prime}|$ in Eq. (13). We end up with a single state $|{{\boldsymbol{n}}}^{\prime\prime}\rangle$ that has one of the quantum numbers $\{Q_{j}\}$. This procedure is repeated, alternating between symmetry-breaking and symmetric collapse bases to explore the different symmetry sectors according to their weights in the GCE until the estimates of observables have reached the desired accuracy.

VIII.2 Suitable symmetry-breaking bases

As described, we suggest to alternate between symmetric and symmetry-breaking collapse bases, where elements of the latter should only have components from a few different symmetry sectors to allow for an efficient simulation.

For a spin-$1/2$ system with conservation of the total magnetization, a simple choice is based on the Sz collapse. For even iteration steps, we can use the $\{\hat{S}^{z}_{i}\}$ eigenbasis. For odd iteration steps, we randomly select $n_{x}$ sites on which the $\hat{S}^{x}_{i}$ eigenbasis is used as the collapse basis, while the $\hat{S}^{z}_{i}$ eigenbasis is used on all other sites (“Sz/Sx”). A new set of $n_{x}$ sites is drawn for every symmetry-breaking collapse. It is straightforward to obtain the symmetric components, $\{|{{\boldsymbol{n}}}^{\prime}_{j}\rangle\}$ in Eq. (11), of the non-symmetric basis states $|{{\boldsymbol{n}}}^{\prime}\rangle$ in MPS form. Note that with the described choice of bases, the parameter $n_{x}$ is directly related to the number of states in the small collection of symmetric METTS described above. Specifically, the states $|{{\boldsymbol{n}}}^{\prime}\rangle$ have $n_{x}+1$ components with different magnetizations that are separately evolved in imaginary time. Correspondingly, the total magnetizations of subsequent symmetric samples ($|\phi_{{\boldsymbol{n}}}\rangle$ and $|\phi_{{{\boldsymbol{n}}}^{\prime\prime}}\rangle$ in Fig. 5) differ by $\Delta M\in\{-n_{x},\dots,n_{x}\}$.

In generalization of this symmetric Sz-Sz/Sx scheme, one can reduce autocorrelations, e.g., by starting from any of the efficiently mixing symmetric block bases discussed in Sec. V and, for odd iteration steps, modify the collapse basis for some randomly selected blocks to allow for changes of the conserved quantities, similar to the role of the $\hat{S}^{x}_{i}$ eigenbasis above.

VIII.3 Quantum number trajectories and convergence

We test the algorithm for the isotropic spin-$1/2$ Heisenberg antiferromagnet, corresponding to $\Delta=1$ in Eq. (10). For collapses, we alternate between 4-site block symmetric Fourier bases (SF4 with blocks shifted by 2 sites in every second use), and nonsymmetric bases, where the $\hat{S}^{x}_{i}$ eigenbasis is used for $n_{x}$ randomly selected sites and the $\hat{S}^{z}_{i}$ eigenbasis for all other sites. We denote this combination of bases by “SF4-Sz/Sx”.

Figure 6 shows trajectories of the total magnetization $M(\nu)=\langle\phi_{{{\boldsymbol{n}}}^{\nu}}|\hat{S}^{z}_{{\mathrm{tot}}}|\phi_{{{\boldsymbol{n}}}^{\nu}}\rangle$ occurring in the Markov chain and the probability distribution $\{p_{M}\}$ of the magnetization in the GCE, obtained by the MPDO approach that was introduced in Ref. Barthel2016-94 . The rows in Fig. 6 show data for different values of $n_{x}$ in the symmetry-breaking collapse. The columns correspond to different inverse temperatures $\beta=1$, $4$, and $8$. Whenever we use the symmetric collapse basis (even steps $\nu$ in the Markov chain), the METTS sample is an $\hat{S}^{z}_{{\mathrm{tot}}}$ eigenstate with magnetization $M(\nu)$. For odd $\nu$, $M(\nu)$ is the $\hat{S}^{z}_{{\mathrm{tot}}}$ expectation value. Note that for $n_{x}=L$, we alternate between the $\{\hat{S}_{i}^{x}\}$ eigenbasis and the SF4 bases (SF4-Sx collapse).

In all cases, the Markov chain appropriately explores the symmetry sectors according to their weights. At lower temperatures, the distribution is peaked around the center ($M=0$), while at higher temperatures, a large number of symmetry sectors contribute significant weight to the GCE. In the zero-temperature limit $\beta\to\infty$, where the system is in its ground state with quantum number $M=0$, the symmetry-breaking collapse forces the resulting superposition to have components with $M\neq 0$. However, in this case, the imaginary-time evolution essentially projects onto the ground state, such that all contributions with $M\neq 0$ are exponentially suppressed. The parameter $n_{x}$ determines the speed with which the Markov chain can explore the symmetry sectors. It limits the maximum change in magnetization that can be achieved in a single symmetry-breaking collapse.

Comparing with quasi-exact purification data Barthel2016-94 , Fig. 7 shows the error of the algorithm as a function of the number of samples. The DMRG truncation weight $\epsilon$ and the Trotter step-size $\Delta\tau$ were again chosen such that the error is dominated by the statistical error. See Ref. Binder2015-92 for a detailed discussion of the interplay of different error sources in the METTS algorithm. We consider the same observables as in Fig. 3, namely the energy per site $\langle\hat{H}/L\rangle$ and the correlators $\langle\hat{S}^{+}_{0}\hat{S}^{-}_{1}\rangle$ and $\langle\hat{S}^{+}_{0}\hat{S}^{-}_{3}\rangle$, all at inverse temperatures $\beta=1$, $4$, and $8$. The error of the energy expectation value decreases with increasing $n_{x}$ until $n_{x}=L/2$. For the nearest-neighbor correlator $\langle\hat{S}^{+}_{0}\hat{S}^{-}_{1}\rangle$ the variations with $n_{x}$ are rather small. For the correlator $\langle\hat{S}^{+}_{0}\hat{S}^{-}_{3}\rangle$ at the highest temperature ($\beta=1$), errors increase with increasing $n_{x}$, and, at the lower temperatures ($\beta=4,8$), the variations with respect to $n_{x}$ are again rather small. These properties can again be explained through the competition between autocorrelations and the spread in the distribution of measurement values.

In general, Fig. 7 shows that also for rather small values of $n_{x}$, the error of our new symmetric METTS algorithm is comparable to the error of the original simulation of the GCE without the use of symmetries, where one alternates between $\{\hat{S}^{z}_{i}\}$ and $\{\hat{S}^{x}_{i}\}$ eigenbases on all sites (Sx-Sz collapse). Hence, the use of symmetries can make METTS simulations significantly more efficient.

IX Conclusions and discussion

We have demonstrated how the METTS algorithm can be modified to allow for the utilization of symmetries. Conceptually, this is straightforward for the canonical ensemble – one just needs to employ collapse bases that respect the symmetries of the system. In practice, it is important to cycle through a set of different collapse bases to reduce autocorrelation times. To this purpose, we have introduced and tested efficiently mixing collapse bases such as Fourier bases and Haar-random bases which involve states that are entangled within blocks of several sites. We have also introduced an efficient algorithm for the simulation of the grand-canonical ensemble under utilization of symmetries, using small collections of symmetric METTS.

Explicitly encoding symmetries in the MPS representation of the quantum states leads to a considerable speedup in the imaginary-time evolution and can hence make the sampling significantly more efficient. For both ensembles, we have discussed the effect of the collapse bases on the convergence of the METTS algorithm. Good bases grant short autocorrelation times in the Markov chain of METTS samples and a narrow distribution of measurement values for the observables of interest.

As demonstrated in this paper, understanding the role of the collapse bases and finding improved bases is a promising route to enhancing the efficiency of METTS simulations.