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.

## I. INTRODUCTION

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 *H*_{t} of the form

where $H\u0304s=Hs\u2297Ib$ corresponds to the Hamiltonian of the bare system, $H\u0304b=Is\u2297Hb$ corresponds to the Hamiltonian of the bare bath, and *H*_{sb} 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 $\beta =(kBT)\u22121$. The states of the total system/bare system/bare bath are then described by the density matrices

where the partition functions are computed by

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 *H*_{sb} ≠ 0_{t}.^{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

Here, tr_{b}(·) 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

where the corresponding partition function is

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

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 *Z*_{t}(*β*) and *Z*_{b}(*β*), 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, *Z*_{t}(*β*) and *Z*_{b}(*β*) 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.

### A. Notation

We denote by $Hs$ and $Hb$ the Hilbert spaces for the system and bath, so that $Ht=Hs\u2297Hb$ 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 *i*th standard basis vector of dimension determined by context and *I*_{t/s/b} and 0_{t/s/b} are identity and zero operators, respectively, on $Ht/s/b$.

## II. BACKGROUND

### A. Typicality

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ödinger^{13} 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 $O\u2208L(Ht)$ is an observable of the total system and $|v\u3009\u2208Ht$ 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

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

Here, $|Ht|\u22121tr(O)$ is the quantum expectation of *O* when the state of the system is described by the density matrix $|Ht|\u22121It$.

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(\rho O)=tr(\rho O\rho )$. 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 $|\omega \u3009\u2254|Ht|\rho |v\u3009$, where |*v*⟩ remains a uniformly chosen random state. Early analyses^{16–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

Thus, if $|v\u3009\u2208Hb$ is a random state chosen uniformly from all states, we have $HA[|v\u3009\u3008v|]=|Hb|\u22121I$, and it is relatively straightforward to see that

Indeed, expand

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

Then, using the linearity of the Hilbert space average and (15), we find that, for any $A\u2208LHt$,

In fact, a simple union bound shows that $(Is\u2297\u3008v|)|Hb|\rho O|Hb|\rho (Is\u2297|v\u3009)$ 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.

### B. Spin systems

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

where we have used the shorthand

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

Without loss of generality, we will take the system to be the first *N*_{s} sites and the bath to be the remaining *N*_{b} = *N* − *N*_{s} sites. Let $Is={1,\u2026,Ns}$ and $Ib={Ns+1,\u2026,N}$. *H*_{t} can then be decomposed according to (1) where

### C. High and low temperature limits

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*(\beta )=Hs+trb(Hsb)/tr(Ib)+O(\beta )$. Since *H*_{sb} accounts for interactions between the system and bath (and, therefore, contains no interactions between spin sites within the bath), we have that tr_{b}(*H*_{sb}) = 0_{s}. This implies that

In this case, $\rho *(\beta )\u2192(2s+1)\u2212NsIs$.

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

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

where *E*_{t} and *E*_{b} are the ground state energies of *H*_{t} and *H*_{b}, respectively.

### D. Past algorithms

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 *H*_{t} in order to compute exp(−*βH*_{t}). However, this quickly becomes intractable unless there are many symmetries present in the system. Moreover, even if exp(−*βH*_{t}) 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(−*βH*_{t}) is a matrix with 2^{20} × 2^{20} 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(−*βH*_{t}) and instead compute quantities like ⟨*v*| exp(−*βH*_{t})|*v*⟩. Moreover, because Krylov subspace methods are “matrix-free,” they only access *H*_{t} 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 *H*_{s}. At low temperature, *ρ**(*β*) depends only on the ground state(s) of *H*_{t} and *H**(*β*) depends only on the ground state energies of *H*_{t} and *H*_{b}.

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(−*βH*_{t}), without ever explicitly constructing the matrix, has not been thoroughly studied.

## III. ALGORITHM

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.

### A. Block stochastic Lanczos quadrature

Let $A\u2208LHt$ and recall that (*I*_{s} ⊗ ⟨*v*|)*A*(*I*_{s} ⊗ |*v*⟩) is an unbiased estimator for $|Hb|\u22121trb(A)$. To improve the estimator’s accuracy, we can average multiple copies. Specifically, given *n*_{v} iid samples {|*v*_{j}⟩} of |*v*⟩ (with HA[|*v*⟩⟨*v*|] = *I*_{b}), we can define the averaged estimator

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 physics^{15,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 *n*_{v} is large. We believe the mathematical correspondence between $|Hb|$ and *n*_{v} is worthy of a more rigorous mathematical analysis. For instance, generalized fundamental thermodynamic relations were recently unraveled by replacing thermodynamic infinite-size limit $(|Hb|\u2192\u221e)$ with multiple-measurement limit (*n*_{v} → ∞).^{33}

Often, *A* = *f*[*H*] for some function $f:R\u2192R$ and $H\u2208LHt$; for instance, *f* = *x* ↦ exp(−*βx*) and *H* = *H*_{t}. If *H* has known diagonalization, then we can easily compute (*I*_{s} ⊗ ⟨*v*|)*f*[*H*](*I*_{s} ⊗ |*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 *I*_{s} ⊗ |*v*⟩ for *k* iterations to obtain a $k|Hs|\xd7k|Hs|$ block-tridiagonal matrix *T*. Given a matrix *V* with *V*^{†}*V* = *I*_{s}, Algorithm 1 computes orthonormal matrices {*Q*_{j}}, *j* = 1, …, *k* + 1 such that $Qi\u2020Qj=\delta i,jIs$ and for all *j* ≤ *k*,

These vectors satisfy a symmetric block-tridiagonal recurrence

where *E*_{k} = |*k*⟩ ⊗ *I*_{s} and

1: procedure BLOCK-LANCZOS(H, V, k) |

2: Q_{1} = V, |

3: for j = 1, 2, …, k do |

4: $Z=HQj\u2212Qj\u22121Bj\u22121\u2020$ |

5: $Aj=Qj\u2020Z$ |

6: Z = Z − Q_{j}A_{j} |

7: Q_{j+1}, B_{j} = qr(Z) |

8: return {Q_{j}}, {A_{j}}, {B_{j}} |

1: procedure BLOCK-LANCZOS(H, V, k) |

2: Q_{1} = V, |

3: for j = 1, 2, …, k do |

4: $Z=HQj\u2212Qj\u22121Bj\u22121\u2020$ |

5: $Aj=Qj\u2020Z$ |

6: Z = Z − Q_{j}A_{j} |

7: Q_{j+1}, B_{j} = qr(Z) |

8: return {Q_{j}}, {A_{j}}, {B_{j}} |

The expression (*I*_{s} ⊗ ⟨*v*|)*f*[*H*](*I*_{s} ⊗ |*v*⟩) can be approximated by

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 2*k* − 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, …, *n*_{v}, denoting by *T*_{j} the resulting block-tridiagonal matrix obtained by block Lanczos run on *H* and *I*_{s} ⊗ |*v*_{j}⟩, our estimator for tr_{b}(*f*[*H*]) is

#### 1. Costs

The accuracy and computational cost of (34) depend on both *n*_{v} 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(|Ht\Vert Hs|)=O(|Hb\Vert Hs|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 *n*_{v} while the storage cost is constant. On the other hand, all terms can be computed entirely in parallel at the cost of *n*_{v} times more storage.

Like other typicality algorithms, our approach requires matrix-products with the total system Hamiltonian *H*_{t} 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.

### B. An algorithm for the mean force Hamiltonian

To compute *H**(*β*), we can use the same test states {|*v*_{j}⟩} 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 *H*_{t} and (*I*_{s} ⊗ |*v*_{j}⟩) and $(Tb)j$ is the tridiagonal matrix produced by Lanczos with *H*_{b} and |*v*_{j}⟩, by (7) we obtain the estimators

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(−*β*(*x* − *E*_{0})), where *E*_{0} 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 *H*_{t}.

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

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

## IV. ERROR ANALYSIS

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

### A. Trace estimators

Tail bounds trace estimators were studied in Refs. 35 and 36, with more recent and refined analyses given in Refs. 37–39; see Ref. 29 for historical context. For constants *ɛ* > 0 and *δ* ∈ (0, 1), and $B\u2208L(Hb)$, these analyses aim to bound the number of samples *n*_{v} required so that

The resulting bounds are typically simple functions depending on the distribution of |*v*⟩, the value of *n*_{v}, 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(\Vert B\Vert 22\Vert Hb|\u22121\epsilon 2)ln(2/\delta )$. Here, $O\u0303$ is equivalent to $O$ but with poly-logarithmic factors in the constituent parameters suppressed for readability; i.e., we say, a variable is $O\u0303(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 $A\u2208LHt$ as

where $Am,n\u2208L(Hb)$. Fix *ɛ* > 0 and *δ* ∈ (0, 1), and suppose that for all *m*, *n*, we have chosen *n*_{v} so that

for some $\epsilon \u0303$ and $\delta \u0303$, the exact values of which will soon become apparent. Applying a union bound over all pairs *m*, *n*, we obtain the bound

Next, note that

and that

For any $X\u2208L(Hs)$, we have that $\Vert X\Vert \u2264|Hs|maxm,n\u27e8m|X|n\u27e9$. Putting everything together, we find that

if we take $\epsilon \u0303=\epsilon /|Hs|$ and $\delta \u0303=\delta /|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 *n*_{v} iid samples, *n*_{v} must scale like $O(\epsilon 2)$.

#### 1. A note on high temperatures

At high temperature, our estimator will compute tr_{b}(exp(−*βH*_{t})) 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*(\beta )=Hs+trb(Hsb)/tr(Ib)+O(\beta )$. Thus, at high temperature, the approximation (36) differs from *H*_{s} by an additive factor $nv\u22121\u2211j=1nv(\u27e8vj|\u2297Is)Hsb(|vj\u27e9\u2297Is)/tr(Ib)$, which should be corrected for. Since HA[(⟨*v*_{j}| ⊗ *I*_{s})*H*_{sb}(|*v*_{j}⟩ ⊗ *I*_{s})] = 0_{s}, 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(−*βH*_{t}) 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 tr_{b}(exp(−*βH*_{t})), many samples are required. In this setting, other approaches such as low-rank approximation or hybrid methods make more sense.^{38,40–43}

### B. Block Gauss quadrature

Recall that $V\u2020p(A)V=E1\u2020p[T]E1$ for all polynomials *p* of degree at most 2*k* − 1. For brevity, let $Error=\Vert V\u2020f[H]V\u2212E1\u2020f[T]E1\Vert $. Then, for any polynomial *p* with deg(*p*) ≤ 2*k* − 1,

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

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 *T*_{j} 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 arithmetic^{46,47} can be leveraged to show that the Lanczos algorithm still works well for the task of applying matrix functions to vectors^{48} and quadratic forms^{49} 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 method^{23}); 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 $E1\u2020f[T]E1$ computed by the block Lanczos algorithm still provides a good approximation to *V*^{†}*f*[*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.

## V. NUMERICAL EXPERIMENTS

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$.

### A. Solvable system

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 (*N*_{s} = 2), and put the magnetic field strength at *h* = 0.3*J*. We then run our algorithm using *k* = 30 Lanczos iterations and *n*_{v} = 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.

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 *n*_{v}. 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 *H*_{s}, while at low temperature, the spectrum converges to a constant *E*_{t} − *E*_{b}.

### B. Varying choice of system spins

We now consider a spin ladder. The edge coupling strength is set to *J*, the run coupling strength is set to −0.45*J*, and the magnetic field strength is set to *J*. We set *N* = 20 and consider systems of size *N*_{s} = 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 *n*_{v} = 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.

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 (*H*_{sb} = 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.

### C. Long range interactions

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|i\u2212j|\u2212\alpha $, *α* ≥ 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 0*J* to 2.2*J*. We run our algorithm with *n*_{v} = 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.

At zero temperature, the von Neumann entropy indicates the presence of a magnetic field-induced quantum phase transition^{50–53} from an entangled state to a pure state of the system of interest around *h* ≈ 1.4*J*. In the exactly solvable chain with *α* = ∞, a similar quantum phase transition from entangled to pure state occurs at *h* = 2*J*.^{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 *T*_{c} ≈ 0.5*J*/*β*, 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 *H*_{sb} 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(*H*_{s}*ρ*_{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

can be interpreted as the “deviation” in the system’s internal energy (state function) by coupling to the bath via *H*_{sb}; 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 *n*_{v} = 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.

Borrowing the idea from solvation thermodynamics^{59} 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.

### D. Strong to weak coupling

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

where *ɛ* ≥ 0 determines the coupling strength and $H\u0304s$, $H\u0304b$, and *H*_{sb} are all as in the previous example. We again run our algorithm with *n*_{v} = 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.

## VI. CONCLUSION

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 estimation^{38,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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

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

### APPENDIX A: DERIVATION OF HIGH AND LOW TEMPERATURE LIMITS

#### 1. High temperature limit

Expanding at *β* = 0, we find

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

Thus, combining these expressions,

where, in the final equality, we have used that

Now, recall that *H*_{sb} accounts for interactions between the system and bath and is therefore a linear combination of terms of the form

where $i+i\u2032=(2s+1)Ns\u22121$ and $j+j\u2032=(2s+1)Nb\u22121$. Applying basic properties of the partial trace, we see that

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 tr_{b}(*H*_{sb}) = 0_{s}.

We have, therefore, established that

#### 2. Low temperature limit

Write the eigenvalue decomposition of *H*_{t} as $Ht=\u2211i=1|Ht|Ei|\psi i\u3009\u3008\psi i|$ for orthogonal eigenvectors |*ψ*_{i}⟩. Then,

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

Note that

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

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

We, therefore, have that

where *E*_{t} and *E*_{b} are the ground state energies of *H*_{t} and *H*_{b}, respectively.

### APPENDIX B: A CONSISTENT PARTIAL TRACE ESTIMATOR

Consider a random pure product state |*v*_{s}⟩ ⊗ |*v*_{b}⟩ where HA[|*v*_{s/b}⟩⟨*v*_{s/b}|] = *I*_{s/b}. We have that

Thus, estimators of this form have the desirable property that

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 |*v*_{s}⟩ ⊗ |*v*_{b}⟩ 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*⟩ = |*v*_{1}⟩ ⊗⋯⊗ |*v*_{N}⟩ 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}

### APPENDIX C: SOLUTION TO THE SPIN CHAIN

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

Next, define

For *N*_{s} = 2, the eigenvalues of *ρ**(*β*) are given by

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

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

### APPENDIX D: PHASE PLOT FOR SOLVABLE SPIN CHAIN

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.

## REFERENCES

While *H*_{t} may be sparse, the matrix exponential is not, in general, sparse.

_{3}Ru

_{2}O

_{7}