According to the concept of typicality, an ensemble average can be accurately approximated by an expectation value with respect to a single pure state drawn at random from a high-dimensional Hilbert space. This random-vector approximation, or trace estimator, provides a powerful approach to, e.g. thermodynamic quantities for systems with large Hilbert-space sizes, which usually cannot be treated exactly, analytically or numerically. Here, we discuss the finite-size scaling of the accuracy of such trace estimators from two perspectives. First, we study the full probability distribution of random-vector expectation values and, second, the full temperature dependence of the standard deviation. With the help of numerical examples, we find pronounced Gaussian probability distributions and the expected decrease of the standard deviation with system size, at least above certain system-specific temperatures. Below and in particular for temperatures smaller than the excitation gap, simple rules are not available.
Methods such as the finite-temperature Lanczos method (FTLM) , , , , , ,  that rest on trace estimators , , , , , , , , ,  and thus – in more modern phrases – on the idea of typicality , , , , approximate equilibrium thermodynamic observables with very high accuracy , . In the canonical ensemble, the observable can be evaluated either with respect to a single random vector ,
or with respect to an average over R random vectors,
where numerator and denominator are averaged with respect to the same set of random vectors. The components of with respect to an orthonormal basis are taken from a Gaussian distribution with zero mean (Haar measure , , ), but in practice other distributions work as well. T, β and denote the temperature, inverse temperature and the Hamiltonian, respectively.
In this work, we discuss the accuracy of (1) and (2), where we particularly focus on the dependence of this accuracy on the system size or, to be more precise, the dimension of the effective Hilbert space spanned by thermally occupied energy eigenstates. While it is well established that the accuracy of both equations increases with the square root of this dimension, we shed light on the size dependence from two less studied perspectives. First, we study the full probability distribution of random-vector expectation values, for the specific example of magnetic susceptibility and heat capacity in quantum spin systems on a one-dimensional lattice. At high temperatures, our numerical simulations unveil that these distributions are remarkably well described by simple Gaussian functions over several orders of magnitude. Moreover, they clearly narrow with the inverse square root of the Hilbert-space dimension towards a δ function. Decreasing temperature at fixed system size, we find the development of broader and asymmetric distributions. Increasing the system size at fixed temperature, however, distributions become narrow and symmetric again. Thus, the mere knowledge of the standard deviation turns out to be sufficient to describe the full statistics of random-vector expectation values–at least at not too low temperatures.
The second central perspective of our work is taken by performing a systematic analysis of the scaling of the standard deviation with the system size, over the entire range of temperature and in various quantum spin models including spin-1/2 and spin-1 Heisenberg chains, critical spin-1/2 sawtooth chains, as well as cuboctahedra with spins 3/2, 2 and 5/2. We show a monotonous decrease of the standard deviation with increasing effective Hilbert-space dimension, as long as the temperature is high compared to some system-specific low-energy scale. Below this scale, the scaling can become unsystematic if only a very few low-lying energy eigenstates contribute. However, when averaging according to (2) over a decent number (∼ 100) of random vectors, one can still determine the thermodynamic average very accurately in all examples considered by us. A quite interesting example constitutes the critical spin-1/2 sawtooth chain, where a single state drawn at random is enough to obtain this average down to very low temperatures.
This paper is organized as follows. In Section 2 we briefly recapitulate models, methods, as well as typicality-based estimators. In Section 3 we present our numerical examples both for frustrated and unfrustrated spin systems. The paper finally closes with a summary and discussion in Section 4.
In this article we study several spin systems at zero magnetic field. They are of finite size and described by the Heisenberg model,
where the sum runs over ordered pairs of spins. Here and in the following operators are marked by a tilde, i.e. denotes the spin-vector operator at site i. denotes the exchange interaction between a spin at site i and a spin at site j. With the sign convention in (3), corresponds to antiferromagnetic interaction.
Numerator and denominator of (2), the latter is the partition function, are evaluated using a Krylov-space expansion, i.e. a spectral representation of the exponential in a Krylov space with as starting vector of the Krylov-space generation, compare , . One could equally well employ Chebyshev polynomials , ,  or integrate the imaginary-time Schrödinger equation with a Runge–Kutta method , , , the latter is used later in this paper as well.
If the Hamiltonian possesses symmetries, they can be used to block-structure the Hamiltonian matrix according to the irreducible representations of the employed symmetry groups , , which yields for the partition function
labels the subspace that belongs to the irreducible representation γ, NL denotes the dimension of the Krylov space and is the n-th eigenvector of in this Krylov space grown from . The energy eigenvalue is . To perform the Lanczos diagonalization for larger system sizes, we use the public code spinpack , .
In our numerical studies we evaluate the uncertainty of a physical quantity by repeating its numerical evaluation NS times. For this statistical sample we define the standard deviation of the observable in the following way:
is either evaluated according to (1) (m = r) or to (2) (m=FTLM), depending on whether the fluctuations of approximations with respect to one random vector or with respect to an average over R vectors shall be investigated.
We consider two physical quantities, the zero-field susceptibility as well as the heat capacity. Both are evaluated as variances of magnetization and energy, respectively, i.e.
We compare our results with the well-established high-temperature estimate
Here E0 denotes the ground-state energy. In general the prefactor α depends on the specific system, its structure and size, as well as on temperature , , but empirically often turns out to be a constant of order α ≈ 1 for high enough temperatures, compare , , . Rigorous error bounds, see Refs. , , share the dependence on , but lead to a prefactor that can be substantially larger than the empirical finding.
3 Numerical Results
We now present our numerical results. First, in the following Section 3.1, the full probability distribution of random-vector expectation values is discussed for shorter spin chains, where this distribution can be easily obtained by generating a large set of different random vectors. In the remainder of Section 3 the size dependence of the standard deviation is investigated for longer spin chains of spin and s = 1, which are treated by Lanczos methods. The interesting behaviour of a quantum critical delta chain is studied as well. Finally, we discuss the dependence of the standard deviation on the spin quantum number for a body of fixed size, the cuboctahedron.
3.1 Distribution of Random-Vector Expectation Values for Smaller Antiferromagnetic Spin-1/2 Chains
As a first step, in order to judge the accuracy of the single-state estimate in (1), it is instructive to study its full probability distribution p, obtained by drawing many [here ] random vectors. To be more precise, we evaluate the numerator of (1) for different random states , while its denominator is calculated as the average over all ,
The advantage of using this equation, instead of (1), is that the mean coincides with (2), the latter should be used to correctly obtain the low-temperature average in system of finite size . However, at sufficiently high temperatures or in sufficiently large systems, one might equally well use (1), as we have checked.
The single results for (9) are then collected into bins of appropriate width in order to form a “smooth” distribution p. While one might expect that p will be approximately symmetric around the respective thermodynamic average, the width of the distribution indicates how reliable a single random vector can approximate the ensemble average.
In this section, we study the probability distribution p (in the following denoted as pχ and pC) for the quantities and , and exemplarily consider the one-dimensional spin-1/2 Heisenberg model with antiferromagnetic nearest-neighbour coupling J > 0 and chain length N. Note that, as discussed in the upcoming Sections 3.2–3.5, details of the model can indeed have an impact on the behaviour of p in certain temperature regimes. Note further, that we focus in this section on small to intermediate system sizes N ≤ 20, where p can be easily obtained by generating a large set of different random vectors and evolving these vectors in imaginary time by, e.g. a simple Runge–Kutta scheme. We have checked that the Runge–Kutta scheme employed in this section has practically no impact on p.
To begin with, in Figure 1(a), pχ is shown for different chain lengths at infinite temperature . For all values of N shown here, we find that pχ is well described by a Gaussian distribution  over several orders of magnitude. While the mean of these Gaussians is found to accurately coincide with the thermodynamic average , we moreover observe that the width of the Gaussians becomes significantly narrower upon increasing N. This fact already visualizes that the accuracy of the estimate in (1) improves for increasing Hilbert-space dimension. In particular, as shown in the inset of Figure 1, the standard deviation scales as , where is the dimension of the Hilbert space. This is in agreement with (8) for α ≈ 1.2 and at β = 0. Note that as pχ is found to be a Gaussian, the width is sufficient to describe the whole distribution (apart from the average).
To proceed, Figure 1(b) again shows the probability distribution pχ, but now for the finite temperature . There are two important observations compared to the previous case of . First, for small N = 12, one clearly finds that pχ now takes on an asymmetric shape and the tails are not described by a Gaussian anymore. Importantly, however, upon increasing the system size N, pχ becomes narrower and eventually turns into a Gaussian again. One may speculate about possible reasons for the observed asymmetry: It might reflect an asymmetry of the distribution, which is already present at and small N, and then increases with increasing β; or it might also stem from the boundedness (positivity) of the observables, although the bounds are still far away for the presented case of in Figure 1(b). While this asymmetry remains to be explored in more detail in future work, it is expected that the Gaussian shape breaks down for small dimensions of the effective Hilbert space . It is worth pointing out that, even for very large dimensions, the very outer tails of the distribution are expected to be of non-Gaussian nature . Yet, these tails are hard to resolve in our numerical simulations, as a huge number of samples would be needed.
As a second difference compared to , we find that although pχ becomes narrower for larger N also at , this scaling is now considerably slower as a function of dimension d (see inset of Figure 1(b)). This is caused by the smaller effective Hilbert-space dimension at . As a consequence, for a fixed value of N, the single-state estimate in (1) becomes less reliable at compared to . However, let us stress that accurate calculations are still possible at T > 0 as long as N is sufficiently large. (Recall, that N ≤ 20 was chosen to be able to generate a large set of random vectors.)
In order to analyze the development of the probability distribution with respect to temperature in more detail, Figure 2(a) shows pχ for various values of in the range , for a fixed small system size N = 12. Note that the qualitative behaviour in principle is independent of N, but better to visualise for small N with a broader pχ. Starting from the high-temperature limit , we find that the maximum of pχ gradually shifts towards smaller values upon decreasing T.
This shift of the maximum is clearly visualised also in Figure 2(b), which shows the same data, but in a different style. Moreover, Figure 2(b) additionally highlights the fact that the probability distribution pχ for a fixed (and small) value of N becomes broader (and asymmetric) for intermediate values of T. Note, that pχ might become narrower again for smaller values of T, see also discussion in Sections 3.2 and 3.3.
Eventually, in Figures 3 and 4, we present analogous results for the full probability distribution p, but now for the heat capacity . Overall, our findings for pC are very similar compared to the previous discussion of pχ. Namely, we find that at , pC is very well described by Gaussians over several orders of magnitude. Moreover, the standard deviation again scales as at . As shown in Figure 3(b) and also in Figure 4, the emerging asymmetry of the probability distribution at small N and finite T is found to be even more pronounced for the heat capacity compared to the previous results for . Interestingly, we find that the maximum of pC, on the contrary, displays only a weak dependence on temperature (at least for the values of shown in Figure 4 – naturally, it is expected to change at lower temperatures and will go to zero at temperature T = 0).
3.2 Larger Antiferromagnetic Spin-1/2 Chains
Using a Krylov-space expansion one can nowadays reach large system sizes of for spins , see e.g. . But as we also perform a statistical analysis we restrict calculations to N ≤ 36 spins.
Following the scaling behaviour of as well as , which is shown in Figures 1 and 3, one expects a very narrow distribution of both quantities for N = 36 compared to e.g. N = 20, as the dimension is times bigger for N = 36, which yields a 256 times narrower distribution. Such a distribution is smaller than the linewidth in a plot.
That the distributions are narrow can be clearly seen by eye inspection in Figure 5 where the light blue curves depict thermal expectation values according to (1). For they fall on top of each other and coincide with the average over R = 100 realisations. Below this temperature the distributions widen, which is magnified by the fact that the real physical quantities susceptibility and heat capacity contain factors of and , respectively.
Their standard deviation is provided in Figure 6. Coming from high temperatures, the typical behaviour (8) switches to a behaviour that in general depends on system (here chain) and size below a characteristic temperature (here ). Nevertheless, the qualitative expectation that the standard deviation shrinks with increasing system size is met down to , below which no definite statement about the dependence on system size can be made. We conjecture that with growing N the increasing density of low-lying states as well as the vanishing excitation gap between singlet ground state and triplet first excited state influence the behaviour at very low temperatures strongly.
3.3 Antiferromagnetic Spin-1 Chains
In order to monitor an example where a vanishing excitation gap cannot be expected, not even in the thermodynamic limit, we study spin-1 chains that show a Haldane gap , , see Figure 7. The scaling formula (8) indeed suggests that for the standard deviations of the larger system with N = 24 should exceed those of the smaller system with N = 20, compare crossing curves of the estimator in Figure 8. However, the actual simulations show that this is not the case. The low-temperature fluctuations in the gap region are smaller for the larger system, at least for the two investigated system sizes.
3.4 Critical Spin-1/2 Delta Chains
As the final one-dimensional example we investigate a delta chain (also called sawtooth chain) in the quantum critical region, i.e. thermally excited above the quantum critical point (QCP) , , . The QCP is met when the ferromagnetic nearest-neighbour interaction J1 and the antiferromagnetic next-nearest neighbour interaction J2 between spins on adjacent odd sites assume a ratio of . At the QCP the system features a massive ground-state degeneracy due to multi-magnon flat bands as well as a double-peak density of states , , . Moreover, the typical finite-size gap is virtually not present at the QCP .
As the QCP does not depend on the size of the system and the structure of the analytically known multi-magnon flat band energy eigenstates does not either, we do not expect to find large finite-size effects when investigating the standard deviation of observables, e.g. of the heat capacity. It turns even out that by eye inspection no fluctuations are visible in Figure 9(a). The figure shows thermal expectation values (1) that virtually fall on top of each other. This means that a single random vector provides the equilibrium thermodynamic functions for virtually all temperatures. When evaluating the standard deviation, dashed curves in Figure 9(b), it turns out that it is unusually small, even for very low temperatures. The estimator (8) to which we compare had to be scaled in this case which might have two reasons. One reason could be that the large ground state degeneracy cannot be fully captured by the Krylov-space expansion and thus the evaluation of the estimator (1) by means of Eq. (4) is inaccurate. The other reason could be that the empirical finding of α ≈ 1 is not appropriate in this special case of a quantum critical system. However, the general rule that trace estimators are more accurate in larger Hilbert spaces is also observed here. The standard deviation of the smaller delta chain with N = 28 is a few times larger than for N = 32.
The result is an impressive example for what it means that a quantum critical system does not possess any intrinsic scale in the quantum critical region , . The only available scale is temperature. This means in particular that the low-energy spectrum is dense and therefore does not lead to any visible fluctuations of the estimated observables.
3.5 Antiferromagnetic Cuboctahedra with Spins 3/2, 2 and 5/2
Our last scaling analysis differs from the previous examples. The cuboctahedron is a finite-size body, that is equivalent to a kagome lattice with N = 12 , , . The structure is shown in Figure 10(b). Here, we vary the spin quantum number, not the size of the system. The dimension of the respective Hilbert spaces grows considerably which leads to the expected scaling (8) above temperatures of . But the low-temperature behaviour, in particular of the heat capacity for temperatures below the crossing of the estimators, eludes the expected order of more accurate results, i.e. smaller fluctuations for larger Hilbert-space dimension.
While the low-temperature behaviour and the standard deviation of the susceptibility are largely governed by the energy gap between singlet ground state and triplet excited state, and this does not vary massively with the spin quantum number, the heat capacity is subject to stronger changes (Fig. 11). When going from smaller to larger spin quantum numbers the strongly frustrated spin system loses some of its characteristic quantum properties while becoming more classical with increasing spin s. In particular, the low-lying singlet states below the first triplet state which dominate the low-temperature heat capacity move out of the gap for larger spin s , .
It may thus well be that the type of Hilbert space enlargement, due to growing system size which leads to the thermodynamic limit or growing spin quantum number which leads to the classical limit, is important for the behaviour of the estimators (1) and (2) at low temperatures.
4 Discussion and Conclusions
To summarize, we have studied the finite-size scaling of typicality-based trace estimators. In these approaches, the trace over the high-dimensional Hilbert space is approximated by either (i) a single random state or (ii) the average over a set of R random vectors. In particular, we have focussed on the evaluation of thermodynamic observables such as the heat capacity and the magnetic susceptibility for various spin models of Heisenberg type. Here, the temperature dependence of these quantities has been generated by means of a Krylov-space expansion where the random states are used as a starting vector for the expansion.
As a first step, we have studied the full probability distribution of expectation values evaluated with respect to single random states. As an important result, we have demonstrated that for sufficiently high temperatures and large enough system sizes (i.e. sufficiently large effective Hilbert-space dimension Zeff), the probability distributions are very well described by Gaussians . In particular, for comparatively high temperatures, our numerical analysis has confirmed that the standard deviation of the probability distribution scales as , and that this width already describes the full distribution.
In contrast, for lower temperatures, we have shown that (i) the probability distributions can become non-Gaussian and (ii) the scaling of can become more complicated and generally depends on the specific model and observable under consideration. While a larger Hilbert-space dimension often leads to an improved accuracy of the random-state approach at low temperatures as well, compare the investigation on kagome lattice antiferromagnets of sizes N = 30 and N = 42 in , we have also provided examples where this expectation can break down for too small Zeff, compare also .
A remarkable example is provided by the spin-1/2 sawtooth chain with coupling-ratio . Due to the (virtually) gapless dense low-energy spectrum at the quantum critical point, we have found that statistical fluctuations remain negligible throughout the entire temperature range with only minor dependence on system size (see also Ref.  for a similar finding in a spin-liquid model).
In conclusion, we have demonstrated that typicality-based estimators provide a convenient numerical tool in order to accurately approximate thermodynamic observables for a wide range of temperatures and models. While in some cases, even a single pure state is sufficient, the accuracy of the results can always be improved by averaging over a set of independently drawn states.
This work was supported by the Deutsche Forschungsgemeinschaft DFG (397067869 (STE 2243/3-1); 355031190 (FOR 2692); 397300368 (SCHN 615/25-1)). Computing time at the Leibniz Center in Garching is gratefully acknowledged. All authors thank Hans De Raedt, Peter Prelovšek, Patrick Vorndamme, Peter Reimann, Jochen Gemmer as well as Katsuhiro Morita for valuable comments.
 P. Prelovšek, in: The Physics of Correlated Insulators, Metals, and Superconductors (Eds. E. Pavarini, E. Koch, R. Scalettar, and R. M. Martin), Forschungszentrum Jülich GmbH Zentralbibliothek, Verlag Jülich, Jülich, Germany 2017.Search in Google Scholar
 G. H. Golub and U. von Matt, Tikhonov regularization for large scale problems, Tech. Rep. Stanford University 1997 technical report SCCM-97-03.Search in Google Scholar
 K. Inoue, Y. Maeda, H. Nakano, and Y. Fukumoto, IEEE Transact. Magn. 55, 1 (2019).Search in Google Scholar
 J. Schulenburg, Spinpack 2.56, Magdeburg University (2017).Search in Google Scholar
 T. Vojta, Ann. Phys. 9, 403 (2000).10.1002/1521-3889(200006)9:6<403::AID-ANDP403>3.0.CO;2-RSearch in Google Scholar
 A. Honecker and M. E. Zhitomirsky, J. Phys.: Conf. Ser. 145, 012082 (2009).Search in Google Scholar
© 2020 Walter de Gruyter GmbH, Berlin/Boston