A more complete analysis of the standard deviation dependence on the number of random samples R is shown in Figure for two different Krylov subspace dimensions, 10 and 20 respectively: this reveals a typical Monte Carlo 1/ √ R scaling where the dimension of the Krylov subspace (i.e. the number of Lanczos quadrature points) has very little effects unless unresonably small subspaces are employed. We note at this point that Lanczos quadrature is extremely efficient in approximating the terms reported in Eq.( ) as it requires half of the points, for given target accuracy, compared to other quadrature strategies such as the ones based on Chebyshev polynomials. This, together with a rigorous mathematical error analysis of the SL trace estimator, is discussed in the original work of Ubaru, Chen and Saad. Figure shows the evolution of the stochastic trace estimation with the number of random samples for a single point estimation where the trace approaches the target value with a standard deviation magnitude becoming smaller with the increase of the number of random samples, in line to what shown in Figure . A further aspect to probe is the standard deviation's dependency on the system size. This is shown in Figure for increasingly larger argon clusters. The trend shows how the error decreases with the increase of the system size thus denoting a desirable statistical self-averaging of errors, in line with the fact that, in the thermodynamic limit, fluctuations die out. We note in passing that a similar behavior was observed in stochastic approaches to compute electronic correlation energies. Further words on this aspect will be spent later on in this section. It is now of interest to analyze the scaling with system size of the proposed SL trace estimation algorithm and to compare it with the standard approach requiring the full eigendecomposition of V. The building of the Krylov subspace only requires matrix-vector products which are performed in a linear-scaling fashion in contrast to the standard full-diagonalization approach which is cubic scaling with the system size as shown in Figure . We stress that the plot was generated by keeping the number of random samples constant, however, in virtue of the self- averaging property of the algorithm the number of random samples could have been notably reduced (for a given target accuracy) along with the system size and this would have contributed to an even better scaling of the SL algorithm's CPU-time with system size.