{'fraction_non_alphanumeric': 0.04161858411519743, 'fraction_numerical': 0.007601224223571321, 'mean_word_length': 4.439819868995633, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Symmetric minimally entangled typical thermal states for canonical and grand-canonical ensembles

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.


I. INTRODUCTION
The density matrix renormalization group (DMRG) is a powerful numerical technique for the simulation of one-dimensional (1D) strongly correlated quantum systems [1][2][3]. Its concise formulation in terms of matrix product states (MPS) [4][5][6] 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 finitetemperature properties. The historically first was the quantum transfer-matrix renormalization group [7][8][9][10] 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 [11][12][13] that can be encoded in matrix product form [14,15]. These matrix product purifications (MPP) were successfully applied to calculate for example finite-temperature correlation and response functions of quantum spin chains [16][17][18]. As, e.g., described in Ref. [19], the employed MPPs are in one-to-one relation with matrix product density operators (MPDO) [20]. 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) [21,22]. The efficiency of the METTS algorithm is limited by the statistical error induced by the sampling [23].
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. [24] 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 quantityQ, the equilibrium state of the system is given by the (here, unnormalized) canonical ensembleˆ Here, the Lagrange multiplier α fixes the expectation value ofQ. In more complex cases with multiple conserved quantitiesQ (j) , one can also consider ensembles like exp[−β(Ĥ Q (1) + α 2Q (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 [22,23]. 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 {|σ i | σ i = 1, . . . , d}. Matrix product states (MPS) for the system have the form with σ := (σ 1 , . . . , σ L ) and D i−1 × D i matrices A σ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 |ψ that are eigenstates of a conserved quantityQ. 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 [25,26]. 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 formQ = iq i withq i |σ i = q(σ i )|σ i [27]. We can construct an MPS (3) with quantum number Q by imposing selection rules on the tensor elements [A σi i ] ai,bi . Specifically, one assigns quantum numbers q(a i ) and q(b i ) to the matrix indices and imposes selection rules like Setting q(a 1 ) = 0 and q(b L ) = Q, this ensures that Q|ψ = Q|ψ .
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.
Here, B := {|n } is an appropriate orthonormal basis for the full system with n|n = δ n,n . 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.
Thermal expectation values Ô β = Tr(ˆ βÔ )/ Tr(ˆ β ) can then be computed by averaging φ n ν |Ô|φ n ν over the Markov chain. If the states |n are (block) product states, the projective measurements can be executed in a sweep through the lattice by doing local projective measurements [21][22][23].
In the following, we discuss how symmetries can be utilized in METTS simulations for canonical and grandcanonical ensembles, i.e., how the conservation ofQ = iq i eigenvalues can be used to substantially reduce computation costs. Please note that Ref. [24] 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 notQ 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 realtime evolutions of the METTS, in which the symmetries are exploited, much more efficient [24]. 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 {|n } to the correct symmetry sector. In particular, one should work with an orthonormal basis of 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Ĥ Q commutes withQ, the symmetry constraints on the MPS also hold during the imaginary-time evolution |n → |φ n . Now, the projective measurements |φ n → |n need to be done such that also |n has quantum number Q. If we use a symmetric collapse basis (i.e., every basis state is aQ eigenstate), we always stay in the same symmetry sector as transition probabilities p n→n = | n |φ n | 2 , in the projective measurements, vanish for states |n with quantum number Q = 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. [30]. 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 {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 {|n } for all odd iteration steps and a second basis {|ñ } for all even iteration steps. Detailed balance is still fulfilled in every second iteration step, as P n ñ p n→ñ pñ →n = P n ñ p n →ñ pñ →n , where |n and |n are from basis 1 and |ñ from basis 2. A simple example is to collapse alternatingly to {Ŝ z i } and {Ŝ x i } eigenstates, respectively, for a spin-1/2 system as described in Ref. [22]. 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 [31] for each iteration step and lattice site [23]. Note that in both cases, these measurements break the symmetry associated with the conservation of total magnetization S z tot = iŜ z i or total particle numberN tot = in 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.

A. Maximally mixing bases
If we use a single collapse basis {|n } such as {Ŝ z i } eigenstates for a spin system or {n i } eigenstates for a system of bosons or fermions, there is no dynamics at all at infinite temperature (β = 0). Starting from an arbitrary initial state |n , transitions to all other basis states are impossible such that the METTS simulation is stuck in the state |n .
Having the METTS dynamics at high temperatures in mind, we can minimize autocorrelation times by switching between collapse bases {|n } and {|ñ } for which the distribution of overlaps | ñ|n | is as flat as possible. In other words, we want that all overlaps | ñ|n | are as small as possible. This guarantees that, at least at high temperatures, transitions to many states are possible and of similar probability. For β = 0 and a Hilbert space of dimension D [32], an optimal combination of bases yields transition probabilities p n→ñ = | ñ|n | 2 = 1/D ∀ n,ñ [33] and also p (2) n→n := ñ p n→ñ pñ →n = 1/D ∀ n,n . More generally, we call a sequence of K bases maximally mixing, if p (K) n→n = 1/D ∀ n,n at infinite temperature, where p (K) are the transition probabilities after K steps, i.e., when having cycled once through all K collapse bases.
We may call a sequence of K bases efficiently mixing if, at infinite temperature, there are many nonzero K-step transition probabilities p (K) n→n of comparable (small) amplitude for every 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.

B. Symmetric collapse bases with efficient mixing
Inspired by the discrete Fourier transform, one can construct maximally mixing bases for any D-dimensional Hilbert space. With a first orthonormal basis {|x | x = 1, . . . , D}, we can choose a second basis as which we call the Fourier basis with respect to {|x }.
Note that we are free to reorder the states |x . 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 {|σ } of {Ŝ z i } eigenstates for the magnetization M subspace H M , we can identify the states |σ with |x , where x = 1, . . . , dim H M , and obtain their Fourier basis {|k } according to Eq. (8). These two bases of H M are maximally mixing. We call the METTS scheme in which one alternates between collapses in this Fourier basis and the {Ŝ z i } eigenbasis "SF-Sz" (Symmetric Fourier w.r.t.Ŝ z i ). It is illustrated in Fig. 1.
However, such a global Fourier transform is technically infeasible because the basis states |k 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 {Ŝ z i } eigenbasis and can construct a Fourier basis for each symmetry sector of that block. This is scheme "SFb-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, . . . , +b/2 and corre- For each of these subspaces, we construct collapse bases according to Eq. (8) with |x referring in this case to the {Ŝ z i } eigenstates of some fixed magnetization on the block (in some i } eigenbasis on each site (Sz). Below, we alternate between the Sz basis and its global symmetric Fourier transform (8) (SF-Sz). These are maximally mixing. In practice, we approximate the SF-Sz scheme by restricting the basis construction to blocks of b sites, which are shifted by b/2 lattice sites for every second block collapse (SFb-Sz). Alternatively, a Haar-random symmetric basis can be chosen on these blocks (SRb-Sz). One can also omit the Sz-collapse in every second step and only use the Fourier basis, again shifting blocks in every second collapse (SFb). The swap collapse randomly partitions the lattice into pairs of sites, where, on each pair, the eigenstates of the swap operator (9) are chosen. The overall collapse basis is the tensor product of these. Again, we approximate the ideal swap collapse by restricting the pairing to blocks of b sites (SWAPb), where we shift by b/2 sites in every second collapse. 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, β = 0. With global maximally mixing bases such as the symmetric FourierŜ z i bases (SF-Sz), every second collapse, the up-spin jumps with equal probability to any new site. In the corresponding block scheme SFb-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 (β → ∞), every METTS is equal to the ground state in 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 "SRb-Sz".
For both the SFb-Sz and SRb-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 SFb and SRb.
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 SWAP = |↑↑ ↑↑| + |↑↓ ↓↑| + |↓↑ ↑↓| + |↓↓ ↓↓|.
Let SWAP i,j act on sites i and j of the lattice. As [SWAP i,j ,Ŝ z tot ] = 0, we can use a symmetric eigenbasis of the swap operator, 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 ("SWAPb" 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 {|n (κ) } (κ = 1, . . . , K) determines the METTS ensembles {|φ 4 + O(β) with probabilities 1 2 , leading to large statistical errors. Note that if our observable of interest happened to be, say,Ŝ 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 {|n } 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 n→n = | n |φ n | 2 = δ n,n . For small values of β, the probability that two subsequent samples are equal is still high because | n |φ n | 2 = δ n,n + O(β). 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 |φ n 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 [34][35][36] with Hamiltonian with varying values of the anisotropy parameter ∆. We focus on symmetric METTS computations for the canonical ensemble (CE) with zero magnetization and also compare the convergence speeds to those of nonsymmetric simulations of the grand-canonical ensemble (GCE) with Ŝ z tot = 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 [14], which has been adapted to also describe canonical ensembles and utilize symmetries [15,37], 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 (Error N 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 [38]. For the METTS simulations, the DMRG truncation weight and the Trotter time step ∆τ were chosen such that the corresponding errors are negligible compared to the statistical error. See, e.g., Ref. [23] for details or Ref. [3] for a general review.
In Fig. 3, the first two columns show results for the isotropic antiferromagnet (∆ = 1 in Eq. (10)) at inverse temperatures β = 2 and 8, and the third column shows results for ∆ = 3 in the gapped Néel phase at inverse temperature β = 8. The rows correspond to different observables: the energy per site Ĥ /L and the correlators Ŝ +

0Ŝ
− 1 and Ŝ + 0Ŝ − 3 . Site i = 0 is located at the center of the lattice, i.e., correlators are evaluated near the center to minimize finite-size effects [39]. At high temperatures (β = 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 Ŝ + 0Ŝ − 3 , 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 Ŝ + 0Ŝ − 3 , while it is broader for our novel block bases SFb and SRb. Still, the latter reduce autocorrelations between samples: The curves corresponding to our new block bases roughly follow the 1/ √ 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 or-  der 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 ∆, the Fourier bases (SFb) 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 {Ŝ z i } and {Ŝ x i } 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 α tuned such that Q gc α = 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 (Rb) 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. [30], METTS have so far been used exclusively to simulate grand-canonical ensembles (2) [21-24, 40, 41], 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 (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.

A. 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 |n with some quantum number Q, we use imaginary-time evolution to obtain the corresponding symmetric METTS sample |φ n , 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 |n only has contributions from a small number of different symmetry sectors. We write the new basis state as a superposition of (normalized) symmetry eigenstates, where |n j is the component with quantum number Q j . Subsequently, the |n j are evolved in imaginary time separately, again exploiting symmetries.
The DMRG time evolution [28,29] usually entails a Trotter decomposition. After each time step s∆τ → (s + 1)∆τ , 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 |n j consistent, we have to determine their relative weights P n j = n j |e −βĤ |n j . These can be obtained by multiplying the renormalization factors r As indicated, we do not actually multiply the r (j) s but rather accumulate the sum of their logarithms because the P n j decay exponentially in the system size and inverse temperature β. 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 −∆τĤ .
With the resulting normalized states |φ n j , the normalized METTS sample is then given by and |φ n j = P −1/2 n j e −βĤ/2 |n j . However, note that we can evaluate any observableÔ without explicitly encoding the superposition |φ n as an MPS according to φ n |Ô|φ n = kj c * k c j φ n k |Ô|φ n j . For symmetryconserving observables [Ô,Q] = 0 we have φ n k |Ô|φ n j = |n tDMRG |φ n ∝ e −βĤ/2 |n φ n |Ô|φ n |n = j α j |n j p n → n = | n |φ n | 2 sy m m e tr y -b re a k in g |n 1 |n 2 |n 3 |φ n 1 |φ n 2 |φ n 3 |φ n = j c j |φ n j φ n |Ô|φ n |n p n → n sy m m e tr ic tDMRG |φ n ∝ e −βĤ/2 |n φ n |Ô|φ n |n = j α j |n j p n → n sy m m e tr y -b re a k in g |n 1 |n 2 |n 3 |φ n 1 |φ n 2 |φ n 3 |φ n = j c j |φ n j φ n |Ô|φ n Figure 5. Illustration of our symmetric METTS algorithm for the GCE using small sets of symmetry eigenstates. An initial symmetric basis state |n is evolved in imaginary time to obtain the METTS |φn . After the evaluation of observables, a suitable symmetry-breaking collapse results in a non-symmetric basis state |n . It is split into its symmetric components {|n j }, which are separately evolved to obtain a small collection of symmetric states {|φ n j }. Observables are evaluated with respect to this collection according to Eq. (13), and a collapse using a symmetric basis yields a new symmetric basis state |n . We keep alternating between symmetry-breaking and symmetric collapse bases to explore different symmetry sectors according to their weights in the GCE.
0 for k = j as Q k = Q j ∀ k = j, and the expression reduces to the simple sum Subsequently, the small collection of symmetric METTS is collectively collapsed to a new single symmetric basis state |n , using, e.g., any of the symmetric collapse bases described in Sec. V. The transition probabilities for this projective measurement can be obtained usingÔ = |n n | in Eq. (13). We end up with a single state |n 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.

B. 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 {Ŝ z i } eigenbasis. For odd iteration steps, we randomly select n x sites on which theŜ x i eigenbasis is used as the collapse basis, while theŜ 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, {|n j } in Eq. (11), of the non-symmetric basis states |n 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 |n have n x + 1 components with different magnetizations that are separately evolved in imaginary time. Correspondingly, the total magnetizations of subsequent symmetric samples (|φ n and |φ n in Fig. 5) differ by ∆M ∈ {−n x , . . . , 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Ŝ x i eigenbasis above.

C. Quantum number trajectories and convergence
We test the algorithm for the isotropic spin-1/2 Heisenberg antiferromagnet, corresponding to ∆ = 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Ŝ x i eigenbasis is used for n x randomly selected sites and theŜ 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 (ν) = φ n ν |Ŝ z tot |φ n ν 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. [15]. The rows in Fig. 6 show data for different values of n x in the symmetry-breaking  10)) with Ŝ z tot = 0 and system size L = 64. The density plots show the distribution {pM } of the total magnetization M , computed with MPDOs as described in Ref. [15]. Lines show exemplary trajectories of M (ν) in Markov chains of symmetric METTS simulations using the SF4-Sz/Sx collapse bases as discussed in Sec. VIII. They are characterized by nx, the number of sites on which the symmetry-breakingŜ x i eigenbasis is used. Note that for even iteration steps ν, the state is an eigenstate ofŜ z tot with quantum number M (ν) while, for odd ν, we use the expectation value M (ν) = φnν |Ŝ z tot |φnν .
collapse. The columns correspond to different inverse temperatures β = 1, 4, and 8. Whenever we use the symmetric collapse basis (even steps ν in the Markov chain), the METTS sample is anŜ z tot eigenstate with magnetization M (ν). For odd ν, M (ν) is theŜ z tot expectation value. Note that for n x = L, we alternate between the {Ŝ x i } 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 β → ∞, 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 = 0. However, in this case, the imaginary-time evolution essentially projects onto the ground state, such that all contributions with M = 0 are exponentially suppressed. The parameter n x determines the speed with which the Markov chain can explore the symmetry sectors. It lim-its the maximum change in magnetization that can be achieved in a single symmetry-breaking collapse.
Comparing with quasi-exact purification data [15], Fig. 7 shows the error of the algorithm as a function of the number of samples. The DMRG truncation weight and the Trotter step-size ∆τ were again chosen such that the error is dominated by the statistical error. See Ref. [23] 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 Ĥ /L and the correlators Ŝ + at the highest temperature (β = 1), errors increase with increasing n x , and, at the lower temperatures (β = 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 {Ŝ z i } and {Ŝ 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.