We describe a numerical algorithm for approximating the equilibrium-reduced density matrix and the effective (mean force) Hamiltonian for a set of system spins coupled strongly to a set of bath spins when the total system (system + bath) is held in canonical thermal equilibrium by weak coupling with a “super-bath”. Our approach is a generalization of now standard typicality algorithms for computing the quantum expectation value of observables of bare quantum systems via trace estimators and Krylov subspace methods. In particular, our algorithm makes use of the fact that the reduced system density, when the bath is measured in a given random state, tends to concentrate about the corresponding thermodynamic averaged reduced system density. Theoretical error analysis and numerical experiments are given to validate the accuracy of our algorithm. Further numerical experiments demonstrate the potential of our approach for applications including the study of quantum phase transitions and entanglement entropy for long range interaction systems.

The equilibrium thermodynamics of quantum systems is a growing area of research,1–3 and in many cases, the systems of interest are open, for instance, because the system of interest is a subsystem of some closed system. This presents a practical barrier to the study of open quantum systems since open quantum systems enjoy many complexities not present in closed systems. Indeed, many thermodynamic properties of even the simplest open quantum systems may be seemingly anomalous, e.g., the specific heat or entropy may be negative.4–6 

In practice, only a relatively small number of quantum systems have Hamiltonians with known analytic diagonalizations, and while numerical techniques exist for certain systems, for instance, a small system coupled to a harmonic reservoir,7–11 efficient algorithms for the case when the total system is comprised entirely of spins are seemingly less available. This paper aims to address this gap by introducing and analyzing an algorithm for computing equilibrium thermodynamic properties of small open quantum spin systems coupled arbitrarily to baths consisting of a small to moderate number of spins.

Consider a total system Hamiltonian Ht of the form

Ht=H̄s+H̄b+Hsb,
(1)

where H̄s=HsIb corresponds to the Hamiltonian of the bare system, H̄b=IsHb corresponds to the Hamiltonian of the bare bath, and Hsb an interaction term accounting for non-negligible interactions between the system and bath.

Throughout, we assume the total system (system + bath) is in a canonical equilibrium state due to weak contact with a “super-bath” that holds the total system at inverse Boltzmann temperature β=(kBT)1. The states of the total system/bare system/bare bath are then described by the density matrices

ρt/s/b(β)=exp[βHt/s/b]Zt/s/b(β),
(2)

where the partition functions are computed by

Zt/s/b(β)=tr(exp[βHt/s/b]).
(3)

Due to system–bath entanglement, the density matrix ρs(β) for the bare system does not describe the equilibrium state of the system when strongly coupled to the bath; i.e., when Hsb ≠ 0t.6 Instead, one must start with the total system density and “trace out” the effects of the bath. The resulting reduced system density matrix, sometimes called the mean force Gibbs state, is given by

ρ*(β)=trb(ρt(β)).
(4)

Here, trb(·) is the partial trace with respect to the bath.

The reduced system density matrix ρ*(β) can be expressed in terms of an effective Hamiltonian H*(β), often called the Hamiltonian of the mean force, by the relation

ρ*(β)=exp[βH*(β)]Z*(β),
(5)

where the corresponding partition function is

Z*(β)=tr(exp[βH*(β)])=Zt(β)Zb(β).
(6)

Thus, the Hamiltonian of the mean force has an explicit formula

H*(β)=1βlntrb(exp[βHt])Zb(β).
(7)

We turn readers to Ref. 6 for a more detailed discussion on the Hamiltonian of the mean force and mean force Gibbs state.

The partition function Z*(β) for the reduced system provides access to thermodynamic quantities of the reduced system. These quantities include the heat capacity, magnetization, susceptibility, etc., and are often functions of the Helmholtz free energy F*(β) = −β−1 ln(Z*(β)) and can, therefore, be written in terms of the difference of the corresponding quantities of the bare system and bare bath. In fact, since Z*(β) only depends on Zt(β) and Zb(β), any quantity depending on Z*(β) can be computed using what have become standard numerical techniques; see Sec. II D.

Other properties of the reduced system cannot alone be derived from the partition functions of the bare system and bath. Indeed, Zt(β) and Zb(β) do not account for the structure of interactions between the system and bath and, therefore, cannot express any quantity, such as the von Neumann entropy,12 which may depend on entanglement between the system and bath. In these cases, one must obtain the Hamiltonian of the mean force H*(β) or the corresponding density matrix ρ*(β). Computing such quantities numerically is the main focus of this paper.

We denote by Hs and Hb the Hilbert spaces for the system and bath, so that Ht=HsHb is the Hilbert space for the total system. For a Hilbert space H, we denote by |H| the dimension of the space and by L(H) the set of self-adjoint linear operators on H. Throughout, |i⟩ is the ith standard basis vector of dimension determined by context and It/s/b and 0t/s/b are identity and zero operators, respectively, on Ht/s/b.

Broadly, quantum typicality refers to the idea that, in many cases, a random state is representative of the overall state of a system. Early concepts of typicality were hinted at by Schrödinger13 and proved rigorously by von Neumann;14 see Ref. 15 for an overview. Mathematically, the notion of typicality can be viewed as concentration of a random variable about its expectation value.

Before we describe the notion of typicality on which our algorithm is based, we introduce a similar, yet mathematically equivalent, form of typicality asserting that the quantum expectation value ⟨v|O|v⟩ of an observable O in a state |v⟩ is overwhelmingly likely to be near to the quantum expectation value of O, at least when |v⟩ is chosen randomly from a suitable ensemble.

More precisely, suppose OL(Ht) is an observable of the total system and |vHt is a random state sampled from the uniform distribution on the set of all states. Then, the “Hilbert space average” (HA) of the density matrix |v⟩⟨v| (with respect to the above distribution) is

HA[|vv|]=|Ht|1It.
(8)

Therefore, using basic properties of the trace and Hilbert space average, the quantum expectation value of O when the system is in state |v⟩ satisfies

HAv|O|v=HAtr(|vv|O)
(9)
=tr(HA[|vv|]O)
(10)
=|Ht|1tr(O).
(11)

Here, |Ht|1tr(O) is the quantum expectation of O when the state of the system is described by the density matrix |Ht|1It.

More generally, if the state of the system is described by an arbitrary density matrix ρ, then the quantum expectation of an observable O is tr(ρO)=tr(ρOρ). Therefore, in order to find an ensemble |ω⟩ so that the quantum expectation value ⟨ω|O|ω⟩ has Hilbert space average equal to tr(ρO), we simply define |ω⟩ by |ω|Ht|ρ|v, where |v⟩ remains a uniformly chosen random state. Early analyses16–18 showed that ⟨ω|O|ω⟩ concentrates about its quantum expectation tr(ρO) by bounding the variance of ⟨ω|O|ω⟩ and applying Chebyshev’s inequality. Subsequent analyses show that ⟨ω|O|ω⟩ is in fact sub-Gaussian and concentrates far more sharply about tr(ρO) than suggested by Chebyshev’s inequality. Mathematically precise bounds are discussed in Sec. IV.

We now introduce the version of typicality on which our algorithm is based. Recall that the partial trace of ρO with respect to the bath Hilbert space can be expressed as

trb(ρO)=i=1|Hb|(Isi|)ρO(Is|i).
(12)

Thus, if |vHb is a random state chosen uniformly from all states, we have HA[|vv|]=|Hb|1I, and it is relatively straightforward to see that

HA(Isv|)|Hb|ρO|Hb|ρ(Is|v)=trb(ρO).
(13)

Indeed, expand

|v=i=1|Hb|i|v|i
(14)

and observe that, for any i,j=1,2,,|Hb|,

HAi|vj|v=|Hb|1,i=j,0,ij.
(15)

Then, using the linearity of the Hilbert space average and (15), we find that, for any ALHt,

HA(Isv|)A(Is|v)
(16)
=i=1|Ht|j=1|Ht|HAi|vj|v(Isi|)A(Is|j)
(17)
=i=1|Hb||Hb|1(Isi|)A(Is|i)
(18)
=|Hb|1trb(A).
(19)

In fact, a simple union bound shows that (Isv|)|Hb|ρO|Hb|ρ(Is|v) also exhibits sub-Gaussian concentration about its mean. From a linear algebraic perspective, this is equivalent to previously defined versions of typicality. Even so, we were unable to find explicit reference to this phenomenon in the literature.

For concreteness, we will consider total systems consisting of N interacting spins sites of spin number s in a magnetic field of strength h pointing in the z-direction. The corresponding Heisenberg Hamiltonian for the total system is

Ht=i,j=1NJi,jσiσj+h2i=1Nσiz,
(20)

where we have used the shorthand

Ji,jσiσj=Ji,jxσixσjx+Ji,jyσiyσjy+Ji,jzσizσjz,
(21)

where Ji,jx/y/z describes the coupling strength between sites i and j in the x/y/z coordinate directions. Here, σix/y/z gives the component spin operator for the ith spin site and acts trivially on the Hilbert spaces associated with other spin sites but as the (2s + 1) × (2s + 1) component spin matrix σx/y/z on the ith spin site. In matrix form, σix/y/z is written

σix/y/z=IIi1 termsσx/y/zIIN(i1) terms.
(22)

Without loss of generality, we will take the system to be the first Ns sites and the bath to be the remaining Nb = NNs sites. Let Is={1,,Ns} and Ib={Ns+1,,N}. Ht can then be decomposed according to (1) where

H̄s=i,jIsJi,jσiσj+h2iIsσiz,
(23)
H̄b=i,jIbJi,jσiσj+h2iIbσiz,
(24)
Hsb=iIs,jIbjIs,iIbJi,jσiσj.
(25)

From the structure of the total system Hamiltonian (20) considered in this paper, we can derive the limits for the mean force Hamiltonian and mean force Gibbs state at high and low temperature. We summarize these limits here and provide proofs in  Appendix A.

In the high temperature limit, as β → 0, it is straightforward, although a bit tedious, to verify that H*(β)=Hs+trb(Hsb)/tr(Ib)+O(β). Since Hsb accounts for interactions between the system and bath (and, therefore, contains no interactions between spin sites within the bath), we have that trb(Hsb) = 0s. This implies that

limβ0H*(β)=Hs.
(26)

In this case, ρ*(β)(2s+1)NsIs.

In the low temperature limit, as β → ∞, we have that ρt(β) is convergent to r1i=1r|ψiψi|, where {|ψi⟩} are the r ground states for the total system This implies that ρ*(β) is also convergent to a fixed density matrix. Specifically,

limβρ*(β)=1ri=1rtrb(|ψiψi|).
(27)

Then, since H*(β) = −β−1 ln[Z*(β)ρ*(β)], it is easy to verify that

limβH*(β)=limβln(Z*(β))I=(EtEb)I,
(28)

where Et and Eb are the ground state energies of Ht and Hb, respectively.

A major challenge to the design of algorithms for quantum spin systems is that problem sizes grow exponentially with the total system size. Perhaps the simplest numerical approach is to apply an exact eigensolver to Ht in order to compute exp(−βHt). However, this quickly becomes intractable unless there are many symmetries present in the system. Moreover, even if exp(−βHt) could somehow be computed efficiently, the cost of storing it would be very high. For instance, if the total system is comprised of 20 spin-12 particles, exp(−βHt) is a matrix with 220 × 220 entries, which would require nearly 8.8 Tbyte of memory to store in a 64 bit double precision format.19 

Over the past several decades, typicality based approaches, such as the finite temperature Lanczos method (FTLM), have become among the most widely used numerical methods for approximating equilibrium thermodynamic properties of closed quantum systems.20–26 Such methods make use of Krylov subspace methods to avoid explicitly forming matrix functions such as exp(−βHt) and instead compute quantities like ⟨v| exp(−βHt)|v⟩. Moreover, because Krylov subspace methods are “matrix-free,” they only access Ht through matrix-vector products. Theoretical analysis of such algorithms is an active area of research.27–29 

Krylov subspace methods can also be used to numerically compute the mean force Hamiltonian and reduced system density matrix at high and low temperature using the limits from Sec. II C. Indeed, at high temperature ρ*(β) is trivial and H*(β) converges to Hs. At low temperature, ρ*(β) depends only on the ground state(s) of Ht and H*(β) depends only on the ground state energies of Ht and Hb.

We emphasize that the task of computing a partial trace of an explicit matrix, even naively, is not particularly difficult. More efficient approaches to this task have also been studied.30 However, to the best of our knowledge, the task of computing the partial trace of implicit matrices such as exp(−βHt), without ever explicitly constructing the matrix, has not been thoroughly studied.

In this section, we describe a numerical method for computing the mean force Hamiltonian and reduced system density matrix. In the case of an empty system, so that Ht=Hb, everything in this section reduces to the well-known stochastic Lanczos quadrature algorithm.28,29,31 A theoretical error analysis is given in Sec. IV.

Let ALHt and recall that (Is ⊗ ⟨v|)A(Is ⊗ |v⟩) is an unbiased estimator for |Hb|1trb(A). To improve the estimator’s accuracy, we can average multiple copies. Specifically, given nv iid samples {|vj⟩} of |v⟩ (with HA[|v⟩⟨v|] = Ib), we can define the averaged estimator

1nvi=1nv(Isvj|)A(Is|vj).
(29)

Physically, we can view this averaging technique as constructing a new total system containing multiple independent copies of the original total system. This point of view has its origins in Gibbs’s 1902 concept of the statistical ensemble.32 While many analyses from physics15,17,18 show that a high dimensional Hilbert space |Hb| can lead to quantum typicality under regular conditions given in their models, in numerical analysis, it is more common to consider the error of estimator (29) when nv is large. We believe the mathematical correspondence between |Hb| and nv is worthy of a more rigorous mathematical analysis. For instance, generalized fundamental thermodynamic relations were recently unraveled by replacing thermodynamic infinite-size limit (|Hb|) with multiple-measurement limit (nv → ∞).33 

Often, A = f[H] for some function f:RR and HLHt; for instance, f = x ↦  exp(−βx) and H = Ht. If H has known diagonalization, then we can easily compute (Is ⊗ ⟨v|)f[H](Is ⊗ |v⟩). Unfortunately, diagonalizing large Hermitian matrices is often exceedingly expensive. In fact, even for highly structured matrices, such as those considered in this paper, exact diagonalization may be too costly.

A natural approach to avoiding such costs is to apply the block Lanczos algorithm, Algorithm 1, to H and Is ⊗ |v⟩ for k iterations to obtain a k|Hs|×k|Hs| block-tridiagonal matrix T. Given a matrix V with VV = Is, Algorithm 1 computes orthonormal matrices {Qj}, j = 1, …, k + 1 such that QiQj=δi,jIs and for all jk,

span{V,HV,,HjV}=span{Q1,Q2,,Qj}.
(30)

These vectors satisfy a symmetric block-tridiagonal recurrence

HQ=QT+Qk+1BkEk,
(31)

where Ek = |k⟩ ⊗ Is and

T=A1B1B1Bk1Bk1Ak,Q=|||Q1Q2Qk|||.
(32)
ALGORITHM 1.

Block Lanczos.

1: procedure BLOCK-LANCZOS(H, V, k
2: Q1 = V
3: forj = 1, 2, …, kdo 
4: Z=HQjQj1Bj1 
5: Aj=QjZ 
6: Z = ZQjAj 
7: Qj+1, Bj = qr(Z
8: return {Qj}, {Aj}, {Bj
1: procedure BLOCK-LANCZOS(H, V, k
2: Q1 = V
3: forj = 1, 2, …, kdo 
4: Z=HQjQj1Bj1 
5: Aj=QjZ 
6: Z = ZQjAj 
7: Qj+1, Bj = qr(Z
8: return {Qj}, {Aj}, {Bj

The expression (Is ⊗ ⟨v|)f[H](Is ⊗ |v⟩) can be approximated by

(1|Is)f[T](|1Is).
(33)

Assuming k|Hs| is small enough so that T can be diagonalized exactly, f[T] can be computed directly. This is a block Gauss quadrature approximation and is exact if f is a polynomial of degree at most 2k − 1.34 See Sec. 6.6 for details.

To obtain our final estimator, we apply this approach to each term in (29). Specifically, for j = 1, …, nv, denoting by Tj the resulting block-tridiagonal matrix obtained by block Lanczos run on H and Is ⊗ |vj⟩, our estimator for trb(f[H]) is

1nvj=1nv(1|Is)f[Tj](|1Is).
(34)

1. Costs

The accuracy and computational cost of (34) depend on both nv and k. Application of the block Lanczos method to each term of (34) uses k block matrix-vector products with blocks of |Hs| vectors along with O(|HtHs|)=O(|HbHs|2) storage. While a block matrix-vector product requires |Hs| times as many operations as a standard inner product, in practice, block matrix-vector products can often be computed nearly as quickly as a single matrix-vector product. If each term in (34) is computed sequentially, the computational cost of computing (34) scales linearly with nv while the storage cost is constant. On the other hand, all terms can be computed entirely in parallel at the cost of nv times more storage.

Like other typicality algorithms, our approach requires matrix-products with the total system Hamiltonian Ht and is, therefore, limited to systems for which it is possible to store several vectors of size |Ht|. Our algorithm is also limited in terms of the system size, since we must store a number of vectors of length |Ht| proportional to the system dimension |Hs|. Even so, our algorithm is significantly more efficient than exact diagonalization based techniques.

To compute H*(β), we can use the same test states {|vj⟩} for estimating the partial trace and the trace over the bath. Specifically, if (Tt)j is the block-tridiagonal matrix produced by the block Lanczos algorithm with Ht and (Is ⊗ |vj⟩) and (Tb)j is the tridiagonal matrix produced by Lanczos with Hb and |vj⟩, by (7) we obtain the estimators

H*(β)1βlnj=1nv(Isvj|)exp[βHt](Is|vj)j=1nvvj|exp[βHb]|vj
(35)
1βlnj=1nv(1|Is)exp[β(Tt)j](|1Is)j=1nv1|exp[β(Tb)j]|1.
(36)

Of course, only the latter estimator is practically computable.

When β is large and the eigenvalues of (Tt)j or (Tb)j are negative, then the computation of the exponential may overflow. To avoid this, we can simply replace x ↦ exp(−βx) with x ↦ exp(−β(xE0)), where E0 is chosen to be smaller than the smallest eigenvalues of (Tt)j and (Tb)j; for instance, it is chosen to be (an estimate of) the smallest eigenvalue of Ht.

It is worth noting that if we use (36) or (35) to compute an approximation to ρ*(β), by (5) we obtain exactly the estimators

ρ*(β)j=1nv(Isvj|)exp[βHt](Is|vj)trj=1nv(vj|Is)exp[βHt](|vjIs)
(37)
j=1nv(1|Is)exp[β(Tt)j](|1Is)trj=1nv(1|Is)exp[β(Tt)j](|1Is).
(38)

Thus, if one requires the eigenvalues of H*(β) and ρ*(β), we suggest first computing the eigenvalues of j=1nv(1|Is)exp[β(Tt)j](|1Is) and then transforming them to obtain the eigenvalues of H*(β) and ρ*(β). This avoids the need for the computation of the matrix logarithm.

In this section, we discuss bounds that provide intuition on how to balance nv and k. For the block-size one case, such bounds have been studied extensively; see Ref. 29 for a recent review.

Tail bounds trace estimators were studied in Refs. 35 and 36, with more recent and refined analyses given in Refs. 3739; see Ref. 29 for historical context. For constants ɛ > 0 and δ ∈ (0, 1), and BL(Hb), these analyses aim to bound the number of samples nv required so that

Probtr(B)|Hb|1nvj=1nvvj|B|vj>ε<δ.
(39)

The resulting bounds are typically simple functions depending on the distribution of |v⟩, the value of nv, and basic properties of B, such as its operator norm, Frobenius norm, or dimension. Roughly speaking, the analyses in Refs. 29 and 40 imply that, for small ɛ, it suffices to take nv=O(B22Hb|1ε2)ln(2/δ). Here, Õ is equivalent to O but with poly-logarithmic factors in the constituent parameters suppressed for readability; i.e., we say, a variable is Õ(h(t)) if, for some k ≥ 0, the variable is O(log(t)kh(t)).

We can leverage such results to provide similar bounds for our partial trace estimator. Toward this end, decompose ALHt as

A=m,n=1|Hs||mn|Am,n,
(40)

where Am,nL(Hb). Fix ɛ > 0 and δ ∈ (0, 1), and suppose that for all m, n, we have chosen nv so that

Probtr(Am,n)|Hb|1nvj=1nvvj|Am,n|vj>ε̃<δ̃
(41)

for some ε̃ and δ̃, the exact values of which will soon become apparent. Applying a union bound over all pairs m, n, we obtain the bound

Probm,n:tr(Am,n)|Hb|1nvj=1nvvj|Am,n|vj>ε̃<|Hs|2δ̃.
(42)

Next, note that

trb(A)=m,n=1|Hs|tr(Am,n)|mn|
(43)

and that

m|(Isvj|)A(Is|vj)|n=(m|vj|)A(|n|vj)=vj|Am,n|vj.
(44)

For any XL(Hs), we have that X|Hs|maxm,nm|X|n. Putting everything together, we find that

Probtrb(A)|Hb|1nvj=1nv(Isvj|)A(Is|vj)>ε<δ
(45)

if we take ε̃=ε/|Hs| and δ̃=δ/|Hs|2. This allows existing bounds for standard trace estimators to be easily carried over to partial trace estimation. Like the basic trace estimators, since (29) is the average of nv iid samples, nv must scale like O(ε2).

1. A note on high temperatures

At high temperature, our estimator will compute trb(exp(−βHt)) efficiently using only a single sample. This is because the randomness in our sample state is averaged over many states of the system, thereby reducing the variance of the output.

This does not necessarily imply an accurate estimate to H*(β) using a single sample. Indeed, recall that H*(β)=Hs+trb(Hsb)/tr(Ib)+O(β). Thus, at high temperature, the approximation (36) differs from Hs by an additive factor nv1j=1nv(vj|Is)Hsb(|vjIs)/tr(Ib), which should be corrected for. Since HA[(⟨vj| ⊗ Is)Hsb(|vj⟩ ⊗ Is)] = 0s, it may be possible to use this difference as a rough indicator of the accuracy of the averaged partial trace estimator.

2. A note on low temperatures

At low temperature, exp(−βHt) is dominated by the ground state, or a few states near the ground state. As such, our approach will have relatively high variance because the randomness in our sample state is averaged over only several states. Thus, while the estimator still provides an unbiased estimate for trb(exp(−βHt)), many samples are required. In this setting, other approaches such as low-rank approximation or hybrid methods make more sense.38,40–43

Recall that Vp(A)V=E1p[T]E1 for all polynomials p of degree at most 2k − 1. For brevity, let Error=Vf[H]VE1f[T]E1. Then, for any polynomial p with deg(p) ≤ 2k − 1,

Error=Vf[H]VVp(A)V+E1p[T]E1E1f[T]E1
(46)
V(f[H]p[H])V+E1(p[T]f[T])E1
(47)
f[H]p[H]+p[T]f[T]
(48)
2maxx[λmin,λmax]|f(x)p(x)|.
(49)

Here, we have used that ‖V‖ ≤ 1 since VV = I. Optimizing over p with deg(p) ≤ 2k − 1, we find that

Error2mindeg(p)2k1maxx[λmin,λmax]|f(x)p(x)|.
(50)

Note that while we could instead apply the regular Lanczos algorithm to each of the vectors in V individually, the resulting algorithm would only be exact for polynomials of degree at most k − 1.

Analytic functions such as the exponential can be approximated by polynomials of degree growing just logarithmically with the desired accuracy.44 This means that k typically does not need to be very large. Then, as long as |Hs| is also not too large, we can directly diagonalize each Tj to compute terms of (34), possibly exploiting the block-tridiagonal structure along the way.

1. Finite precision arithmetic

In finite precision arithmetic, the output of the (block) Lanczos algorithm may be significantly different than what would be obtained in exact arithmetic. In particular, the columns of Q may lose orthogonality. This has led to some hesitance to use Lanczos based approaches without using costly explicit reorthogonalization.21,23,28,45

Careful analysis of the Lanczos algorithm in finite precision arithmetic46,47 can be leveraged to show that the Lanczos algorithm still works well for the task of applying matrix functions to vectors48 and quadratic forms49 even in finite precision arithmetic. In effect, these analyses show that Lanczos performs at least as well as explicit polynomial methods (for instance, the kernel polynomial method23); see Ref. 29 for a discussion and comparison.

While we are aware of no similarly rigorous analyses for the block Lanczos algorithm, we believe it is reasonable that similar results hold, albeit with possibly worse dependencies on certain parameters. Our numerical experiments suggest the iterate E1f[T]E1 computed by the block Lanczos algorithm still provides a good approximation to Vf[H]V, even when orthogonality of the Lanczos vectors is lost. A rigorous analysis of the block Lanczos algorithm in floating point arithmetic is needed in order to make any definitive statement.

In this section, we provide several numerical examples to demonstrate the accuracy and flexibility of our approach. In all cases, we consider isotropic XY spin systems; i.e., Ji,jx=Ji,jy and Ji,jz=0.

We begin with the simple nearest neighbor spin chain with connection strength J. We set N = 18, take the system to be the 1st and 2nd spin in the chain (Ns = 2), and put the magnetic field strength at h = 0.3J. We then run our algorithm using k = 30 Lanczos iterations and nv = 100 samples.

Figure 1 shows the median, 10% quantile, and 90% quantile of 100 independent runs of our algorithm with the above parameters. Because the spin chain is solvable analytically, we are able to compare our results against the exact eigenvalues of ρ*(β) and H*(β), which were computed in Ref. 5 and are summarized in  Appendix C.

FIG. 1.

Eigenvalues of H*(β) and ρ*(β) for a spin chain of length 18 with system taken as the first two sites. Legend: Algorithm median (—), algorithm 10%–90% quantiles (gray filled rectangle), exact solution (⁃ ⁃), direct numerical computation of high and low temperature limits (-------).

FIG. 1.

Eigenvalues of H*(β) and ρ*(β) for a spin chain of length 18 with system taken as the first two sites. Legend: Algorithm median (—), algorithm 10%–90% quantiles (gray filled rectangle), exact solution (⁃ ⁃), direct numerical computation of high and low temperature limits (-------).

Close modal

Note that the quality of the approximation of the eigenvalues of ρ*(β) is accurate visually except at lower temperatures where there is higher variance (although the median outputs of our algorithm still agree very well with the true values). The higher variance at low temperature is expected based on our analysis above and could be decreased by increasing nv. However, at low temperature, the eigenvalues of ρ*(β) can be easily obtained by directly computing the ground state of the total system.

For all temperatures observed, we see a good agreement between our algorithm’s approximation to the eigenvalues of H*(β) and the exact eigenvalues of H*(β). As expected, at high temperature, we observe that the spectrum of H*(β) matches the spectrum of Hs, while at low temperature, the spectrum converges to a constant EtEb.

We now consider a spin ladder. The edge coupling strength is set to J, the run coupling strength is set to −0.45J, and the magnetic field strength is set to J. We set N = 20 and consider systems of size Ns = 2 of spins connected by a rung. There are ten such systems, although only five are unique due to symmetry.

In Fig. 2, we show the temperature dependence of the mean force Hamiltonian’s eigenvalues for these five choices of system computed using k = 30 and nv = 50. Because the bare system is the same in all cases, the high temperature behavior is the same. However, the low temperature limit and the qualitative behavior at intermediate temperature are different.

FIG. 2.

Eigenvalues of H*(β) for a spin ladder with ten rungs with various choices of system (shown as circled sites in depiction of system configuration above images). Legend: algorithm (—), direct numerical computation of high and low temperature limits (-------).

FIG. 2.

Eigenvalues of H*(β) for a spin ladder with ten rungs with various choices of system (shown as circled sites in depiction of system configuration above images). Legend: algorithm (—), direct numerical computation of high and low temperature limits (-------).

Close modal

We observe the eigenvalues of the mean force Hamiltonian appear to “cross,” implying the occurrence of degenerate energy levels in the effective Hamiltonian at certain temperatures. More precisely, if the system starts from high enough temperature, we can guarantee all energy levels are not degenerate. Then, when the temperature decreases to a certain point, degeneracy appears, but it immediately disappears right after passing that temperature. Eventually, those energy levels converge to a single level in the zero temperature limit. We call this phenomenon “temperature-induced degeneracy.” This phenomenon is impossible for weak interaction systems (Hsb = 0) since only strong interaction systems have a temperature-dependent effective Hamiltonian. We suggest that there might be a certain unspecified type of symmetry induced by strong interactions at certain temperature so that the degeneracy appears due to that symmetry. We believe that specifying that strong interaction-induced symmetry and its connection with degeneracy is worthy of a further work.

Next, we turn our attention to a spin chain with long range interactions in the presence of magnetic fields of varying strength. Specifically, for a spin chain with N = 16 spins, we take Ji,jx=Ji,jy=J|ij|α, α ≥ 0. While α = ∞ gives the solvable system studied above, to the best of our knowledge, this system is not exactly solvable for arbitrary α. Throughout, the system is taken as the first 2 spins and the bath as the latter 14.

1. von Neumann entropy

For this experiment, we set α = 1 and vary h from 0J to 2.2J. We run our algorithm with nv = 400 and k = 60 and note that each value of h requires an entirely independent run of our algorithm.

Phase plots showing the relationship between the system’s von Neumann entropy − tr(ρ*(β) ln(ρ*(β))), temperature, and magnetic field strength are given in Fig. 3. The zero temperature limit is computed by using a black box eigensolver to find the ground state of the total system (under the assumption of a single ground state). The relative smoothness between consecutive values of magnetic field strength provides some indication of the variance of the output produced by our algorithm.

FIG. 3.

Relationship between von Neumann entropy, temperature, and magnetic field strength in spin chain with long range power law interactions. Here, the system is taken as the first two spins and the bath as the remaining spins. Maximum value (at a fixed temperature) shows as black “×.”

FIG. 3.

Relationship between von Neumann entropy, temperature, and magnetic field strength in spin chain with long range power law interactions. Here, the system is taken as the first two spins and the bath as the remaining spins. Maximum value (at a fixed temperature) shows as black “×.”

Close modal

At zero temperature, the von Neumann entropy indicates the presence of a magnetic field-induced quantum phase transition50–53 from an entangled state to a pure state of the system of interest around h ≈ 1.4J. In the exactly solvable chain with α = ∞, a similar quantum phase transition from entangled to pure state occurs at h = 2J.5 This implies a dependence of the location of the quantum phase transition on the interaction range. This dependence is studied for a closely related model in Ref. 52.

Since zero temperature is not experimentally realizable, practical studies of quantum phase transitions require observations to be made at finite temperature.51,54 For α = 1, we find a region corresponding to low entanglement extending beyond T = 0 when h is sufficiently large. Despite some finite sample size noise in the output of our algorithm, it is clear that the staircase-like behavior observed at zero temperature extends to finite temperature. These same phenomena are present in the α = ∞ case, and a comparison of the α = 1 phase plots in Fig. 3 with the α = ∞ phase plots in Fig. 6 in the Appendix shows a dependence on the interaction range.

Based on the above finding, we believe that our algorithm, which enables the temperature dependence of the von Neumann entropy (and similar quantities) to be studied at finite temperature, has the potential to inform the study of quantum phase transitions.

In addition to the magnetic field-induced quantum phase transition, we also observe a temperature-induced phase transition by observing the strength of the magnetic field at which the von Neumann entropy is maximal (at a fixed temperature). The values of these maxima are shown as black “×” in Fig. 3 and show the emergence of a critical phenomenon: The maximal von Neumann entropy at high temperature is obtained at zero magnetic field strength. However, for temperatures below a critical temperature Tc ≈ 0.5J/β, the maximal von Neumann entropy occurs at nonzero magnetic field strength.

2. Deviation in internal energy

From (1), it is not instantly clear how to split the energy due to Hsb into the system of interest and the bath.55 Here, we adopt a difference of state functions (equilibrium averages of fluctuating observables) before and after coupling to answer this question. Before coupling, tr(Hsρs) is the equilibrium average of the internal energy of the system. On the other hand, after coupling, the relevant equilibrium average for the internal energy of the system is tr(H*(β)ρ*(β)). Therefore, the difference

tr(H*(β)ρ*(β))tr(Hsρs)
(51)

can be interpreted as the “deviation” in the system’s internal energy (state function) by coupling to the bath via Hsb; i.e. the interaction energy.

We note that the state functions of internal energy given in Eq. (51) are based on the approach by mean energy rather than the approach by partition function. For strongly coupled systems, these two approaches lead to different thermodynamic results,6,56–58 which are not within the scope of this work. On the other hand, the deviation defined by the difference in Eq. (51) is more relevant to the solvation energy generated in a transfer process of taking a solute from a vacuum to a solution.59 

To study this quantity, we set h = 0 and vary α, the parameter controlling the interaction range. We run our algorithm with nv = 100 and k = 60. As seen in Fig. 4, the deviation is dependent on α. Specifically, we observe that at low temperatures, shorter range interactions (larger α) correspond to a higher deviation, while at high temperatures, longer range interactions (lower α) correspond to a higher deviation. Moreover, we observe that the deviation of internal energy for longer range interactions changes sign from positive to negative in certain temperature regions.

FIG. 4.

Difference of tr(H*(β)ρ*(β)) and tr(Hsρs(β)) for a spin chain with varying interaction decay rates. In all cases, the system is taken as the first two spins and the bath as the remaining spins. In order, purple to yellow, α = 0, 0.1, 0.3, 0.5, 1, 2, ∞.

FIG. 4.

Difference of tr(H*(β)ρ*(β)) and tr(Hsρs(β)) for a spin chain with varying interaction decay rates. In all cases, the system is taken as the first two spins and the bath as the remaining spins. In order, purple to yellow, α = 0, 0.1, 0.3, 0.5, 1, 2, ∞.

Close modal

Borrowing the idea from solvation thermodynamics59 and hybridization processes,60 our observation suggests that the coupling process for the system with the longer range interactions with its coupled bath can be either an endothermic process (positive regions) or an exothermic process (negative regions), and it is determined by the temperature! This further implies that the system and its bath switches their interactions from being attractive to begin repelling or vice versa at certain critical temperatures. Our simulations also indicate the existence of a critical point αc ∈ (1, 2) such that deviation is negative at all temperatures if α > αc, i.e., the coupling process for the system is always exothermic when the range of interactions is too short.

We now study the effect of the coupling strength on the deviation of the long range spin chain from the previous example when α = 1 and h = 0. Specifically, we consider the Hamiltonian

Ht=H̄s+H̄b+εHsb,
(52)

where ɛ ≥ 0 determines the coupling strength and H̄s, H̄b, and Hsb are all as in the previous example. We again run our algorithm with nv = 100 and k = 60.

In Fig. 5, we observe that as the coupling strength decreases, the deviation also decreases. In particular, in the limit ɛ = 0, the deviation is zero. This aligns with the limits described in Sec. II C. While the deviation depends on ɛ, we observe a change in sign for all ɛ > 0. This suggests the existence of temperature-induced sign changes in the energy deviation depends on the range, rather than the strength, of interactions.

FIG. 5.

tr(H*(β)ρ*(β)) and tr(Hsρs(β)) for a spin chain with long range power law interactions and varying system-bath coupling strengths. In all cases, the system is taken as the first two spins and the bath as the remaining spins. In order, purple to yellow, ɛ = 0 to 1 in increments of 0.1.

FIG. 5.

tr(H*(β)ρ*(β)) and tr(Hsρs(β)) for a spin chain with long range power law interactions and varying system-bath coupling strengths. In all cases, the system is taken as the first two spins and the bath as the remaining spins. In order, purple to yellow, ɛ = 0 to 1 in increments of 0.1.

Close modal

We have introduced a numerical algorithm, based on the concept of partial typicality, for computing the mean force Hamiltonian and average reduced system density matrix of strongly coupled spin systems. Numerical experiments on a solvable system indicate that our algorithm can produce highly accurate results. This is complemented by a theoretical analysis of the behavior of the algorithm. Further experiments on a range of systems that are not exactly solvable demonstrate the flexibility and power of the algorithm. We hope that this algorithm will enable further study of the thermodynamic properties of open quantum systems.

There are a number of concrete directions for future work. First, it would be interesting to extend methods that balance low-rank approximation and stochastic trace estimation38,40,43 to the task of computing partial traces of matrix functions. In particular, the recent approach of Ref. 43 demonstrates the effectiveness of hybrid approaches for matrix function trace estimation problems in quantum physics. Generalizing the approach to computing partial traces would enable higher quality approximations at low temperature. It would also be interesting to study how variants of the finite temperature Lanczos method, such as the microcanonical Lanczos method (MCLM)61 and low temperature Lanczos method (LTLM),45 can be extended to the setting of this paper. Finally, we believe it should be relatively straightforward to use the same techniques used in this paper, typicality for partial traces and block Krylov subspace methods, to compute dynamical quantities.

The authors thank Michele Campisi, Hong Qian, Peter Talkner, Lowell Thompson, Yao Wang, Yijing Yan, and Ying-Jen Yang for feedback and suggestions.

This work was supported by the National Science Foundation under Grant No. DGE-1762114. Any opinions, findings, and conclusions, or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article.

1. High temperature limit

Expanding at β = 0, we find

lntrb(exp[βHt])
(A1)
=lntrb(IsIbβH+O(β2))
(A2)
=lntr(Ib)Isβtrb(H)+O(β2)
(A3)
=lnIsβtrb(Ht)tr(Ib)+O(β2)+ln[tr(Ib)Is]
(A4)
=βtrb(Ht)tr(Ib)+O(β2)+ln(tr(Ib))Is.
(A5)

Here, we have used the property that ln[AB] = ln[A] + ln[B] if A and B commute. Similarly, we can expand

ln(tr(exp[βHb]))=ln(tr(IbβHb+O(β2)))
(A6)
=βtr(Hb)tr(Ib)+O(β2)+ln(tr(Ib)).
(A7)

Thus, combining these expressions,

H*(β)=1βlntrb(exp[βHt])tr(exp[βHb])
=1βlntrb(exp[βHt])
(A8)
+1βlntr(exp[βHb])Is
(A9)
=trb(Ht)tr(Ib)tr(Hb)tr(Ib)Is+O(β)
(A10)
=Hs+trb(Hsb)tr(Ib)+O(β),
(A11)

where, in the final equality, we have used that

trb(Ht)=trb(HsIb+IsHb+Hsb)
(A12)
=tr(Ib)Hs+tr(Hb)Is+trb(Hsb).
(A13)

Now, recall that Hsb accounts for interactions between the system and bath and is therefore a linear combination of terms of the form

Iiσx/y/zIiIbIsIjσx/y/zIj
(A14)
=(Iiσx/y/zIi)(Ijσx/y/zIj),
(A15)

where i+i=(2s+1)Ns1 and j+j=(2s+1)Nb1. Applying basic properties of the partial trace, we see that

trb(Iiσx/y/zIi)(Ijσx/y/zIj)
(A16)
=trbIjσx/y/zIj(Iiσx/y/zIi)
(A17)
=trb(Ij)trb(σx/y/z)trb(Ij)(Iiσx/y/zIi)
(A18)
=0s.
(A19)

Here, we have used the fact that tr(σx/y/z) = 0 for any spin number s. From numerical computation of the reduced density matrix, the trace is linear, so we in fact have that trb(Hsb) = 0s.

We have, therefore, established that

limβ0H*(β)Hs.
(A20)

2. Low temperature limit

Write the eigenvalue decomposition of Ht as Ht=i=1|Ht|Ei|ψiψi| for orthogonal eigenvectors |ψi⟩. Then,

ρ(Ht)=exp[βHt]=i=1|Ht|exp(βEi)|ψiψi|i=1|Ht|exp(βEi).
(A21)

When β → ∞, only the terms corresponding to the ground state (smallest absolute eigenvalue) remain.

Note that

H*(β)=1βln[Z*(β)ρ*(β)]
(A22)
=1βln(Z*(β))Isβ1ln[ρ*(β)].
(A23)

Since ρ*(β) is convergent to a fixed density matrix, β−1 ln[ρ*(β)] → 0s as β → ∞. By definition,

1βln(Z*(β))=1βln(Zt(β)/Zb(β))
(A24)
=1βln(Zt(β))+1βln(Zb(β)).
(A25)

In the low temperature limit, the partition functions for the bare system and bare bath are dominated by their respective ground state energies so that

1βln(Zt/b(β))=1βln(tr(exp(βHt/b)))
(A26)
=λmax(Ht/b)+O(β1).
(A27)

We, therefore, have that

limβH*(β)=(EtEb)Is,
(A28)

where Et and Eb are the ground state energies of Ht and Hb, respectively.

Consider a random pure product state |vs⟩ ⊗ |vb⟩ where HA[|vs/b⟩⟨vs/b|] = Is/b. We have that

vs|(Isvb|)A(Is|vb)|vs=(vs|vb|)A(|vs|vb)
(B1)
=(Ibvs|)vb|A|vb|(Ibvs|).
(B2)

Thus, estimators of this form have the desirable property that

HA(vs|vb|)A(|vs|vb)|vs=vs|trb(A)|vs,
(B3)
HA(vb|vs|)A(|vb|vs)|vb=vb|trs(A)|vb.
(B4)

That is, we obtain approximations of the trace and partial trace that are consistent in the sense that our estimate of the trace of the partial traces is equal.

Expressions of the form |vs⟩ ⊗ |vb⟩ are called rank-one vectors and have been studied for the task of trace estimation with the goal of reducing the amount of randomness required.62 It is conceivable that there are situations in which using a pure product state |v⟩ = |v1⟩ ⊗⋯⊗ |vN⟩ would be desirable. While this produces unbiased estimators for arbitrary partial traces, the number of such estimators that must be averaged to reach a fixed error is exponential in N.63,64

We summarize the relevant quantities from Ref. 5. For k = 1, …, N, define

λk(N)=h2JcoskπN+1,
(C1)
Nk(N)=1+exp(βλk(N))1.
(C2)

Next, define

σ1xσ2x=4N+1k=1NsinkπN+1sin2kπN+1Nk(N),
(C3)
σjz=4N+1k=1NsinjkπN+12Nk(N),
(C4)
σ1zσ2z=σ1zσ2zσ1xσ2x2,
(C5)
δ=4σ1xσ2x2+(σ1zσ2z)2.
(C6)

For Ns = 2, the eigenvalues of ρ*(β) are given by

p1=(1+σ1z+σ2z+σ1zσ2z)/4,
(C7)
p2=(1δσ1zσ2z)/4,
(C8)
p3=(1+δσ1zσ2z)/4,
(C9)
p4=(1σ1zσ2z+σ1zσ2z)/4.
(C10)

Moreover, the partition function of a length N chain is given by

ZN=exp(βNh/2)k=1N1+exp(βλk(N)).
(C11)

We have Z*(β)=ZN/ZNb, and the eigenvalues of H*(β) are given by

hj=1βln(Z*(β)pj),j=1,2,3,4.
(C12)

Phase plots showing the relationship between the von Neumann entropy −tr(ρ*(β) ln(ρ*(β))), temperature, and magnetic field strength for the exactly solvable spin chain with N = 16 and the first two spins taken as the system are given in Fig. 6.

FIG. 6.

Relationship between von Neumann entropy, temperature, and magnetic field strength in spin chain with short range nearest neighbor interactions. Here, the system is taken as the first two spins and the bath as the remaining spins.

FIG. 6.

Relationship between von Neumann entropy, temperature, and magnetic field strength in spin chain with short range nearest neighbor interactions. Here, the system is taken as the first two spins and the bath as the remaining spins.

Close modal
1.
J.
Gemmer
,
M.
Michel
, and
G.
Mahler
, “
Quantum thermodynamics: Emergence of thermodynamic behavior within composite quantum systems
,” in
Lecture Notes in Physics
, 2nd ed. (
Springer
,
2009
).
2.
S.
Vinjanampathy
and
J.
Anders
, “
Quantum thermodynamics
,”
Contemp. Phys.
57
,
545
579
(
2016
).
3.
R.
Alicki
and
R.
Kosloff
, “
Introduction to quantum thermodynamics: History and prospects
,” in
Fundamental Theories of Physics
(
Springer International Publishing
,
2018
), pp.
1
33
.
4.
G.-L.
Ingold
,
P.
Hänggi
, and
P.
Talkner
, “
Specific heat anomalies of open quantum systems
,”
Phys. Rev. E
79
,
061105
(
2009
).
5.
M.
Campisi
,
D.
Zueco
, and
P.
Talkner
, “
Thermodynamic anomalies in open quantum systems: Strong coupling effects in the isotropic XY model
,”
Chem. Phys.
375
,
187
194
(
2010
).
6.
P.
Talkner
and
P.
Hänggi
, “
Colloquium: Statistical mechanics and thermodynamics at strong coupling: Quantum and classical
,”
Rev. Mod. Phys.
92
,
041002
(
2020
).
7.
N.
Makri
and
D. E.
Makarov
, “
Tensor propagator for iterative quantum time evolution of reduced density matrices. I. Theory
,”
J. Chem. Phys.
102
,
4600
4610
(
1995
).
8.
N.
Makri
and
D. E.
Makarov
, “
Tensor propagator for iterative quantum time evolution of reduced density matrices. II. Numerical methodology
,”
J. Chem. Phys.
102
,
4611
4618
(
1995
).
9.
A.
Gelzinis
and
L.
Valkunas
, “
Analytical derivation of equilibrium state for open quantum system
,”
J. Chem. Phys.
152
,
051103
(
2020
).
10.
A. S.
Trushechkin
,
M.
Merkli
,
J. D.
Cresser
, and
J.
Anders
, “
Open quantum system dynamics and the mean force Gibbs state
,”
AVS Quantum Sci.
4
(
1
),
012301
(
2022
).
11.
Y.-F.
Chiu
,
A.
Strathearn
, and
J.
Keeling
, “
Numerical evaluation and robustness of the quantum mean force Gibbs state
,”
Phys. Rev. A
106
,
012204
(
2022
).
12.
J.
von Neumann
and
R.
Beyer
,
Mathematical Foundations of Quantum Mechanics, Goldstine Printed Materials
(
Princeton University Press
,
1955
).
13.
E.
Schrödinger
, “
Energieaustausch nach der Wellenmechanik
,”
Ann. Phys.
388
,
956
968
(
1927
).
14.
J.
von Neumann
, “
Beweis des Ergodensatzes und desH-Theorems in der neuen Mechanik
,”
Z. Phys.
57
,
30
70
(
1929
); arXiv:abs/1003.2133 (in English).
15.
S.
Goldstein
,
J. L.
Lebowitz
,
C.
Mastrodonato
,
R.
Tumulka
, and
N.
Zanghì
, “
Normal typicality and von Neumann’s quantum ergodic theorem
,”
Proc. R. Soc. A
466
,
3203
3224
(
2010
).
16.
P.
Reimann
, “
Typicality for generalized microcanonical ensembles
,”
Phys. Rev. Lett.
99
,
160404
(
2007
).
17.
C.
Bartsch
and
J.
Gemmer
, “
Dynamical typicality of quantum expectation values
,”
Phys. Rev. Lett.
102
,
110403
(
2009
).
18.
S.
Sugiura
and
A.
Shimizu
, “
Thermal pure quantum states at finite temperature
,”
Phys. Rev. Lett.
108
,
240401
(
2012
).
19.

While Ht may be sparse, the matrix exponential is not, in general, sparse.

20.
J.
Skilling
, “
The eigenvalues of mega-dimensional matrices
,” in
Maximum Entropy and Bayesian Methods
(
Springer Netherlands
,
1989
), pp.
455
466
.
21.
J.
Jaklič
and
P.
Prelovšek
, “
Lanczos method for the calculation of finite-temperature quantities in correlated systems
,”
Phys. Rev. B
49
,
5065
5068
(
1994
).
22.
R.
Schnalle
and
J.
Schnack
, “
Calculating the energy spectra of magnetic molecules: Application of real- and spin-space symmetries
,”
Int. Rev. Phys. Chem.
29
,
403
452
(
2010
).
23.
A.
Weiße
,
G.
Wellein
,
A.
Alvermann
, and
H.
Fehske
, “
The kernel polynomial method
,”
Rev. Mod. Phys.
78
,
275
306
(
2006
).
24.
J.
Schnack
,
J.
Richter
, and
R.
Steinigeweg
, “
Accuracy of the finite-temperature Lanczos method compared to simple typicality-based estimates
,”
Phys. Rev. Res.
2
,
013186
(
2020
).
25.
H.
Schlüter
,
F.
Gayk
,
H.-J.
Schmidt
,
A.
Honecker
, and
J.
Schnack
, “
Accuracy of the typicality approach using Chebyshev polynomials
,”
Z. Naturforsch., A
76
,
823
834
(
2021
).
26.
F.
Jin
,
D.
Willsch
,
M.
Willsch
,
H.
Lagemann
,
K.
Michielsen
, and
H.
De Raedt
, “
Random state technology
,”
J. Phys. Soc. Jpn.
90
,
012001
(
2021
).
27.
I.
Han
,
D.
Malioutov
,
H.
Avron
, and
J.
Shin
, “
Approximating spectral sums of large-scale matrices using stochastic Chebyshev approximations
,”
SIAM J. Sci. Comput.
39
,
A1558
A1585
(
2017
).
28.
S.
Ubaru
,
J.
Chen
, and
Y.
Saad
, “
Fast estimation of tr(f(A)) via stochastic Lanczos quadrature
,”
SIAM J. Matrix Anal. Appl.
38
,
1075
1099
(
2017
).
29.
T.
Chen
,
T.
Trogdon
, and
S.
Ubaru
, “
Randomized matrix-free quadrature for spectrum and spectral sum approximation
,” arXiv:2204.01941 [math.NA] (
2022
).
30.
J.
Maziero
, “
Computing partial traces and reduced density matrices
,”
Int. J. Mod. Phys. C
28
,
1750005
(
2017
).
31.
Z.
Bai
,
G.
Fahey
, and
G.
Golub
, “
Some large-scale matrix computation problems
,”
J. Comput. Appl. Math.
74
,
71
89
(
1996
).
32.
J.
Gibbs
,
Elementary Principles in Statistical Mechanics: Developed with Especial Reference to the Rational Foundation of Thermodynamics
(
C. Scribner’s Sons
,
1902
).
33.
Z.
Lu
and
H.
Qian
, “
Emergence and breaking of duality symmetry in generalized fundamental thermodynamic relations
,”
Phys. Rev. Lett.
128
(
15
),
150603
(
2022
).
34.
G.
Golub
and
G.
Meurant
,
Matrices, Moments and Quadrature with Applications
, Princeton Series in Applied Mathematics (
Princeton University Press
,
2009
).
35.
H.
Avron
and
S.
Toledo
, “
Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix
,”
J. ACM
58
,
1
34
(
2011
).
36.
F.
Roosta-Khorasani
and
U.
Ascher
, “
Improved bounds on sample size for implicit matrix trace estimators
,”
Found. Comput. Math.
15
,
1187
1212
(
2014
).
37.
A.
Cortinovis
and
D.
Kressner
, “
On randomized trace estimates for indefinite matrices with an application to determinants
,”
Found. Comput. Math.
22
,
875
(
2021
).
38.
R. A.
Meyer
,
C.
Musco
,
C.
Musco
, and
D. P.
Woodruff
, “
Hutch++: Optimal stochastic trace estimation
,” in
Symposium on Simplicity in Algorithms (SOSA)
(
Society for Industrial and Applied Mathematics
,
2021
), pp.
142
155
.
39.
T.
Chen
,
T.
Trogdon
, and
S.
Ubaru
, “
Analysis of stochastic Lanczos quadrature for spectrum approximation
,” in
Proceedings of the 37th International Conference on Machine Learning
[
Proceedings of Machine Learning Research (PMLR)
,
2021
]; arXiv:2105.06595 [cs.DS].
40.
D.
Persson
,
A.
Cortinovis
, and
D.
Kressner
, “
Improved variants of the Hutch++ algorithm for trace estimation
,”
SIAM J. Matrix Anal. Appl.
43
(
3
),
1162
1185
(
2022
).
41.
L.
Lin
, “
Randomized estimation of spectral densities of large matrices made accurate
,”
Numer. Math.
136
,
183
213
(
2016
).
42.
A. S.
Gambhir
,
A.
Stathopoulos
, and
K.
Orginos
, “
Deflation as a method of variance reduction for estimating the trace of a matrix inverse
,”
SIAM J. Sci. Comput.
39
,
A532
A558
(
2017
).
43.
T.
Chen
and
E.
Hallman
, “
Krylov-aware stochastic trace estimation
,” arXiv:2205.01736 [math.NA] (
2022
).
44.
L. N.
Trefethen
,
Approximation Theory and Approximation Practice
, Extended Edition (
SIAM
,
2019
).
45.
M.
Aichhorn
,
M.
Daghofer
,
H. G.
Evertz
, and
W.
von der Linden
, “
Low-temperature Lanczos method for strongly correlated systems
,”
Phys. Rev. B
67
,
161103
(
2003
).
46.
C. C.
Paige
, “
Error analysis of the Lanczos algorithm for tridiagonalizing a symmetric matrix
,”
IMA J. Appl. Math.
18
,
341
349
(
1976
).
47.
C. C.
Paige
, “
Accuracy and effectiveness of the Lanczos algorithm for the symmetric eigenproblem
,”
Linear Algebra Appl.
34
,
235
258
(
1980
).
48.
C.
Musco
,
C.
Musco
, and
A.
Sidford
, “
Stability of the Lanczos method for matrix function approximation
,” in
Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA’18
(
Society for Industrial and Applied Mathematics
,
2018
), pp.
1605
1624
.
49.
L. A.
Knizhnerman
, “
The simple Lanczos procedure: Estimates of the error of the Gauss quadrature formula and their applications
,”
Comput. Math. Math. Phys.
36
,
1481
1492
(
1996
).
50.
A. W.
Rost
,
R. S.
Perry
,
J.-F.
Mercure
,
A. P.
Mackenzie
, and
S. A.
Grigera
, “
Entropy landscape of phase formation associated with quantum criticality in Sr3Ru2O7
,”
Science
325
,
1360
1363
(
2009
).
51.
T.
Werlang
,
C.
Trippe
,
G. A. P.
Ribeiro
, and
G.
Rigolin
, “
Quantum correlations in spin chains at finite temperatures and quantum phase transitions
,”
Phys. Rev. Lett.
105
,
095702
(
2010
).
52.
T.
Koffel
,
M.
Lewenstein
, and
L.
Tagliacozzo
, “
Entanglement entropy for the long-range Ising chain in a transverse field
,”
Phys. Rev. Lett.
109
,
267203
(
2012
).
53.
O.
Breunig
,
M.
Garst
,
A.
Klümper
,
J.
Rohrkamp
,
M. M.
Turnbull
, and
T.
Lorenz
, “
Quantum criticality in the spin-1/2 Heisenberg chain system copper pyrazine dinitrate
,”
Sci. Adv.
3
,
eaao3773
(
2017
).
54.
S. L.
Sondhi
,
S. M.
Girvin
,
J. P.
Carini
, and
D.
Shahar
, “
Continuous quantum phase transitions
,”
Rev. Mod. Phys.
69
,
315
333
(
1997
).
55.
C.
Jarzynski
, “
Stochastic and macroscopic thermodynamics of strongly coupled systems
,”
Phys. Rev. X
7
,
011008
(
2017
).
56.
M. F.
Gelin
and
M.
Thoss
, “
Thermodynamics of a subensemble of a canonical ensemble
,”
Phys. Rev. E
79
,
051121
(
2009
).
57.
U.
Seifert
, “
First and second law of thermodynamics at strong coupling
,”
Phys. Rev. Lett.
116
,
020601
(
2016
).
58.
J.-T.
Hsiang
and
B.-L.
Hu
, “
Quantum thermodynamics at strong coupling: Operator thermodynamic functions and relations
,”
Entropy
20
,
423
(
2018
).
59.
A. Y.
Ben-Naim
,
Solvation Thermodynamics
(
Springer Science and Business Media
,
2013
).
60.
H.
Gong
,
Y.
Wang
,
H.-D.
Zhang
,
Q.
Qiao
,
R.-X.
Xu
,
X.
Zheng
, and
Y.
Yan
, “
Equilibrium and transient thermodynamics: A unified dissipaton-space approach
,”
J. Chem. Phys.
153
,
154111
(
2020
).
61.
M. W.
Long
,
P.
Prelovšek
,
S.
El Shawish
,
J.
Karadamoglou
, and
X.
Zotos
, “
Finite-temperature dynamical correlations using the microcanonical ensemble and the Lanczos algorithm
,”
Phys. Rev. B
68
,
235106
(
2003
).
62.
Z.
Bujanovic
and
D.
Kressner
, “
Norm and trace estimation with random rank-one vectors
,”
SIAM J. Matrix Anal. Appl.
42
,
202
223
(
2021
).
63.
R.
Vershynin
, “
Concentration inequalities for random tensors
,”
Bernoulli
26
,
3139
(
2020
).
64.
S.
Bamberger
,
F.
Krahmer
, and
R.
Ward
, “
The Hanson-Wright inequality for random tensors
,” arXiv:2106.13345 [math.PR] (
2021
).