# Neural Network Thermodynamics

Dennis Wagner  Department of Physics and Astronomy, University of Manitoba,
Winnipeg R3T 2N2, Canada  University Wuppertal, Wuppertal, Germany  Andreas
Klümper  University Wuppertal, Wuppertal, Germany  Jesko Sirker
<EMAIL_ADDRESS>Department of Physics and Astronomy, University of
Manitoba, Winnipeg R3T 2N2, Canada  Manitoba Quantum Institute, University of
Manitoba, Winnipeg R3T 2N2, Canada

###### Abstract

We present three different neural network algorithms to calculate
thermodynamic properties as well as dynamic correlation functions at finite
temperatures for quantum lattice models. The first method is based on
purification, which allows for the exact calculation of the operator trace.
The second one is based on a sampling of the trace using minimally entangled
states, whereas the third one makes use of quantum typicality. In the latter
case, we approximate a typical infinite-temperature state by wave functions
which are given by a product of a projected pair and a neural network part and
evolve this typical state in imaginary time.

##  I Introduction

While being able to obtain precise results for the ground states of quantum
lattice Hamiltonians is essential to understand different states of matter and
quantum phase transitions between such states, phase transitions can also
occur at finite temperatures $T$ for dimensions $d\geq 2$. Furthermore,
experiments are always performed at $T>0$. The influence of thermal
fluctuations on the physical properties of a quantum system are therefore also
of fundamental interest. Experiments on solid state systems, for example,
often probe thermodynamic quantities such as the specific heat, the
susceptibility, or the compressibility. Also of interest are linear response
functions at finite temperatures such as the spin structure factor, measured
in neutron scattering experiments, or the spin-lattice relaxation rate,
measured in nuclear magnetic resonance.
Except for exactly solvable quantum models such as Gaussian fermionic and
bosonic models and Bethe ansatz integrable models in one dimension—which share
the common property of having an infinite set of local conservation
laws—approximate analytical or numerical methods are needed to study the
thermodynamic properties of a quantum system. For a closed system, we are
usually interested in calculating the partition function $Z$ which gives
access to the free energy $f$ as well as in expectation values of operators
$A$. These quantities are obtained in the canonical ensemble as

$Z=\text{tr}\,\text{e}^{-\beta H},\;f=-\frac{T}{N}\ln Z,\;\langle
A\rangle=\frac{\text{tr}\\{A\text{e}^{-\beta H}\\}}{Z}\ ,$ (1)

where $T=1/\beta$ denotes the temperature, $H$ is the Hamiltonian, $N$ the
number of lattice sites, and we set $k_{B}=1$. Important numerical methods to
obtain approximations of the partition function and expectation values
include, in particular, various Quantum Monte Carlo (QMC) techniques [1, 2].
While QMC can be applied to problems in any dimension, it suffers from a sign
problem in the sampling weights for fermionic models and frustrated spin
models [3, 4]. In such cases, it is usually not possible to reach low
temperatures. Furthermore, dynamical quantities at finite temperatures can
only be obtained in imaginary time (Matsubara frequency) with the analytical
continuation of the numerical data to real times being an ill-posed problem
with no unique solution [5]. An alternative approach to calculate the
thermodynamic properties of quantum systems are matrix product state/density-
matrix renormalization group (MPS/DMRG) algorithms [6, 7]. This approach was
originally invented to obtain the ground state properties of one-dimensional
quantum systems [6]. Two different ideas were then pursued to extend such
algorithms to finite temperatures: On the one hand, one can map any one-
dimensional quantum lattice model at finite temperatures to a two-dimensional
classical model. For the classical model, one can then define a so-called
quantum transfer matrix (QTM) which evolves along the spatial direction and
whose largest eigenvalue and corresponding eigenstate completely determine the
properties of the system in the thermodynamic limit. The QTM can then be
approximated as a matrix product operator [8, 9, 10, 11]. In this approach, it
is also possible to calculate dynamical correlation functions directly in real
time [12, 13].
An alternative approach is based on introducing an ancilla site for each real
site of the lattice and to obtain the partition function by purification [14].
Geometrically, this corresponds to replacing the cylinder which represents the
partition function in the two-dimensional classical model by two sheets
(formed by the real sites and ancilla sites) connected at the two edges by a
state which is a direct product of maximally entangled local states between a
site and its corresponding ancilla. Similar to the QTM algorithm, one can also
include a real-time evolution and thus access dynamical correlations at finite
temperatures [12, 15, 16, 17, 18]. While fermionic and frustrated spin models
in one dimension do not pose any difficulties for MPS-type algorithms in
contrast to QMC, they can only be applied to small clusters in dimensions
higher than one in which case calculations are usually restricted to the
ground state [19, 20].
Recently, variational QMC approaches, where the parameters of the variational
wave function are optimized using a neural network, have attracted renewed
attention [21, 22, 23, 24] as a numerical approach possibly able to overcome
some of these issues. Variational QMC approaches have, at least in principle,
the advantage that they are applicable in any dimension, that sign problems
can be circumvented, and that real-time dynamics can be studied directly.
However, the choice of the variational wave function biases the obtained
results. It is expected that by using ansatz wave functions based on neural
networks with many variational parameters, the system is able to learn by
itself which parameters are most important instead of fixing the form of the
wave function by physical intuition. Mathematically, the basis for this
approach are universal approximation theorems which guarantee that any smooth
real function can be approximated to any desired precision by a neural network
with a non-polynomial activation function (for example, the sigmoid function)
[25, 26]. Such approximation theorems can also be generalized to complex-
valued neural networks [27] which are needed to approximate complex wave
functions. In practice, some form of gradient descent in the parameter space
of the variational wave function is then used to optimize the ground state
energy. The variational principle guarantees that the obtained energy is
always larger than or equal to the true ground state energy. This approach has
been applied successfully to investigate the ground state properties of a
number of different low-dimensional quantum lattice models. However, for
fermionic and frustrated spin systems where a Marshall sign rule [28] does not
exist, learning the sign structure has remained challenging at least in ansatz
wave functions which are explicitly translationally invariant [29, 30, 31].
In this paper, we want to expand such neural network approaches to allow for
the calculation of static and dynamic finite-temperature properties. We will
discuss three possible approaches: Similar to MPS/DMRG algorithms, we can
introduce an additional ancilla site for each real lattice site and obtain
traces via purification. The second approach is based on sampling the trace
using minimally entangled states [32]. The third approach is based on the
notion of quantum typicality. Properly implemented, a single pure state can
reproduce the expectation values of observables in the Gibbs ensemble for
large system sizes [33, 34, 35]. We note that the second approach has also
very recently been studied in the neural network context in Ref. [36]. In all
cases, we start out by constructing an appropriate infinite-temperature state
$|\Psi(\infty)\rangle$ and then adjust the parameters of the neural network by
solving the Schrödinger equation for the imaginary time evolution

$|\Psi(\beta)\rangle=\exp(-\beta H/2)|\Psi(\infty)\rangle\,.$ (2)

We note that the task here is technically more demanding than in the ground
state case. While any gradient descent through parameter space which is
successively lowering the energy of the system is acceptable if we are only
interested in obtaining an approximation for the ground state, we want to
follow here as precisely as possible the thermodynamic path so that the state
does describe the physical properties of the system at any given temperature
$T=1/\beta$.
Our paper is organized as follows: In Sec. II, we present the neural network
algorithm based on purification and test the algorithm by calculating
thermodynamic and dynamic properties at finite temperatures for the Heisenberg
model. In Sec .III, we instead sample the trace using minimally entangled
states. Thermodynamic properties are computed by evolving these states in
imaginary time. In Sec. IV, we then introduce, as a third alternative, an
algorithm based on quantum typicality and compare the obtained results with
those obtained by purification and the sampling of the trace. The last section
is devoted to a short summary and a critical discussion of the power of neural
network algorithms for thermodynamic simulations and of some of the remaining
open issues.

##  II Purification

In order to exactly calculate the partition function $Z$ as well as
thermodynamic expectation values, a trace over all basis states has to be
performed. Since we do not want to diagonalize the Hamiltonian itself but
rather express wave functions in the Hilbert space by a variational ansatz, we
cannot perform the trace directly. Let us assume that $|\\{s\\}\rangle$ is a
local basis for the quantum lattice model we want to consider with dimension
$d$. Then, one way to obtain thermodynamic quantities solely based on wave
functions and the action of operators on these wave functions is to introduce
additional ancilla degrees of freedom ${s^{\prime}}$ by

$|\\{s\\}\rangle\to|\\{s\\}\rangle\otimes|\\{s^{\prime}\\}\rangle\equiv|\\{s,s^{\prime}\\}\rangle\
.$ (3)

While the Hamiltonian will act only on the real lattice sites, the state of
the ancillas remains unchanged. We can therefore think of the ancillas as an
‘index’ keeping track of the initial state of the real sites before we have
acted on this state with some operator. In particular, we will consider the
initial state

$|\Psi(\infty)\rangle=\bigotimes_{i=1}^{N}\sum_{s=-s^{\prime}}|s,s^{\prime}\rangle_{i}\
,$ (4)

where $N$ is the number of lattice sites. $|\Psi(\infty)\rangle$ represents
the (not normalized) infinite-temperature state and an expectation value or a
dynamical correlation function can be obtained by

$\displaystyle\langle\hat{O}_{1}(0)\hat{O}_{2}(t)\rangle_{T=\infty}$
$\displaystyle=$
$\displaystyle\frac{\langle\Psi(\infty)|(\hat{O}_{1}(0)\otimes\mathbb{1})(\hat{O}_{2}(t)\otimes\mathbb{1})|\Psi(\infty)\rangle}{\langle\Psi(\infty)|\Psi(\infty)\rangle}$

with $\langle\Psi(\infty)|\Psi(\infty)\rangle=d^{N}$. A finite temperature
state is then obtained by imaginary time evolution

$|\Psi(\beta)\rangle=(\exp(-\beta H/2)\otimes\mathbb{1})|\Psi(\infty)\rangle$
(6)

and the partition function is given by
$Z(\beta)=\langle\Psi(\beta)|\Psi(\beta)\rangle$. Note that the two
sheets—consisting of the real spins and the ancilla spins, respectively—become
connected when performing a partial trace over the ancillas and thus together
form a cylindrical geometry. In the numerical calculations, we therefore find
it advantageous to move half of the imaginary time evolution to the ancillary
spins and to apply $\exp(-\beta H/4)$ on both sheets. This makes the problem
more symmetric which is helpful for an efficient variational calculation.
Similarly, the energy per site is then obtained in a more symmetric manner as

$e=\frac{1}{2N}\frac{\langle\Psi(\beta)|(H\otimes\mathbb{1}+\mathbb{1}\otimes
H)|\Psi(\beta)\rangle}{\langle\Psi(\beta)|\Psi(\beta)\rangle}\,.$ (7)

The description of the purification algorithm so far is general and can be
performed with any variational ansatz wave function, provided that the initial
infinite-temperature state (4) can be represented accurately. As a proof of
principle, we consider in the following a simple implementation using a
modified restricted Boltzmann machine (mRBM). We note that the performance of
the algorithm discussed below might possibly be improved by using deep neural
networks (DNN) as, for example, convolutional neural networks (CNN), along
similar lines as in ground state algorithms [24]. Here the focus is on
demonstrating how thermodynamic calculations can be performed in principle and
what the fundamental limitations are. For convenience, we label both the real
and the ancilla spins as $s_{i}$ with $i=1,\cdots,N$ for the real spins and
$i=N+1,\cdots,2N$ for the ancillas. For the components of the wave function,
we then choose the following ansatz

$\displaystyle\Psi(s)=\langle s|\Psi\rangle$ (8) $\displaystyle=$
$\displaystyle\exp\left(-\sum_{i=1}^{2N}a_{i}s_{i}\right)\prod_{j=1}^{M}\cosh\left(b_{j}+\sum_{i=1}^{2N}W_{ij}s_{i}\right)$
$\displaystyle\times$
$\displaystyle\exp\left(-\sum_{i=1}^{2N}c_{i}(s_{i}+s_{(i+N)\textrm{mod}(2N)})^{2}\right)\
.$

Here $\boldsymbol{a},\boldsymbol{b}$ and $\boldsymbol{c}$ are complex vectors
with $2N$ and $M$ components, respectively, while $\mathbf{W}$ is a complex
matrix of dimension $2N\times M$. The second line in Eq. (8) is the standard
RBM ansatz for a system with $2N$ sites [21]. We modify this ansatz by
multiplying with the term in the last line, which couples every spin with its
respective ancilla. The (not normalized) infinite-temperature state
$|\Psi(\infty)\rangle$ can now formally be obtained by setting
$\boldsymbol{a},\boldsymbol{b},\boldsymbol{W}=0$ and $c_{i}\to+\infty$,
$i=1,\cdots,2N$. This way, the term in the last line exponentially suppresses
all components where real site and ancilla site do not have opposite signs and
gives all components where the signs are opposite equal weight in the wave
function. In practice, however, this limit cannot be performed exactly because
numerical gradient methods—required to change the parameters of the wave
function during the imaginary time evolution—will then not be able to find the
proper descent path. We therefore set $c_{i}=c=\mbox{const}>0$ and choose the
real and imaginary parts of $\boldsymbol{a}$ and $\boldsymbol{W}$ from a
random normal distribution centered around zero with a very small width $w\sim
10^{-2}$ and keep $\boldsymbol{b}=0$.
In order to time evolve this approximation of the infinite-temperature state,
we use the time-dependent variational principle (TDVP) which is a simple
consequence of Eq. (6). Taking into account that the temperature is encoded in
the variational parameters
$\mathbb{W}=\mathbb{W}[\beta]\equiv\\{\boldsymbol{c}[\beta],\boldsymbol{a}[\beta],\boldsymbol{b}[\beta],\boldsymbol{W}[\beta]\\}$
we obtain

$\frac{d}{d\beta}|\Psi(\mathbb{W}[\beta])\rangle=\sum_{\alpha}\frac{\partial\mathbb{W}_{\alpha}}{\partial\beta}|\partial_{\mathbb{W}_{\alpha}}\Psi(\mathbb{W}[\beta])\rangle=-\frac{H}{2}|\Psi(\mathbb{W}[\beta])\rangle\,.$
(9)

We can now multiply with the vectors
$\langle\partial\mathbb{W}_{\alpha^{\prime}}\Psi|$ from the left to obtain

$\sum_{\alpha}\underbrace{\frac{\langle\partial_{\mathbb{W}_{\alpha^{\prime}}}\Psi|\partial_{\mathbb{W}_{\alpha}}\Psi\rangle}{\langle\Psi|\Psi\rangle}}_{\boldsymbol{G}}\frac{\partial\mathbb{W}_{\alpha}}{\partial\beta}=-\frac{1}{2}\underbrace{\frac{\langle\partial_{\mathbb{W}_{\alpha^{\prime}}}\Psi|H|\Psi\rangle}{\langle\Psi|\Psi\rangle}}_{\boldsymbol{H}}\
,$ (10)

where we have defined the matrix $\boldsymbol{G}$ and the vector
$\boldsymbol{H}$. Since the parameters $\mathbb{W}_{\alpha}$ are not
necessarily linearly independent, the matrix $\boldsymbol{G}$ is, in general,
not regular. We can, however, always compute the Moore-Penrose pseudo-inverse
with a small cutoff which we denote by $\boldsymbol{G}^{-1}$. This leads to
the following matrix equation describing the time evolution of the state
$|\Psi(\mathbb{W}[\beta]\rangle$ via a change of the variational parameters
$\mathbb{W}$

$\frac{\partial\mathbb{W}}{\partial\beta}=-\frac{1}{2}\boldsymbol{G}^{-1}\boldsymbol{H}\Rightarrow\mathbb{W}[\beta+\Delta\beta]\approx\mathbb{W}[\beta]-\frac{1}{2}\Delta\beta\boldsymbol{G}^{-1}\boldsymbol{H},$
(11)

where $\Delta\beta$ is the imaginary time step (learning rate) and on the
r.h.s. of Eq. (11) we have written, for simplicity, the approximate solution
of the differential equation in the form of a Newton forward method. In the
numerical calculations presented below, we will use instead the Heun method
which is more stable, in particular, if the states are Monte Carlo sampled.
The elements of the matrix $\boldsymbol{G}$ and the vector $\boldsymbol{H}$
can be expressed in terms of the probability distribution $|\Psi(s)|^{2}$, for
example,

$\frac{\langle\partial_{\mathbb{W}_{\alpha^{\prime}}}\Psi|\partial_{\mathbb{W}_{\alpha}}\Psi\rangle}{\langle\Psi|\Psi\rangle}=\frac{\sum_{s}|\Psi(s)|^{2}\left(\frac{\partial_{\mathbb{W}_{\alpha^{\prime}}}\Psi(s)}{\Psi(s)}\right)^{*}\frac{\partial_{\mathbb{W}_{\alpha}}\Psi(s)}{\Psi(s)}}{\sum_{s}|\Psi(s)|^{2}}\
.$ (12)

###  II.1 Adaptive Heun Method

The adaptive Heun method is used throughout this paper for the imaginary time
evolution described in Eqs. (6) and (11). In the following, we will use $\tau$
to denote the imaginary time instead of $\beta$. The reason for this change in
notation will become clear in Sec. II.2. The adaptive Heun method is an
iterative method where in every time step the error of the integration is
estimated and the length of the time step $\Delta\tau$ adjusted in a way that
the error stays constant. We describe this predictor-corrector method in the
following for a second order integration for which the standard notation for
the derivative $\dot{\mathbb{W}}=f(\mathbb{W})$ and the exact solution
$\mathbb{W}(\tau+\Delta\tau)$ at time $\tau+\Delta\tau$, where $\mathbb{W}$ is
a $P$-component vector and $P$ the number of parameters, is used [24]. We
define

$\displaystyle k_{1}$ $\displaystyle=f(\mathbb{W}_{n})$ $\displaystyle k_{2}$
$\displaystyle=f(\mathbb{W}_{n}+\Delta\tau\cdot k_{1})$
$\displaystyle\mathbb{W}_{n+1}$
$\displaystyle=\mathbb{W}_{n}+\frac{\Delta\tau}{2}\left(k_{1}+k_{2}\right),$
(13)

which leads to a cubic error in the integration step

$\mathbb{W}_{n+1}=\mathbb{W}(\tau+\Delta\tau)+\alpha(\Delta\tau)^{3},$ (14)

where $\alpha$ is a constant vector. This is now compared to the error for two
time steps of length $\frac{\Delta\tau}{2}$

$\mathbb{W}^{\prime}_{n+1}=\mathbb{W}(\tau+\Delta\tau)+\underbrace{2\alpha\left(\frac{\Delta\tau}{2}\right)^{3}}_{\delta},$
(15)

where $\delta$ is the integration error. The difference of these two results

$||\mathbb{W}_{n+1}-\mathbb{W}^{\prime}_{n+1}||=\frac{3}{4}||\alpha||(\Delta\tau)^{3}\equiv
3\delta$ (16)

allows to adjust the time step

$\Delta\tau^{\prime}=\Delta\tau\left(\frac{\varepsilon}{\delta}\right)^{\frac{1}{3}}$
(17)

to obtain a desired accuracy $\varepsilon$. The norm in Eq. (16) is defined by
$||x||=\frac{1}{P}\sqrt{x^{\dagger}\mathbf{G}x}$. Note that the tolerance
$\varepsilon$ and the initial step size $\Delta\tau$ are in principle
parameters which need to be specified at the beginning of the calculation. The
adjusted step size $\Delta\tau^{\prime}$ is always used as the proposed step
size $\Delta\tau$ in the next step. By applying a pretempering at $\tau=0$ of
a few steps in which parameters are not changed, $\Delta\tau^{\prime}$
converges to a value that is used as the initial step size. Thus, the only
remaining parameter to be fixed is the desired tolerance $\varepsilon$ which
defines the accuracy of the method. In the following, we will denote the
number of iterative steps in the Heun algorithm by $R$.

###  II.2 Beta Modification

The variational imaginary-time evolution in Eq. (9) needs a more thorough
consideration that takes into account the projection onto the variational
manifold, which is a submanifold of the projective Hilbert space
$\mathcal{P}(\mathcal{H})$ representing all physical states. We start from the
projected Schrödinger equation by using the projector
$\mathtt{P}_{\Psi_{\mathbb{W}}(\tau)}$ which restricts the change of the
variational state to be in the tangent space
$\mathcal{T}_{\Psi_{\mathbb{W}}(\tau)}\mathcal{M}$ of the variational manifold
$\mathcal{M}$ [37, 36]:

$\displaystyle\mathtt{P}_{\Psi_{\mathbb{W}}(\tau)}\frac{d|\Psi_{\mathbb{W}}(\tau)\rangle}{d\tau}=$
$\displaystyle-\frac{1}{2}\mathtt{P}_{\Psi_{\mathbb{W}}(\tau)}H|\Psi_{\mathbb{W}}(\tau)\rangle$
(18) $\displaystyle=$
$\displaystyle-\frac{1}{2}\mathtt{P}_{\Psi_{\mathbb{W}}(\tau)}(H-E(\tau))|\Psi_{\mathbb{W}}(\tau)\rangle$

This is accomplished by using the orthogonal projector

$\mathtt{P}_{\Psi_{\mathbb{W}}(\tau)}=\frac{2(\tilde{\boldsymbol{G}}^{-1})^{{}^{\mu\eta}}}{\langle\Psi_{\mathbb{W}}(\tau)|\Psi_{\mathbb{W}}(\tau)\rangle}|V_{\mathbb{W},\mu}(\tau)\rangle\langle
V_{\mathbb{W},\eta}(\tau)|\ ,$ (19)

where
$|V_{\mathbb{W},\mu}(\tau)\rangle=\mathtt{Q}_{\Psi_{\mathbb{W}}(\tau)}\partial_{\mathbb{W}_{\mu}}|\Psi_{\mathbb{W}}(\tau)\rangle$
is a local basis of the tangent space
$\mathcal{T}_{\Psi_{\mathbb{W}}(\tau)}\mathcal{M}$ at the state
$|\Psi_{\mathbb{W}}(\tau)\rangle$ and

$\tilde{\boldsymbol{G}}=\frac{\langle
V_{\mathbb{W},\mu}(\tau)|V_{\mathbb{W},\eta}(\tau)\rangle}{\langle\Psi_{\mathbb{W}}(\tau)|\Psi_{\mathbb{W}}(\tau)\rangle}\
.$ (20)

The projector

$\mathtt{Q}_{\Psi_{\mathbb{W}}(\tau)}=\left(\mathds{1}-\frac{|\Psi_{\mathbb{W}}(\tau)\rangle\langle\Psi_{\mathbb{W}}(\tau)|}{\langle\Psi_{\mathbb{W}}(\tau)|\Psi_{\mathbb{W}}(\tau)\rangle}\right)$
(21)

allows one to define unique representatives of the tangent space. Due to the
projection, the time traveled in the manifold and the imaginary time are not
equal. This is quantified by the following ratio [36]

$\displaystyle\begin{split}\gamma=&\frac{||\mathtt{P}_{\Psi_{\mathbb{W}}(\tau)}(H-E(\tau))|\Psi_{\mathbb{W}}(\tau)\rangle||^{2}}{||(H-E(\tau))|\Psi_{\mathbb{W}}(\tau)\rangle||^{2}}\\\
=&-\frac{1}{\sigma^{2}(\tau)}\mathtt{P}_{\Psi_{\mathbb{W}}(\tau)}\left(\frac{dE(\tau)}{d\tau}\right)\
,\end{split}$ (22)

where the two quantities of interest on the right hand side are easily
calculated at every time step

$\displaystyle E(\tau)=$
$\displaystyle\frac{\langle\Psi_{\mathbb{W}}(\tau)|H|\Psi_{\mathbb{W}}(\tau)\rangle}{\langle\Psi_{\mathbb{W}}(\tau)|\Psi_{\mathbb{W}}(\tau)\rangle}$
(23) $\displaystyle\sigma^{2}(\tau)=$
$\displaystyle\frac{\langle\Psi_{\mathbb{W}}(\tau)|(H-E(\tau))^{2}|\Psi_{\mathbb{W}}(\tau)\rangle}{\langle\Psi_{\mathbb{W}}(\tau)|\Psi_{\mathbb{W}}(\tau)\rangle}\
.$ (24)

For greater clarity, the action of the projector in Eq. (22) is calculated
step by step

$\displaystyle\begin{split}&\frac{dE(\tau)}{d\tau}=\frac{1}{2}\left[\frac{\langle\Psi_{\mathbb{W}}(\tau)|H\left(\mathtt{P}_{\Psi_{\mathbb{W}}(\tau)}\frac{d}{d\tau}|\Psi_{\mathbb{W}}(\tau)\rangle\right)}{\langle\Psi_{\mathbb{W}}(\tau)|\Psi_{\mathbb{W}}(\tau)\rangle}+c.c.\right]\\\
=&-\frac{1}{2}\left[\frac{\langle\Psi_{\mathbb{W}}(\tau)|H\left(\mathtt{P}_{\Psi_{\mathbb{W}}(\tau)}(H-E(\tau))|\Psi_{\mathbb{W}}(\tau)\rangle\right)}{\langle\Psi_{\mathbb{W}}(\tau)|\Psi_{\mathbb{W}}(\tau)\rangle}+c.c.\right]\\\
=&-\frac{\langle\Psi_{\mathbb{W}}(\tau)|(H-E(\tau))\mathtt{P}_{\Psi_{\mathbb{W}}(\tau)}^{2}(H-E(\tau))|\Psi_{\mathbb{W}}(\tau)\rangle}{\langle\Psi_{\mathbb{W}}(\tau)|\Psi_{\mathbb{W}}(\tau)\rangle}\
.\end{split}$ (25)

In the final line, we employed Eq. (18) and leveraged a property of projection
operators, namely,
$\mathtt{P}_{\Psi_{\mathbb{W}}(\tau)}^{2}=\mathtt{P}_{\Psi_{\mathbb{W}}(\tau)}$.
The imaginary time $\tau$ is therefore adjusted and the modified inverse
temperature $\beta(\tau)$ computed by integrating Eq. (22)

$\gamma=\frac{d\beta}{d\tau}\quad\Rightarrow\quad\beta(\tau)=-\int_{0}^{\tau}\frac{1}{\sigma^{2}}\frac{dE}{d\tau^{\prime}}\Big{|}_{\mathbb{W}(\tau^{\prime})}d\tau^{\prime}\
.$ (26)

We emphasize that the temperature adjustment does not fully correct the errors
due to the restriction of the time evolution to the variational manifold
$\mathcal{M}$. The proper state in the full Hilbert space is not known and, in
general, different from the variational state.

###  II.3 Results for Small Systems

To investigate the performance of the proposed algorithm, we use the spin-1/2
Heisenberg chain with Hamiltonian

$H=\sum_{j=1}^{N-1}\left\\{-\frac{1}{2}(S^{+}_{j}S^{-}_{j+1}+h.c.)+S^{z}_{j}S^{z}_{j+1}\right\\}$
(27)

as a test case. We use open boundary conditions. We note that while for
periodic boundary conditions the number of variational parameters in the
neural network can be reduced due to translational invariance, this does
create a new numerical problem for systems where there is no Marshall sign
rule as, for example, for the Hamiltonian in Eq. (27) with a sign flip in the
first term. This, however, is a separate issue which we do not wish to discuss
here. To show that the purification algorithm can indeed properly describe the
thermodynamics of the spin-1/2 Heisenberg chain, we start by considering small
systems where we can work with the full Hilbert space and can directly compare
to exact diagonalization results. These calculations will also already
indicate some issues and differences to zero-temperature calculations. In Fig.
1, we show results for the isotropic Heisenberg model with $N=7$ sites using
the mRBM in the purification algorithm, where we keep all basis states and use
no beta modification. Note that due to the doubling of the number of sites
when adding the ancillas, this already corresponds to a Hilbert space
dimension of $2^{14}$.

Figure 1: Inner and free energy for a Heisenberg chain with $N=7$ sites, a
desired tolerance of $\epsilon=10^{-9}$ and $R=2\cdot 10^{3}$ iterative steps
obtained using the purification algorithm without the beta modification. In
the mRBM ansatz, $M=28$ hidden units are used. The insets show the absolute
errors for $c=0.3,0.4,0.5$.

The results show that for $M=28$ hidden units in the mRBM ansatz (8), the
exact results for the energy $e(T)$ and the free energy $f(T)$ are reproduced
with absolute errors $\sim 10^{-3}-10^{-4}$. Note that we fix the
normalization by setting $\langle\Psi(\infty)|\Psi(\infty)\rangle=2^{N}$ which
gives the proper free energy at infinite temperatures. If the parameter $c$ in
the mRBM ansatz is chosen too small, then we already start with a state whose
energy is strongly deviating from zero—the infinite-temperature value of the
energy for the chosen form of the Hamiltonian in Eq. (27)—and thus a bad
approximation of the infinite temperature state. The course of the imaginary
time in the Heun algorithm is different depending on the respective initial
state and the topology of the area that is traversed in the time evolution
(see Fig. 2).

Figure 2: Imaginary time $\tau$ in the Heun algorithm versus number of
iterative steps for a Heisenberg chain with $N=7$ sites. The parameters are
the same as in Fig. 1. The beta modification is depicted in the inset.

The effect of the beta modification depicted in the lower right inset is
negligible for temperatures larger than the finite size gap and is smaller
than the inherent error of the approach. Thus, it is not visible in Fig. 1.
The lower inset in Fig. 1 shows that the free energy is obtained with a
similar accuracy. For larger systems, it will no longer be possible to use the
full Hilbert space. Instead, a Monte Carlo sampling using the Metropolis
algorithm is employed. To better understand potential additional issues which
may arise due to the sampling, we start by showing the results for the same
system size, $N=7$, in Fig. 3.

Figure 3: Inner and free energy for a Heisenberg chain with $N=7$ sites using
the purification algorithm and a Monte Carlo sampling of the Hilbert space
($M=28$, $4\cdot 10^{4}$ samples, $\epsilon=10^{-7}$, $R=2\cdot 10^{3}$). The
results with beta modification are shown as solid lines, the ones without as
dashed lines. The error in the energy for the case $c=0.2$ is not shown in the
upper inset.

Here, we find that the results for the energy have a similar accuracy as those
for the full Hilbert space simulations if sufficiently many Monte Carlo
samples are used and provided that the initial state is carefully chosen. The
latter turns out to be more intricate here. First, a certain amount of noise
in the $a$ and $W$ parameters is required for the imaginary time evolution to
work properly. We have chosen to draw the real and imaginary parts of both
sets of parameters from independent normal distributions with a width
$w=0.02$. Second, the initial parameters $c_{i}=c=\mbox{const}>0$ have to be
chosen just large enough to obtain
$\langle\Psi(\infty)|H|\Psi(\infty)\rangle\sim 10^{-3}-10^{-4}$. If the
initial $c$ parameters are chosen too large, then the algorithm remains stuck
close to $|\Psi(\infty)\rangle$ for a large number of steps. Additionally, the
sampling fluctuations become large and the algorithm shows instabilities for
low temperatures. By decreasing the tolerance in the adaptive Heun algorithm
and using the beta modification, acceptable results can still be obtained, but
come with the cost of an increase in the number of iterations (see Fig. 4).

Figure 4: Imaginary time $\tau$ in the Heun algorithm versus number of
iterative steps for a Heisenberg chain with $N=7$ sites. The parameters are
the same as in Fig. 3. The beta modification is depicted in the inset.

Furthermore, the approximation of the infinite-temperature state is bounded by
the sampling error. Thus, it is not feasible in this case to use values
$c>0.4$. Making $c$ too small, on the other hand, leads to an initial error
given by the non-zero initial energy which remains present down to low
temperatures.
In practice, we proceed as follows: With $w$ small and fixed we obtain results
for a series of calculations with increasing initial values for $c$ and
different values for the desired tolerance $\epsilon$. The results are in
general comparable for a certain range of initial values for $c$. Furthermore,
we note that the error is dominated by the sampling error which is why an
increase in the number of hidden units $M$ does not have a measurable effect.
This can be deduced as well from the beta modification having no substantial
effect up until low temperatures. Making $M$ too large, generally seems to
lead to additional instabilities. For the considered model, $M\in[N,4N]$ seems
to be optimal. While we find that this approach to find the right initial
parameters and the right number of hidden units does work, we note that it
requires to run multiple simulations with different initial parameters $c$ and
$\epsilon$. We note also, that without being able to compare with other
methods it might not always be easy to decide what the optimal initial
parameters are. Furthermore, we note that the free energy is more sensitive
than the energy to the quality of approximation of the initial state and to
the accuracy of the time evolution (see Fig. 3).

###  II.4 Results for System Sizes Beyond Exact Diagonalization

In this section, we want to briefly show that the purification algorithm can
be scaled up and used to access the thermodynamics of quantum lattice models
for system sizes outside the reach of exact diagonalization. We want to stress
again that the purpose of this paper is a proof of principle demonstration of
using neural networks to investigate thermodynamic properties and the
discussion of certain generic obstacles. There are a number of directions in
which the algorithms presented here could possibly be improved, including but
not limited to the use of more complicated network structures such as CNN’s or
modifications to the imaginary time evolution by using, for example, implicit
iterative methods.
As an example, we present in Fig. 5 results obtained using the purification
algorithm for a Heisenberg chain with 30 sites and compare the results to QMC
results for the energy and the free energy, where the latter is obtained by
the respective integration of the former. We used the stochastic series
expansion code available at [38] which is described in Ref. [39].

Figure 5: Inner and free energy for an open Heisenberg chain with $N=30$ sites
using $4\cdot 10^{4}$ samples and $R=2\cdot 10^{3}$ iteration steps. We choose
$c=0.34$ and $\epsilon=10^{-6}$. Results without the beta modification are
represented by dashed lines. The results for $M=60$ are omitted in the main
plot for visual clarity.

As for the smaller system sizes considered in the previous section, a post-
selection of data is still required and can be done in a systematic manner. We
find again that the parameter range $c=0.3...0.4$ gives initial energies
$|e|\sim 10^{-3}$ and allows for an accurate imaginary time evolution. The
algorithm is executed for different tolerances. For $c=0.34$ and $M=60$, the
results for the energy appear to be converged sufficiently for the beta
modified results to stay within the mentioned accuracy for all temperatures.
Increasing the number of hidden units results in a small improvement for low
temperatures, but instabilities that occur earlier in imaginary time $\tau$.
However, the algorithm is able to proceed faster to reach lower temperatures.
We want to emphasize that the beta modification depicted in Fig. 6 saves us
from increasing the number of hidden units to an extent that is not
numerically feasible anymore. The results without beta modification shown as
dashed lines are about one order of magnitude less accurate than the modified
data. The beta modification is thus a useful tool to achieve a good accuracy
without having to increase the number of hidden units to the point where
numerical instabilities become more and more problematic. The lowest
temperatures reached are limited by the progress of the Heun method and the
beta modification, which at some point result in a temperature which is not
decreasing anymore. To overcome this, increasing the number of hidden units
and the number of sampled states is ultimately inevitable. To obtain a better
approximation of the ground state, the regularization has to be adjusted by
adding a diagonal matrix with entries being around $10^{-4}-10^{-6}$ (see
Ref.[21] and [30]) to the matrix $\boldsymbol{G}$ explained in Eq. (10) . This
is however hardly applicable for the computation of finite temperature results
since it distorts the evolution itself. We find that the free energy is more
sensitive to the sampling error and that the beta modification is not
applicable for low temperatures leading to the large deviation observed in
Fig. 5. An additional issue are jumps in the inner and free energy at
temperatures close or below the finite size gap. The simulation has to be
stopped in this case. Related problems were also observed in Ref. [40], where
it was suggested that they can be mediated by modifying the respective loss
function. This is also the reason why the free energy can not be obtained by
integrating the energy respectively.

Figure 6: Imaginary time $\tau$ in the Heun algorithm versus number of
iterative steps for an open Heisenberg chain with $N=30$ sites. The parameters
are the same as in Fig. 5. The effect of the beta modification is depicted in
the inset.

###  II.5 Correlation Functions

To check how well the proper thermodynamics is approximated beyond the energy,
we compute also the correlation functions $C^{d}_{\parallel}=\langle
S_{k}^{z}S_{k+d}^{z}\rangle$ and $C^{d}_{\perp}=-\frac{1}{2}\langle
S_{k}^{+}S_{k+d}^{-}+S_{k}^{-}S_{k+d}^{+}\rangle$. In Fig. 7, the results for
$k=10$ and $d=2$ are compared to the QMC results and agree to an accuracy of
order $10^{-3}$. The deviations observed in Fig. 7 are dominated by the Monte
Carlo sampling error. We note that the SU(2) symmetry is broken, in general,
but that it is restored at low temperatures.

Figure 7: Correlation functions $C^{2}_{\parallel,\perp}$ for a chain of
length $N=30$ obtained using the purification method are compared to the QMC
result. The parameters are the same as in Fig. 5 and $M=60$. The inset shows
$|C^{2}_{\parallel}-\frac{1}{2}C^{2}_{\perp}|$ and the deviations of the two
correlation functions from the QMC result.

###  II.6 Real-time Dynamics at Finite Temperatures

Among the most interesting questions related to using neural networks as
approximators for quantum wave functions, is whether or not this approach will
lead to major progress in the calculation of dynamical correlation functions.
Numerically, real-time evolution is a difficult problem because QMC does only
give direct access to imaginary times and analytical continuations of
numerical data to real times are not unique. In time-dependent DMRG
algorithms, on the other hand, the growth in time of the entanglement entropy
often severely limits the accessible simulation times. DMRG is based on a
truncation of the Hilbert space. If the reduced density matrix $\rho_{A}$ has
a maximum dimension $\chi$, then the maximal entanglement entropy which can be
faithfully represented is
$S^{\textrm{max}}_{\textrm{ent}}=-\text{tr}\rho_{A}\ln\rho_{A}=\ln\chi$. If
the entanglement entropy grows linearly in time—which is the generic case
following a quantum quench from a product state and for equilibrium
correlations at high temperatures, see below—leading to a volume law
$S_{\textrm{ent}}\sim N$ at long times, then one needs
$\chi\sim\min\\{\exp(t),\exp(N)\\}$ states in the truncated Hilbert space.
This exponential barrier often prevents to fully address dynamical problems.
In neural network algorithms, on the other hand, the Hilbert space is not
explicitly truncated and the hidden units typically are fully connected to all
units in the input layer. This means that in principle, states with large
entanglement can already be represented by neural networks with a small number
of hidden units. Indeed, it was shown in Ref.[41] that a volume law entangled
state can be represented with $M\sim N$ hidden units only. This, however, does
not answer the question whether or not the dynamics of a quantum system can be
efficiently simulated in general using a neural network representation of the
wave function. While it is an interesting theoretical result that specific
volume-law entangled states exist which do have an efficient representation in
this formalism, much more is required to investigate the long-time dynamics of
a quantum system. Here, typically all the eigenstates of the system will
contribute with most of them—in a clean, non-disordered system—showing volume-
law entanglement. Only if a large fraction of these states can be efficiently
represented by a given network can we hope to obtain accurate results for
dynamical quantities at long times.
In the following, we want to numerically investigate how many hidden units are
required to obtain accurate results for the spin-spin autocorrelation function

$\langle
S^{z}_{n}(0)S^{z}_{n}(t)\rangle=\frac{1}{Z}\text{tr}\left\\{S^{z}_{n}\text{e}^{itH}S^{z}_{n}\text{e}^{-itH}\right\\}$
(28)

at site $n$ of a Heisenberg chain at infinite temperatures using the
purification algorithm. In addition, we also calculate the entanglement
entropy

$S=-\text{tr}\rho_{A}\ln\rho_{A}\ ,$ (29)

where $\rho_{A}$ is the reduced density matrix for a subsystem $A$ consisting
of half of the chain of real spins and ancilla spins. We begin with results
for $N=4$ shown in Fig. 8.

Figure 8: Top: Autocorrelation function $\langle
S^{z}_{n}(0)S^{z}_{n}(t)\rangle$ for $n=N/2$ at infinite temperatures for a
Heisenberg chain with OBC and $N=4$ sites. Bottom: Corresponding entanglement
entropy for cutting the chain of real and ancilla spins into two equal halves.
The RBM results for different $M$ are compared to the exact result.

While the results are numerically stable even for $M=8$ hidden units, we see
that we need at least $M=16$ to accurately approximate the exact result. From
the panel showing the entanglement entropy, it is clear that while $M=8$ is
sufficient to represent certain states with volume-law entanglement [41], it
does not capture the full entanglement structure of the time evolved state for
times $Jt\gtrsim 2$.
Next, we show the same results for $N=6$ in Fig. 9.

Figure 9: Same as Fig. 8 for a chain with $N=6$ sites. The results for
$M=24,48$ are only shown up to $Jt\sim 13$.

In this case, at least $M=96$ hidden units are required. We note again that
already for $M=12,24,48$ the entanglement entropy reaches values at long times
of the order of the maximum entanglement in the time evolved state. However,
the proper time evolution is only captured for $M\geq 96$. Lastly, we also
show in Fig. 10 results for $N=8$.

Figure 10: Same as Fig. 8 for a chain with $N=8$ sites. The corresponding
entanglement entropies have not been calculated.

Note that due to the ancillas, this already corresponds to a system with $16$
sites and we do keep the full Hilbert space. Since the calculation of the
entanglement entropy is numerically expensive, we only show the correlation
function. From the results it is clear that even $M=256$ hidden units are not
sufficient.
To summarize, the number of hidden units $M$ required to accurately simulate
the infinite-temperature autocorrelation (28) increases much faster than
linear. Assuming that $M=512$ hidden units give reasonably converged results
for $N=8$, the increase follows in the best case scenario a power law $M\sim
N^{\alpha}$ with $\alpha\gtrsim 4$ but is more consistent with an exponential
scaling. While a potential power-law increase in the hidden units with system
size would make this approach in theory more powerful than time-dependent
DMRG, it might be very difficult in practice to benefit from this scaling (if
it is not exponential to begin with as seems more likely): while the number of
states kept in a time-dependent DMRG calculation is only limited by the
available memory and the compute time to perform a single-value decomposition
or a diagonalization allowing to reach truncated dimensions of order $\chi\sim
10^{5}-10^{6}$, time-dependent neural networks for $N\sim 10^{2}$ with $M\sim
10^{8}$ hidden units might be difficult to optimize given that the Hilbert
space needs to be sampled. We also note that in DMRG the dimension $\chi$ can
be adaptively increased with simulation time. Based on the neural network
results shown in Figs. 8, 9, 10 it seems, on the other hand, that the
calculation becomes inaccurate already at short times if an insufficient
number of hidden units is used. Furthermore, there does not seem to be an easy
way to adaptively increase the number of hidden units during a calculation.
Finally, we note that the strong, likely exponential increase of the number of
required hidden units is consistent with a recent publication where similar
behavior was observed for quench dynamics [42].

##  III Sampling

To compute the ensemble average of a quantum mechanical operator A (see Eq.
(1)) at infinite temperature, one needs to compute the expectation value with
respect to every vector of an orthonormal set
$\\{|\phi_{n}\rangle\\}_{n=1}^{2^{N}}$

$\langle
A\rangle_{\beta=0}=\frac{\text{tr}[A]}{2^{N}}=\frac{1}{2^{N}}\sum_{n=1}^{2^{N}}\langle\phi_{n}|A|\phi_{n}\rangle\
.$ (30)

The idea is now that instead of working with a full orthonormal set we can
also sample the trace using random states. These are constructed by choosing
$2^{N}$ random vector components $c_{n}$ from a symmetric probability
distribution ([32]) to form the infinite-temperature states

$|\Psi(\beta=0)\rangle=\sum_{n=1}^{2^{N}}c_{n}|n\rangle\ .$ (31)

It is essential that the vector components are uncorrelated. The ensemble
average is then approximated by using $S$ realizations of those states.

$\frac{1}{S}\sum_{s=1}^{S}\langle\Psi_{s}(\beta=0)|A|\Psi_{s}(\beta=0)\rangle\approx
E(|c_{n}|^{2})\ \text{tr}(A)$ (32)

For normalized states, the expectation value of any $c_{n}$ is
$E(|c_{n}|^{2})=\frac{1}{2^{N}}$ and thus the correct value of the ensemble
average at infinite temperature is obtained. The error goes in general like
$\frac{1}{\sqrt{S}}$. An upper bound is computed in Ref. [32] in terms of the
operator $A$. To include this idea into a variational algorithm using the RBM
ansatz, the parameters of the NN have to be chosen properly. States that
fulfill the above criteria and are efficiently represented by the RBM are, for
example, the so called minimally entangled typical thermal states (METTS)
[43]. The word ‘typical’ will later be defined and used in Sec. IV in the
sense of quantum typicality. To avoid confusion, it is dropped in this work
for the notation of these states.

###  III.1 Minimally Entangled Thermal States (METS)

The Minimally Entangled Thermal States (METS) were first discussed as a tool
to perform imaginary time evolution in a standard DMRG algorithm [43]. In Ref.
[36], the approach was then generalized for the NN ansatz. A random METS is
defined as

$|\Psi_{0}\rangle=\bigotimes_{j=1}^{N}\left(\frac{\xi_{j,\uparrow}|\uparrow\rangle_{j}+\xi_{j,\downarrow}|\downarrow\rangle_{j}}{\sqrt{|\xi_{j,\uparrow}|^{2}+|\xi_{j,\downarrow}|^{2}}}\right),$
(33)

where the vector $\boldsymbol{\xi}$ is chosen according to a standard complex
normal distribution $\mathbb{C}N(\mu=0,C=\mathbb{1})$ with $\mu$ the mean and
$C$ the covariance matrix of the distribution

$P_{0}(\boldsymbol{\xi})=\left(\frac{1}{\pi}\right)^{2N}\exp\left(-\sum_{j=1}^{N}\sum_{\sigma=\uparrow,\downarrow}|\xi_{j,\sigma}|^{2}\right)\,.$
(34)

This is equivalent to choosing the real and imaginary part according to a real
normal distribution $N(\mu=0,\sigma=1/2)$ with $\sigma$ being the variance. A
METS state can be obtained in the RBM, Eq. (8) with $\boldsymbol{c}=0$, up to
an irrelevant overall constant by setting

$a_{j}=-\frac{1}{2}\ln\left(\frac{\xi_{j,\uparrow}}{\xi_{j,\downarrow}}\right)\
,$ (35)

$\boldsymbol{W}=0$ and $\boldsymbol{b}=0$, see Ref. [36]. These states will be
referred to as METS-NN from here on. Numerically, one needs to add some noise
to the NN parameters $\boldsymbol{W}$ and $\mathbf{b}$ as in the purification
method. We decided to choose the real and imaginary parts of both of them from
a normal distribution centered around zero with $\sigma_{NN}=1/N$ so that the
variance of the wave function components stays constant for variable system
size.

###  III.2 Results for the Energy

The set of METS-NN states at infinite temperature is evolved by the adaptive
Heun method as explained in Sec. II.1. The imaginary time of the algorithm
$\tau$ is afterwards modified as explained in Sec. II.2. The results for the
energy of a Heisenberg chain of length $N=30$ and $S=30$ states are depicted
in Fig. 11. The accuracy is here also of order $10^{-3}$ which is expected
from the sampling error. We also note that the temperatures reached using the
purification method are lower than the ones reached here. This is depicted in
Fig. 12 and due to the already mentioned convergence issues at low
temperatures (see Sec. II.4). Compared to the purification method, this
problem is compounded here because one needs to simulate $\sim 30$ states to
low temperatures without any such instabilities occurring. Additionally, the
state in which these instabilities occur first defines the minimum temperature
that can be reached. To allow for the computation of the free energy by
integrating the energy results respectively, as it was done for the QMC
results in Fig. 5, the evolution has to reach lower temperatures. This can
only be achieved by increasing the number of hidden units drastically.

Figure 11: Inner energy of the METS-NN ansatz compared to the QMC solution for
an open Heisenberg chain of length $N=30$ using $4\cdot 10^{4}$ Metropolis
samples, $R=2\cdot 10^{3}$ and $\epsilon=10^{-6}$. The average is taken over
$S=30$ initial states. The dashed lines show the data without beta
modification. Figure 12: Imaginary time $\tau$ in the Heun algorithm versus
number of iterative steps for the METS-NN ansatz for an open Heisenberg chain
of length $N=30$ and parameters defined in Fig. 11. The beta modification is
shown in the inset.

###  III.3 Correlation Functions

For the sampling method, it is checked as well how accurate thermodynamics is
approximated beyond the energy by computing spin-spin correlators like
$C^{d}_{\parallel}$ and $C^{d}_{\perp}$ for $k=10$ (see Sec. II.5). These are
not obtained as an average over local operators such as the energy where an
average over $N-1$ bond energies is taken. Thus, the METS-NN algorithm results
in larger errors for single spin-spin correlators. In Fig. 13, the results for
$C^{2}_{\parallel}$ and $C^{2}_{\perp}$ are compared to the QMC results. It is
simply no longer sufficient to use $\sim 30$ states only. Instead, to obtain a
result that is consistent within the Monte Carlo sampling error, $\sim 1000$
states are required. Note that a preselection of states would break the
relation in Eq. (32). We emphasize that for any specific single state
$C^{d}_{\parallel}$ and $C^{d}_{\perp}$ vary, with neither of them being more
accurate in general. The results consequently also violate $SU(2)$ invariance.

Figure 13: Same as Fig. 7 using the METS-NN states with the parameters
specified in Fig. 11.

##  IV Quantum Typicality

The third method investigated in this paper promises to provide one state that
is sufficient for the computation of finite temperature expectation values for
large system sizes. The definition of thermal pure quantum states (TPQS) or
quantum typical states comes without the equal a priori probability postulate
of statistical mechanics and replaces the ensemble average by pure state
expectation values (see Ref. [44] and [45]). Quantum typical states were first
discussed for the canonical situation of a small system S coupled to an
environment E. In the standard microcanonical formulation, the coupled system
is assumed to be in an equiprobable state $\Omega$ corresponding to a definite
energy E. In the TPQS formulation, the coupled system is instead assumed to be
in a pure state $\rho$. It is shown in Ref. [44] and [45] using Levy’s Lemma
that the mixed states of the system S—obtained in both cases by tracing out
the environment E—are nearly indistinguishable $\rho_{S}\approx\Omega_{S}$.
This is founded on properties of high-dimensional normed spaces [46] and
especially the theorem of measure concentration [47]. How this idea can be
used to construct a TPQS at arbitrary temperature is discussed in Refs. [35,
48], culminating in the following general definition for a TPQS: Let
$\\{|n\rangle\\}_{s=1}^{K}$ be an orthonormal basis of the Hilbert space
$\mathcal{H}$ and

$|\Psi(\beta=0)\rangle=\sum_{n=1}^{K}c_{n}|n\rangle$ (36)

a normalized pure state with components
$c_{n}=\frac{x_{n}+iy_{n}}{\sqrt{2}}\in\mathbb{C}$, where the real parameters
$x_{n}$ and $y_{n}$ are drawn from the unit normal distribution. This state
can be evolved in imaginary time

$|\Psi(\beta)\rangle=\exp\left(\frac{-\beta H}{2}\right)|\Psi(0)\rangle\ .$
(37)

The ensemble average and the pure state expectation value of an operator $A$
are then defined by

$\langle A\rangle_{\beta}^{\textrm{ens}}=\frac{\text{tr}[Ae^{-\beta
H}]}{Z}\quad,\quad\langle
A\rangle_{\beta}^{\textrm{TPQS}}=\langle\Psi|A|\Psi\rangle\ ,$ (38)

respectively. The operator A is assumed to be a low-degree polynomial
($\text{degree}\leq m$) of local operators [49]. If $\forall\epsilon>0$ a
function $\eta_{N}(\epsilon)$ exists with
$\eta_{N}(\epsilon)\xrightarrow[]{N\rightarrow\infty}0$ such that the
probability

$P\left(|\langle A\rangle_{\beta}^{TPQS}-\langle
A\rangle_{\beta}^{ens}|\geq\epsilon\right)\leq\eta_{N}(\epsilon)\ ,$ (39)

then $|\Psi\rangle$ is called a TPQS [35]. One can show that the right hand
side is for normally distributed $c_{n}$ bounded from above by

$\eta_{N}(\epsilon)\leq\frac{N^{2m}e^{-\mathcal{O}(N)}}{\epsilon^{2}}\ .$ (40)

This means that $\eta_{N}(\epsilon)$ is decreasing exponentially in system
size $N$. A similar result is obtained in Ref. [33] and even extended to
dynamical correlation functions [34]. For large system sizes, it is therefore
sufficient to consider only a few or even a single pure state. In Ref. [50],
it was investigated thoroughly which minimum requirements the randomly chosen
wave function components of the state $\Psi(\beta=0)$ need to fulfill for the
state to be a TPQS resulting in the following two constraints:

  * •

The components $c_{n}$ are uncorrelated.

$\overline{c_{m}^{*}c_{n}}=\overline{|c_{n}|^{2}}\delta_{n,m}$ (41)

  * •

The mixed state density matrix $\rho$ has low purity

$\text{tr}\rho^{2}=\sum_{n}\rho_{n}^{2}\ll 1\ ,$ (42)

where

$\rho_{n}=\overline{|c_{n}|^{2}/\sqrt{\sum|c_{n}|^{2}}}\ .$ (43)

This constraint guarantees that averages are not dominated by just a few large
$\rho_{n}$.

Some probability distributions like the normal and the uniform distribution
that fulfill these two conditions have been discussed in Refs. [35, 49]. In
this work, we are, however, interested in using a variational ansatz like the
NN ansatz to compute expectation values of operators at finite temperatures.
Unfortunately, a NN ansatz in general and the RBM ansatz

$\Psi_{\mathbb{W}}(\mathbf{s})=\
\exp\left(-\sum_{i=1}^{N}a_{i}s_{i}\right)\prod_{j=1}^{M}\cosh\left(b_{j}+\sum_{i=1}^{N}W_{ij}s_{i}\right)$
(44)

in particular are problematic starting points. The non-linear activation
function leads, in general, to a small number of amplitudes $|\Psi(s)|$
dominating the wave function, thus violating the second condition stated above
for a typical state. Additionally, it was proven in Ref. [51] and discussed in
Ref. [52] and [53] that TPQS have volume-law entanglement ruling out the METS-
NN as possible candidates.

###  IV.1 Pair-Product (PP) Wave Function

To fulfill all the mentioned restrictions for a possible variational TPQS, the
NN (RBM) is complemented by a Pair-Product (PP) wave function component

$|\Phi_{PP}\rangle=P_{G}^{\infty}\left[\sum_{i,j=0}^{N}\sum_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}F_{ij}^{\sigma,\sigma^{\prime}}c_{i,\sigma}^{\dagger}c_{j,\sigma^{\prime}}^{\dagger}\right]^{\frac{N_{e}}{2}}|0\rangle\
,$ (45)

where $P_{G}^{\infty}=\prod_{i}^{N}(1-n_{i,\uparrow}n_{i,\downarrow})$ is the
Gutzwiller factor prohibiting double occupancy, $c_{i,\sigma}^{\dagger}$
($c_{i,\sigma}$) the creation (annihilation) operator of a fermion with spin
$\sigma$, $F_{ij}^{\sigma,\sigma^{\prime}}$ are complex parameters,
$n_{i,\sigma}=c_{i,\sigma}^{\dagger}c_{i,\sigma}$ the occupation number
operator and $N_{e}$ is the number of electrons which has to be even (see
Refs. [54, 30]). For the Heisenberg model, this translates to $N_{e}=N$. It
should be mentioned that the PP wave function can be generalized to an odd
number of electrons $N_{e}$ [55], however, since we are interested in
expectation values in the thermodynamic limit, it is more convenient to
consider spin chains of even length $N$. The PP-ansatz stems from more general
considerations of variational wave functions for lattice models in the
fermionic formulation [56]. The two approaches are easily combined since
$s_{i}^{z}=\frac{1}{2}(n_{i,\uparrow}-n_{i,\downarrow})$ and there is no
double occupation allowed here. This leads to the following PP-NN-ansatz for
the overall wave function

$|\Psi\rangle=\sum_{\mathbf{s}}\underbrace{\Psi_{\mathbb{W}}(\mathbf{s})\Phi_{PP}(\mathbf{s})}_{\Psi(s)}|\mathbf{s}\rangle\
.$ (46)

This ansatz was first considered in Ref. [30] to improve the accuracy of the
ground state of the Heisenberg model and later used to compute the phase
diagram of the $J_{1}-J_{2}$ square lattice Heisenberg model [31]. The overlap
of the PP wave function with a real space configuration
$|\boldsymbol{\sigma}\rangle=c_{1,\sigma_{1}}\cdot\cdot\cdot
c_{N,\sigma_{N}}|0\rangle$ is given by the Pfaffian

$\Phi_{PP}(\boldsymbol{\sigma})=(N_{e}/2)!\ Pf(X[\boldsymbol{\sigma}])$ (47)

with the skew-symmetric matrix

$X_{ij}^{\sigma_{i}\sigma_{j}}=F_{ij}^{\sigma_{i}\sigma_{j}}-F_{ji}^{\sigma_{j}\sigma_{i}}\
,$ (48)

and

$Pf(X[\boldsymbol{\sigma}])=\frac{1}{2^{N/2}(N/2)!}\sum_{\alpha\in
S_{N}}\text{sgn}(\alpha)\prod_{i=1}^{N/2}X_{\alpha(2i-1)\alpha(2i)}[\boldsymbol{\sigma}]\
.$ (49)

The real and imaginary parts of the initial variational parameters
${X_{ij}^{\sigma_{i}\sigma_{j}}=F_{ij}^{\sigma_{i}\sigma_{j}}-F_{ji}^{\sigma_{j}\sigma_{i}}}$
are chosen randomly from a normal distribution centered around zero. This
means that there are $4$ parameter arrays ($X^{\uparrow\uparrow}$,
$X^{\uparrow\downarrow}$, $X^{\downarrow\uparrow}$,
$X^{\downarrow\downarrow}$) each containing $(N^{2}-N)/2$ different parameters
and that the initial values for the amplitudes
$\Phi_{PP}(\boldsymbol{\sigma})$ are also randomly distributed around zero. We
also consider a state, called PP-NN-local in the following, where the
projected pairs are limited to nearest neighbors only
($F_{ij}^{\sigma,\sigma^{\prime}}=0$ if $|i-j|>1$).
If one chooses the initial parameters of the RBM to be zero, or at least for
numerical convenience so small that the variance of the NN part of the wave
function is small compared to the one of the PP part, the infinite-temperature
state is solely described by the PP part. The correlation of the wave
functions of two different spin configurations in the PP part is computed to
be
$\overline{\Phi_{PP}(\boldsymbol{\sigma}_{k})\Phi_{PP}(\boldsymbol{\sigma}_{m})}=0$
since two configurations differ in every term of the Pfaffian by at least one
random variable $X_{ij}^{\sigma_{i}\sigma_{j}}$. As long as the number of
parameters is smaller than the Hilbert space dimension, non-zero higher order
correlations will exist. However, we expect that such correlations have no
influence on the expectation values of local observables. I.e., we expect that
the difference in expectation values obtained by an actual TPQS given, for
example, by normally distributed $c_{n}$ and those obtained by the variational
PP-NN ansatz in Eq. (46) will be small if a sufficient number of parameters is
kept.

###  IV.2 The Infinite-Temperature TPQS

First, we check if the states chosen according to the PP-NN ansatz, see Eqs.
(46)-(49), are good approximations of an infinite-temperature TPQS. As
comparison, the normally distributed states chosen according to Eq. (36) and
the METS-NN states chosen according to the NN ansatz, see Eqs. (33)-(35), are
used. Second, the entanglement entropy of those states as a measure of the
approximation complexity is also computed.

####  IV.2.1 Comparison of Typicality Features

It is not a priori clear to which extent the conditions for a TPQS, discussed
at the beginning of Sec. IV, are fulfilled by the PP-NN ansatz. Additionally,
apart from the entanglement entropy, it is not clear why the METS-NN states do
not describe a TPQS. We investigate first to what extent Eq. (41) is fulfilled
compared to a TPQS with normally distributed random numbers for the wave
function components - from here on referred to as the Gaussian case. It is
emphasized that in the PP-NN case the real and imaginary parts of all the NN
parameters are not chosen as described in Sec. III.1 but from a normal
distribution centered around zero with $\sigma_{NN}=1/N$ while those of the PP
parameters are always chosen from a normal distribution centered around zero
with $\sigma_{PP}=({N/2})!\cdot\frac{2^{N/2}}{(N!)^{2/N}}$ to keep the
variance of both wave function components (see Eq. (46)) constant for growing
system size. This way, the NN parameters in the PP-NN ansatz do not contribute
to the properties of the infinite-temperature state but are only relevant for
the numerics of the imaginary time evolution. Thus, those parameters can be
neglected in the following discussion. To check the validity of Eq. (41), the
square root of the mean of the absolute value of the pairwise Pearson
correlation coefficient matrix is compared in Fig. 14 for varying number of
states $n$. The full Hilbert space is used for a spin chain of length $N=12$.
Of course, the diagonal elements are excluded. The results for the Gaussian
case show as expected a $n^{-1/2}$ behavior as is checked by a respective fit.
While this also holds for the PP ansatz up to an overall constant, the PP-
local and the METS-NN ansatz differ significantly from this behavior.

Figure 14: Square root of the mean of the absolute of the pairwise Pearson
correlation coefficient matrix of the wave function components for $n=10^{3}$
different states and a spin chain of length $N=12$.

To check what this means for the distribution of expectation values, the
standard deviation of $\langle S_{0}^{z}\rangle$ is computed for $n=10^{3}$
sets of parameters $\mathbb{W}$ and varying length $N$ in Fig. 15. The
expectation values are obtained for every set of parameters by using $10^{4}$
Metropolis samples. Since the Gaussian case reaches values affected by the
sampling error, the full Hilbert space is used. While the Gaussian case shows
exponential decay to the infinite-temperature result, the data for the PP
ansatz can be fitted very well by a function that is proportional to
$e^{-a\sqrt{N}}$, where $a$ is a fitting parameter. However, the METS-NN and
the PP-local ansatz show no decay with system size at all since they are
product states and thus only depend on the parameters connected to the 1st
spin. This is as well as the result for the Gaussian case analytically
verified by computing the probability distribution of $\langle
S_{0}^{z}\rangle$ which is in all those three cases given by

$2\langle S_{0}^{z}\rangle=\frac{A}{A+B}-\frac{B}{A+B}\ ,$ (50)

where

$A=\sum_{k}^{K/2}|z_{k}^{\uparrow}|^{2}\qquad
B=\sum_{k}^{K/2}|z_{k}^{\downarrow}|^{2}$ (51)

and the real and imaginary part of $z_{k}^{\uparrow}$ and $z_{k}^{\downarrow}$
are chosen randomly from a normal distribution. The actual variance of this
normal distribution is irrelevant in this case because it cancels out and is
thus chosen to be one. Therefore, $A$ and $B$ follow a $\chi_{K}^{2}$
distribution. Each term in Eq.(50) is computed to follow a probability
distribution $p(y)$ with $y\in[0,1]$. These, the respective value for $K$ and
the resulting variances are noted in the following table and plotted as solid
lines in Fig. 15.

Ansatz | K | $p(y)$ | $\sigma[\langle S_{0}^{z}\rangle]$
---|---|---|---
Gaussian | $2^{N}$ | $\frac{\Gamma(2^{N})}{\Gamma^{2}(2^{N-1})}y^{2^{N-1}-1}(1-y)^{2^{N-1}-1}$ | $\sqrt{\frac{2^{N-1}+1}{2^{N+1}+2}-\frac{1}{4}}$
PP-local | 4 | $6y(1-y)$ | $\frac{1}{\sqrt{20}}$
METS-NN | 2 | $1$ | $\frac{1}{\sqrt{12}}$

The promising relations for the Gaussian case in Eqs. (39) and (40) thus do
not translate to the discussed variational wave functions with a finite number
of parameters. While the exponential decay is weakened for the PP ansatz,
there is no decay at all with the METS-NN and the PP-local ansatz. For
operators like the Hamiltonian in Eq.(27), there is still a decay of the
standard deviation in all cases for larger spin chains due to the benefit of
averaging over local spin-spin operators along the chain. Naturally, the
results for the mean of the expectation value converge in all cases to the
infinite temperature value with increasing number of parameter sets (see
Ref.[32] and the discussion in Sec.III.1).

Figure 15: Standard deviation of the operator $\langle S_{0}^{z}\rangle$
computed using $10^{4}$ Metropolis samples and $n=10^{3}$ different sets of
parameters chosen as discussed above.

####  IV.2.2 Entanglement Entropy of the Infinite-Temperature TPQS

An important difference between the ansatz wave functions is observed when
considering the entanglement entropy. While the PP-NN ansatz shows a volume
law scaling due to the PP part of the wave function, the METS-NN ansatz shows
no entanglement because the wave function is a product state. The different
scaling with the length of the chain is illustrated for the full Hilbert space
in Fig. 16. This is very relevant since states with volume law entanglement
can in general not be efficiently represented. The RBM is only capable of
doing so for some specific states [41]. Thus, the METS-NN ansatz defines an
initial state that is much easier to approximate. To allow for a similar state
in the PP-NN ansatz, the PP-NN-local ansatz is used. This ansatz shows an
entanglement entropy that is only nonzero if there is a PP pair that crosses
the cut in the middle of the spin chain. This is only the case for chain
lengths for which $N/2$ is odd. By increasing the support of the PP ansatz,
one can obtain states whose properties lie in between the blue and magenta
curves in Fig. 15 and Fig. 16.

Figure 16: Entanglement Entropy of the PP, the PP-local and the METS-NN
initial states compared to the Gaussian case for different chain lengths of an
open Heisenberg chain.

###  IV.3 The variational imaginary time evolution

To compare the PP-NN and the PP-NN-local ansatz to the previous methods, we
average in each case over $30$ initial states and time evolve these according
to Eq. (11). The averaging is still useful due to the fact that the error is
still dominated by the lack of true quantum typicality over the sampling error
(see Fig. 15). We do so again for a Heisenberg chain of length $N=30$. The
results are shown in Fig. 17 and compared to the same QMC results as before.
The number of samples and the number of iteration steps $R$ is adjusted due to
the larger computational cost for the updates of the PP parameters. For the
computation of the Pfaffian after a spin flip, an effective update rule based
on Cayley’s identity is used [54].

Figure 17: Inner energy of the PP-NN and the PP-NN-local ansatz compared to
the QMC solution for an open Heisenberg chain of length $N=30$ using $10^{4}$
Metropolis samples, $R=10^{3}$ and $\epsilon=10^{-6}$. The average is taken
over $30$ initial states. The dashed lines show the data without beta
modification. The PP-NN-local results with beta modification are omitted
because they are not distinguishable from the QMC results on this scale.

The energy values for high temperatures are for the PP-NN ansatz more accurate
as is expected from the discussion Sec. IV.2.1. However, during the course of
the evolution, the PP-NN ansatz states tend to deviate further from the QMC
results. The maximum deviation is roughly at the point where the heat capacity
is maximal. The most elaborate ansatz with the largest number of parameters
produces the most inaccurate results. This is a consequence of the large
entanglement of the initial states $|\Psi(0)\rangle$ in the PP-NN ansatz.
Therefore, it is more difficult to approximate these states and the ones in
their vicinity, although the number of parameters in the PP-NN ansatz is much
larger compared to the METS-NN ansatz. It is stressed that this is not due to
the smaller number of samples, as is evident from the comparison to the PP-NN-
local ansatz with zero entanglement entropy. This results in a much better
approximation. The difference to the METS-NN ansatz is relatively small.
Summarizing, there is no numerical benefit of using TPQS’s in the PP-NN
ansatz. The theoretical benefit is eaten up by the larger entanglement
entropy, leading to a trajectory in the variational manifold that includes
states that are harder to represent. It might be possible to find an ansatz
that can represent the highly entangled states on the trajectory in special
cases, for example for a certain model. While it seems unlikely that this can
be done in general, deep neural networks could possibly lead to a better
representation of such states and thus alleviate this issue to some degree
[57].
The modification of the imaginary time shown in Fig. 18 illustrates the issue
with the variational TPQS in a different way. Especially for high
temperatures, the imaginary time evolution in the PP-NN ansatz shows a very
large beta modification, which means that the proper evolution and the
evolution projected onto the variational manifold deviate strongly. Although
this feature turns somewhat around for lower temperatures, the PP-NN ansatz
states appear to be not a good starting point due to their large initial
entanglement. On the other hand, the error of the time evolution after the
beta modification is applied to the PP-NN-local states is of the order
$10^{-3}$ which is thus again consistent with the Monte Carlo sampling error.

Figure 18: Imaginary time $\tau$ in the Heun algorithm versus number of
iterative steps for the PP-NN and the PP-NN-local ansatz for an open
Heisenberg chain of length $N=30$ and parameters defined in Fig. 17. The beta
modification is shown in the inset.

##  V Conclusions

We have shown that the energy of an open Heisenberg chain can be computed for
a large temperature range with an accuracy limited by the sampling error
(errors of the order of $10^{-3}$ for an $N=30$ site chain) with both the
purification and the sampling approach. The results for the free energy are
more sensitive and can be computed with the purification method with an
accuracy which is roughly one order of magnitude worse than for the energy.
While we have been concerned with a proof-of-principle investigation of these
algorithms here, there are no fundamental obstacles in scaling up both
algorithms to larger system sizes.
In the quantum typicality approach, there are two competing requirements that
could not be fulfilled: The state should be quantum typical and all states on
the trajectory defined by the evolution in imaginary time should be
efficiently represented by the variational ansatz. While the PP-NN ansatz
fulfills the former in a weaker sense, it fails to meet the latter to an
extent that the beta modification cannot fully correct. As discussed in Sec.
IV, it seems questionable if deep network structures can fully cure these
issues. Additionally, the computational cost of the combined PP-NN ansatz is
of course considerably higher compared to the variational wave functions used
in the sampling method. Not only because of the much larger number of
parameters, but also because of the computation of the spin update for the
matrices $X$ in Eq.(49) which is, despite using sophisticated updating
techniques [54], more computationally time-consuming than the update of the NN
part.
The two most promising methods have distinct up- and downsides with respect to
the required computation time. In the purification method, it is sufficient,
after making sure that the parameter $c$ is within a certain range, to
simulate a single state that is part of a Hilbert space of size $2^{2N}$. In
the sampling approach, on the other hand, it is necessary to simulate $\gtrsim
30$ states in a much smaller Hilbert space of size $2^{N}$. Using the fact
that these states can be time evolved independently of each other, the
computational time is significantly lower when using a high-performance
cluster. However, the purification method gives a very accurate approximation
of the infinite-temperature state and even allows for an accurate computation
of correlation functions. In the sampling method, on the other hand, the
number of states needs to be increased considerably to also obtain accurate
results for local correlation functions.
From the attempt to calculate dynamical correlation functions one can learn
that while certain specific highly entangled states can be efficiently
represented as a neural network, this is not true for generic highly entangled
states. We observe that the number of hidden units has to be increased in a
way which is consistent with an exponential scaling with system size. It thus
appears unlikely that neural networks will significantly improve our abilities
to calculate dynamical correlation functions, at least not in one dimension.
The most fruitful applications of the algorithms discussed in this work might
be to the thermodynamics of two- and three-dimensional frustrated spin models.
Here, the standard QMC approach fails and DMRG calculations are limited to
small clusters. In these cases, it seems possible that a neural-network ansatz
using either the purification or the sampling method might outperform both QMC
and DMRG and might thus provide further insights into the physics of
frustrated spin systems.

## References

  * Gubernatis  _et al._ [2016] J. Gubernatis, N. Kawashima, and P. Werner,  _Quantum Monte Carlo Methods: Algorithms for Lattice Models_ (Cambridge University Press, 2016).
  * Sandvik [2019] A. Sandvik,  _Many-Body Methods for Real Materials_, edited by E. Pavarini, E. Koch, and S. Zhang (Forschungszentrum Julich, 2019).
  * Nightingale and Umrigar [1999] M. P. Nightingale and C. J. Umrigar,  _Quantum Monte Carlo methods in physics and chemistry_ (Springer Dordrecht, 1999).
  * Becca and Sorella [2017] F. Becca and S. Sorella,  _Quantum Monte Carlo Approaches for Correlated Systems_ (Cambridge University Press, 2017).
  * Silver  _et al._ [1990] R. N. Silver, D. S. Sivia, and J. E. Gubernatis, Maximum-entropy method for analytic continuation of quantum Monte Carlo data, Phys. Rev. B 41, 2380 (1990).
  * Schollwöck [2005] U. Schollwöck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011), january 2011 Special Issue.
  * Wang and Xiang [1997] X. Wang and T. Xiang, Transfer-matrix density-matrix renormalization-group theory for thermodynamics of one-dimensional quantum systems, Phys. Rev. B 56, 5061 (1997).
  * Bursill  _et al._ [1996] R. J. Bursill, T. Xiang, and G. A. Gehring, The density matrix renormalization group for a quantum spin chain at non-zero temperature, Journal of Physics: Condensed Matter 8, L583 (1996).
  * Sirker and Klümper [2002a] J. Sirker and A. Klümper, Thermodynamics and crossover phenomena in the correlation lengths of the one-dimensional t-J model, Phys. Rev. B 66, 245102 (2002a).
  * Sirker and Klümper [2002b] J. Sirker and A. Klümper, Temperature-driven crossover phenomena in the correlation lengths of the one-dimensional t-J model, Europhys. Lett. 60, 262 (2002b).
  * Sirker and Klümper [2005] J. Sirker and A. Klümper, Real-time dynamics at finite temperature by the density-matrix renormalization group: A path-integral approach, Phys. Rev. B 71, 241101 (2005).
  * Enss and Sirker [2012] T. Enss and J. Sirker, Light cone renormalization and quantum quenches in one-dimensional Hubbard models, New Journal of Physics 14, 023008 (2012).
  * Feiguin and White [2005] A. E. Feiguin and S. R. White, Finite-temperature density matrix renormalization using an enlarged Hilbert space, Phys. Rev. B 72, 220401 (2005).
  * Karrasch  _et al._ [2012] C. Karrasch, J. H. Bardarson, and J. E. Moore, Finite-temperature dynamical density matrix renormalization group and the Drude weight of spin-$1/2$ chains, Phys. Rev. Lett. 108, 227206 (2012).
  * Karrasch  _et al._ [2013] C. Karrasch, J. H. Bardarson, and J. E. Moore, Reducing the numerical effort of finite-temperature density matrix renormalization group calculations, New Journal of Physics 15, 083031 (2013).
  * Babenko  _et al._ [2021] C. Babenko, F. Göhmann, K. K. Kozlowski, J. Sirker, and J. Suzuki, Exact real-time longitudinal correlation functions of the massive XXZ chain, Phys. Rev. Lett. 126, 210602 (2021).
  * Göhmann  _et al._ [2022] F. Göhmann, K. K. Kozlowski, J. Sirker, and J. Suzuki, Spin conductivity of the XXZ chain in the antiferromagnetic massive regime, SciPost Phys. 12, 158 (2022).
  * White and Scalapino [1998] S. R. White and D. J. Scalapino, Density matrix renormalization group study of the striped phase in the 2D $\mathit{t}-\mathit{J}$ model, Phys. Rev. Lett. 80, 1272 (1998).
  * Zheng  _et al._ [2017] B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, Stripe order in the underdoped region of the two-dimensional Hubbard model, Science 358, 1155 (2017).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Carleo  _et al._ [2019] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Machine learning and the physical sciences, Rev. Mod. Phys. 91, 045002 (2019).
  * Carleo  _et al._ [2018] G. Carleo, Y. Nomura, and M. Imada, Constructing exact representations of quantum many-body systems with deep neural networks, Nat. Comm. 9, 5322 (2018).
  * Schmitt and Heyl [2020] M. Schmitt and M. Heyl, Quantum many-body dynamics in two dimensions with artificial neural networks, Phys. Rev. Lett. 125, 100503 (2020).
  * Cybenko [1989] G. Cybenko, Approximation by superpositions of a sigmoidal function., Math. Control Signal Systems 2, 303–314 (1989).
  * Hornik [1991] K. Hornik, Approximation capabilities of multilayer feedforward networks, Neural Networks 4, 251 (1991).
  * Voigtlaender [2023] F. Voigtlaender, The universal approximation theorem for complex-valued neural networks, Applied and Computational Harmonic Analysis 64, 33 (2023).
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Two soluble models of an antiferromagnetic chain, Annals of Physics 16, 407 (1961).
  * Szabó and Castelnovo [2020] A. Szabó and C. Castelnovo, Neural network wave functions and the sign problem, Phys. Rev. Research 2, 033075 (2020).
  * Nomura  _et al._ [2017] Y. Nomura, A. Darmawan, Y. Yamaji, and M. Imada, Restricted Boltzmann machine learning for solving strongly correlated quantum systems, Phys. Rev. B 96, 205152 (2017).
  * Nomura and Imada [2021] Y. Nomura and M. Imada, Dirac-type nodal spin liquid revealed by refined quantum many-body solver using neural-network wave function, correlation ratio, and level spectroscopy, Phys. Rev. X 11, 031034 (2021).
  * Hams and Raedt [2000] A. Hams and H. D. Raedt, Fast algorithm for finding the eigenvalue distribution of very large matrices, Phys. Rev. E 62, 4365 (2000).
  * Gemmer  _et al._ [2009] J. Gemmer, M. Michel, and G. Mahler, Typicality of observables and states, in  _Quantum Thermodynamics: Emergence of Thermodynamic Behavior Within Composite Quantum Systems_ (Springer Berlin Heidelberg, 2009) pp. 85–93.
  * Steinigeweg  _et al._ [2014] R. Steinigeweg, J. Gemmer, and W. Brenig, Spin-current autocorrelations from single pure-state propagation, Phys. Rev. Lett. 112, 120601 (2014).
  * Sugiura [2017] S. Sugiura,  _Formulation of Statistical Mechanics Based on Thermal Pure Quantum States_ (Springer Singapore, 2017).
  * Hendry  _et al._ [2022] D. Hendry, H. Chen, and A. Feiguin, Neural network representation for minimally entangled typical thermal states, Phys. Rev. B 106, 165111 (2022).
  * Hackl  _et al._ [2020] L. Hackl, T. Guaita, T. Shi, J. Haegeman, E. Demler, and J. I. Cirac, Geometry of variational methods: dynamics of closed quantum systems, SciPost Phys. 9, 48 (2020).
  * Sandvik, A. W. [2012] Sandvik, A. W., Quantum Monte Carlo methods at work for novel phases of matter (2012), http://physics.bu.edu/~sandvik/trieste12/index.html, Last accessed on 2023-10-13.
  * Sandvik [2010] A. W. Sandvik, Computational studies of quantum spin systems, AIP Conference Proceedings 1297, 135 (2010).
  * Zhang  _et al._ [2023] H. Zhang, R. J. Webber, M. Lindsey, T. C. Berkelbach, and J. Weare, Understanding and eliminating spurious modes in variational Monte Carlo using collective variables, Phys. Rev. Res. 5, 023101 (2023).
  * Deng  _et al._ [2017] D.-L. Deng, X. Li, and S. Das Sarma, Quantum entanglement in neural network states, Phys. Rev. X 7, 021021 (2017).
  * Lin and Pollmann [2022] S.-H. Lin and F. Pollmann, Scaling of neural-network quantum states for time evolution, physica status solidi (b) 259, 2100172 (2022).
  * White [2009] S. R. White, Minimally entangled typical quantum states at finite temperature, Phys. Rev. Lett. 102, 190601 (2009).
  * Short and Winter [2006] S. P. A. Short and A. Winter, Entanglement and the foundations of statistical mechanics, Nature Phys 2, 754–758 (2006).
  * Popescu  _et al._ [2005] S. Popescu, A. J. Short, and A. Winter, The foundations of statistical mechanics from entanglement: Individual states vs. averages, arXiv:quant-ph/0511225 (2005).
  * Milman and Schechtman [1986] V. Milman and G. Schechtman,  _Asymptotic Theory of Finite Dimensional Normed Spaces_ (Springer Berlin, Heidelberg, 1986).
  * Matoušek [2002] J. Matoušek,  _Lectures on Discrete Geometry_ (Springer New York, 2002).
  * Sugiura and Shimizu [2012] S. Sugiura and A. Shimizu, Thermal pure quantum states at finite temperature, Phys. Rev. Lett. 108, 240401 (2012).
  * Sugiura and Shimizu [2013] S. Sugiura and A. Shimizu, Canonical thermal pure quantum state, Phys. Rev. Lett. 111, 010401 (2013).
  * P.Reimann [2007] P.Reimann, Typicality for generalized microcanonical ensembles, Phys. Rev. Lett. 99, 160404 (2007).
  * Page [1993] D. N. Page, Average entropy of a subsystem, Phys. Rev. Lett. 71, 1291 (1993).
  * Bianchi  _et al._ [2022] E. Bianchi, L. Hackl, M. Kieburg, M. Rigol, and L. Vidmar, Volume-law entanglement entropy of typical pure quantum states, PRX Quantum 3, 030201 (2022).
  * Nakagawa  _et al._ [2018] Y. O. Nakagawa, M. Watanabe, H. Fujita, and S. Sugiura, Universality in volume-law entanglement of scrambled pure quantum states, Nature Communications 9, 1635 (2018).
  * Tahara and Imada [2008] D. Tahara and M. Imada, Variational Monte Carlo method combined with quantum-number projection and multi-variable optimization, J. Phys. Soc. Jpn. 77 (2008).
  * Bajdich  _et al._ [2006] M. Bajdich, L. Mitas, G. Drobný, and L. K. Wagner, Pfaffian pairing wave functions in electronic-structure quantum Monte Carlo simulations, Phys. Rev. Lett. 96, 130201 (2006).
  * Misawa  _et al._ [2019] T. Misawa, S. Morita, K. Yoshimi, M. Kawamura, Y. Motoyama, K. Ido, T. Ohgoe, M. Imada, and T. Kato, mvmc—open-source software for many-variable variational Monte Carlo method, Computer Physics Communications 235, 447 (2019).
  * Gao and Duan [2017] X. Gao and L.-M. Duan, Efficient representation of quantum many-body states with deep neural networks, Nature Communications 8, 662 (2017).