In classical statistical mechanics, thermal expectation values are computed by averaging over some probability distribution. In the case of the canonical ensemble, this looks like a sum over all states \(s\) weighted by the exponential of the energy \(E_s\): $$ \langle\mathcal O\rangle_{\text{classical}} = Z^{-1} \sum_s e^{-E_s / T} \mathcal O_s \text. $$ The normalization of \(Z^{-1}\) is irrelevant to the discussion here — in fact, I’ll leave it off all future expressions.

How does this generalize to quantum mechanics? The usual way is to replace the energies with a Hermitian operator (the Hamiltonian), the observable with another Hermitian operator, and the sum with a trace. An expectation value in the quantum canonical ensemble looks like $$ \langle \mathcal O\rangle_{\text{quantum}} = \mathrm{Tr}\;e^{-H / T} \mathcal O $$ To see that this is in fact a generalization of the classical ensemble above, consider the case where the Hamiltonian and observable commute. A basis exists where they are both diagonal, and in that basis the quantum expectation value takes exactly the form of the classical expectation value above.

So far, so familiar. But is this the only sensible generalization of classical statistical mechanics? In a sense, yes. First let’s go back to the classical case. The exponential of the energy \(e^{-E/T}\) is referred to as the Boltzmann factor. If we know the expectation values of a bunch of operators \(\mathcal O_1,\ldots,\mathcal O_n\), then as long as these expectation values are physically consistent, a suitable Boltzmann factor can always be constructed to satisfy them. (Actually, that’s the definition of physically consistent.) Furthermore, if the system is finite and we know the expectation values of enough operators, the Boltzmann factor is actually uniquely determined! (Well, up to a measure-zero set, or something like that.)

That should demotivate us from looking for new formulations of classical statistical mechanics. Similar logic holds for the quantum case. In a quantum system, we’re also constrained by linearity, but the structure of the argument is the same. Given enough (physically realizable) expectation values, the density matrix \(e^{-H / T}\) is uniquely determined. As a result, there’s not much point looking for different formulations of quantum statistical mechanics, either.

Let’s do it anyway.

In a system with \(Q\) qubits, the density matrix is a bulky object with \(2^{2Q}\) entries. I don’t want to deal with that, I want to work with states — nice sleek vectors with a measly \(2^{Q}\) entries. The classical ensembles give probability distributions on classical states; why can’t a quantum ensemble give a probability distribution on quantum (pure) states?

Well, in fact, the density matrix already does. By diagonalizing \(H\), you can see that \(e^{-H / T}\) naturally yields a probability distribution on the set of eigenstates of the Hamiltonian. This isn’t very nice, though, because it sort of presupposes that we know the eigenbasis of the Hamiltonian. For any interesting problem, I don’t.

Fortunately, it’s pretty clear that this probability distribution is far from unique. First, note that every probability distribution on states will induce some density matrix. This follows from the discussion above; alternatively, the density matrix is given explicitly by $$ \rho[p] = \int d |\psi\rangle\; p(\psi)\; |\psi\rangle\langle\psi| \text. $$ Now count dimensions. The space of density matrices is \(2^{2Q}\) dimensional, give or take. The space of possible probability distributions on the continuous space of states? Ah… many more dimensions there.

So, many more probability distributions compatible with the thermal density matrix should exist. As far as I know, finding natural such probability distributions is largely unexplored territory. (There’s one induced by ergodicity and the time-evolution operator \(e^{-i H t}\), although I’m not even sure if it’s unique.) Nevertheless, one such (approximate) construction is given by Sugiura and Shimizu. Start with a uniform (i.e., \(SU(2^Q)\)-symmetric) distribution on the space of states, corresponding to an infinite temperature ensemble. Now take one state, and apply the operator \(e^{-H / (2T)}\). Lastly, if you think of states as being vectors in a Hilbert space, then normalize the result; if you think of states as rays, this is of course unnecessary.

In the thermodynamic limit, expectation values with respect to this distribution match the canonical quantum ensemble described above. In fact, in that limit, a single quantum state sampled from this ensemble yields (almost certainly) all the correct expectation values. A careful proof is given in the paper; a nice physical argument follows just as it would for a classical system. In the infinite-volume limit, we can consider multiple sub-systems each individually in the thermodynamic limit and all uncoupled from each other. Therefore, what appears to be a single thermodynamic “sample” actually contains infinitely many, equally thermodynamic samples.