11institutetext: Department of Physics, University of Malta, Msida MSD 2080,
Malta 11email<EMAIL_ADDRESS>
# Variational Quantum Algorithms for Gibbs State Preparation

Mirko Consiglio  0000-0001-8730-0206

###### Abstract

Preparing the Gibbs state of an interacting quantum many-body system on noisy
intermediate-scale quantum (NISQ) devices is a crucial task for exploring the
thermodynamic properties in the quantum regime. It encompasses understanding
protocols such as thermalization and out-of-equilibrium thermodynamics, as
well as sampling from faithfully prepared Gibbs states could pave the way to
providing useful resources for quantum algorithms. VQA show the most promise
in efficiently preparing Gibbs states, however, there are many different
approaches that could be applied to effectively determine and prepare Gibbs
states on a NISQ computer. In this paper, we provide a concise overview of the
algorithms capable of preparing Gibbs states, including joint Hamiltonian
evolution of a system–environment coupling, quantum imaginary time evolution,
and modern variational quantum algorithms utilizing the Helmholtz free energy
as a cost function, among others. Furthermore, we perform a benchmark of one
of the latest variational Gibbs state preparation algorithms, developed by
Consiglio et al. [16], by applying it to the spin 1/2 one-dimensional $XY$
model.

###### Keywords:

Quantum Computing Quantum Algorithms Quantum Thermodynamics

PQC

    parametrized quantum circuit
VQA

    variational quantum algorithm
NISQ

    noisy intermediate-scale quantum
BFGS

    Broyden–Fletcher–Goldfarb–Shanno
SPSA

    simultaneous perturbation stochastic approximation
TFD

    thermofield double
QAOA

    quantum approximate optimization algorithm
QITE

    quantum imaginary time evolution
METTS

    minimally entangled typical thermal states
ADAPT

    adaptive derivative assembled problem-tailored
QBM

    quantum Boltzmann machine
CBM

    classical Boltzmann machine

##  1 Introduction

Gibbs states (also known as thermal states) can be used for quantum simulation
[12], quantum machine learning [36, 6], quantum optimization [63], and the
study of open quantum systems [52]. Moreover, semi-definite programming [9],
combinatorial optimization problems [63], and training quantum Boltzmann
machines [36], can be tackled by sampling from well-prepared Gibbs states.
Nevertheless, the preparation of an arbitrary quantum state on a quantum
computer, is a QMA-hard problem [78]. Preparing Gibbs states, specifically at
low temperatures, could be as hard as finding the ground-state of that
Hamiltonian [2].

From a physical point of view, A Gibbs state is the quantum state that is at
thermal equilibrium with the surrounding environment, in the canonical
ensemble. Let’s consider a Hamiltonian $\mathcal{H}$, describing $n$
interacting qubits. The Gibbs state at inverse temperature $\beta\equiv
1/(k_{\text{B}}T)$, where $k_{\text{B}}$ Boltzmann constant and $T$ is the
temperature, is defined as

$\rho(\beta,\mathcal{H})=\frac{e^{-\beta\mathcal{H}}}{\mathcal{Z}(\beta,\mathcal{H})},$
(1)

where the partition function $\mathcal{Z}(\beta,\mathcal{H})$ is

$\mathcal{Z}(\beta,\mathcal{H})=\Tr{e^{-\beta\mathcal{H}}}=\sum_{i=0}^{d-1}e^{-\beta
E_{i}}.$ (2)

Here the dimension $d=2^{n}$, while $\\{E_{i}\\}$ are the eigenenergies of
$\mathcal{H}$ (with $\\{\ket{E_{i}}\\}$ denoting the corresponding
eigenstates), i.e. $\mathcal{H}\ket{E_{i}}=E_{i}\ket{E_{i}}$.

While it may be a straightforward procedure to allow a system to naturally
thermalize with its environment, to obtain a Gibbs state, there is a great
significance in being able to obtain a Gibbs of any arbitrary Hamiltonian.
Since most physical systems are bound to specific types of Hamiltonian, a
quantum computer can — at least in principle — simulate any arbitrary quantum
system. Moreover, while a Gibbs distribution can be prepared on a classical
computer, or even have a density matrix describing the Gibbs state stored on
it, is not equivalent to having an actual Gibbs state on a quantum computer.
Since we can continue to evolve, probe and simulate the Gibbs state to study
quantum systems. As Richard Feynman aptly said [71], we require a quantum
system to be able to effectively simulate, another, quantum system.

The objective of this study is twofold: firstly, to offer a comprehensive
overview of Gibbs state preparation algorithms in Section 2, and to explore a
noteworthy application of Gibbs states in quantum Boltzmann machines.
Secondly, to present the recent variational Gibbs state preparation algorithm
by Consiglio et al. [16] in Section 3. In Section 4, we apply the VQA to
prepare Gibbs states of the $XY$ model, while also qualitatively investigating
the scalability of the algorithm. Finally in Section. 5, we draw our
conclusions.

##  2 Overview of Gibbs State Preparation Algorithms

The first computational technique one would usually resort to prepare a Gibbs
state would be exact diagonalisation [33], however, since the full energy
spectrum is needed, this method quickly becomes infeasible once the system
size grows. A different approach is based on thermal pure quantum states [64,
65], which significantly reduces computational overhead when compared with
exact diagonalisation. Nevertheless, since it requires preparing and evolving
random quantum states, it is still limited in terms of system size. On the
other hand, a quantum algorithm based on evaluating thermal pure shadows that
operates using a quantum computer was proposed [18].

Initially, the first algorithms for preparing Gibbs states on a quantum
computer were based on the idea of coupling the system to a register of
ancillary qubits, and letting the system and environment evolve under a joint
Hamiltonian, simulating the physical process of thermalization [70, 59, 46].
Davies [19, 20] showed that, in the limits of small time-steps or weak system-
bath interactions, the dynamics are described by a Lindblad master equation —
the Davies generator — that thermalizes the system in the long-time limit.
Nevertheless, implementing time-evolution on a NISQ devices is currently
impractical, since it requires long coherence times, and a significant number
of ancilla qubits and/or precise qubit reset capabilities. One could resort to
implementing the Davies generator directly on the quantum computer, however
this would once again require applying quantum phase estimation [55] or other
quantum subroutines related to the quantum Fourier transform [10].
Furthermore, Ref. [5] and Ref. [35] studied thermalization and Davies
generators from the perspective of mathematical physics, showing rapid
convergence to the Gibbs state and Gibbs distribution, respectively, under
certain conditions. Ref. [8] devised a method for preparing the Gibbs state of
a local Hamiltonian, assuming that the state has an exponential decay of
correlations, and that it satisfies the approximate Markov condition uniformly
on non-contractible regions. Ref. [83] also developed a quantum algorithm,
based on a local quantum Markov process that can be used to sample from the
Gibbs distribution.

Temme et al. [69] proposed the quantum Metropolis algorithm, that is capable
of preparing Gibbs states through random walks, which is inspired from the
Metropolis–Hastings algorithm [47, 28]. A quadratic speedup is achieved over
the direct implementation of the Metropolis–Hastings algorithm by using
Grover’s algorithm or related quantum algorithmic techniques [52, 11]. Using
modern quantum algorithmic approaches, further polynomial speedups could be
achieved, such as using dimension reduction to sample from Gibbs states [7],
which requires the use of a precise quantum phase estimation subroutine, among
others [82, 51, 14, 79]. A similar approach is the one in Ref. [56], which
requires the use of a circuit very similar to that used in Shor’s algorithm
[60], which can only be feasibly implemented on fault-tolerant quantum
devices. Ref. [22] showed how to adapt “Coupling from the Past”-algorithms
proposed by Ref. [54], to sample from the exact Gibbs distribution, rather
than an approximation, using the quantum Metropolis algorithm. Nevertheless,
compared with the local updates of the classical Metropolis algorithm, all of
these algorithms require the implementation of large, global quantum circuits
across the whole system, making them infeasible for NISQ devices. A recent
review of current quantum Gibbs-sampling algorithms can be found in Ref. [10].

A promising approach for NISQ devices is VQAs, where a hybrid
quantum–classical approach of minimizing an objective function, using a
parametrized quantum circuit (PQC) as a variational ansatz, leads to the
preparation of the Gibbs state. Some approaches utilize a physically-inspired
objective function, the Helmholtz free energy, such in Refs. [80, 43, 13, 21,
58, 16]. Others employ an engineered cost function [53, 57], while others use
quantum approximate optimization algorithm (QAOA)-based approaches [80, 84].
Alternative variational approaches consist of using: truncated Taylor series
to evaluate an approximation of the free energy [76]; adaptive derivative
assembled problem-tailored (ADAPT)-VQA applied to a similar cost function to
the free energy [77]; and a VQA based on McLachlan’s variational principle to
initialize and evolve a thermofield double (TFD) state [43, 80, 84, 53, 57,
37]. A TFD state is defined as a pure state written in basis vectors
$\ket{i}_{A}$ of system $A$, and $\ket{\tilde{i}}_{B}$ of system $B$:

$\ket{\textrm{TFD}(\beta,\mathcal{H})}=\sum_{i=0}^{d-1}\sqrt{\frac{e^{-\beta
E_{i}}}{\mathcal{Z}(\beta,\mathcal{H})}}\ket{i}_{A}\otimes\ket{\tilde{i}}_{B}.$
(3)

Tracing out system $A$ ($B$) yields the Gibbs state on system $B$ ($A$). TFD
states also find applications in black hole theory [32, 41] and teleportation
through traversable wormholes [24, 42, 25].

Alternative algorithms prepare thermal states through quantum imaginary time
evolution (QITE), which typically require, either starting from a maximally
mixed state [45, 23, 85], or a maximally entangled state [81], and carrying
out imaginary time evolution of time $\tau$, resulting in preparing the Gibbs
state at inverse temperature $\beta$, that is directly proportional to $\tau$.
Refs. [48, 66, 26] utilize minimally entangled typical thermal states (METTS),
in combination with QITE or measurement-based variational QITE to compute
thermal averages and correlation functions of Hamiltonians. Ref. [68]
developed a linear scaling QITE algorithm, called Fast QITE, that can also
compute thermal averages. Ref. [75] also employed variational QITE, that is
inspired from the measurement-based approach [48], however, it is adapted to
be ansatz-based QITE. Ref. [61] utilizes random quantum circuits with
intermediate measurements to impose QITE, while Ref. [62] proposed a
fragmented QITE algorithm, for sampling from the Gibbs distribution.

Variational ansätze based on multi-scale entanglement renormalization [58] and
product spectrum ansatz [43] have also been proposed in order to prepare Gibbs
states. Other VQA approaches rely on variational quantum simulation [74, 27]
and variational autoregressive networks [39]. The hybrid quantum–classical
algorithms proposed in Ref. [40] differ from VQA, in that they compute
microcanonical and canonical properties of many-body systems, using: filtering
operators, similar to quantum phase estimation; and Monte Carlo simulations.
Ref. [15] employs Green’s functions to sample from the Gibbs distribution.
Finally, Ref. [29] proposed quantum-assisted simulation to prepare thermal
states, which, contrary to VQAs, does not require a hybrid quantum–classical
feedback loop.

###  2.1 Quantum Boltzmann Machines

One pertinent application of Gibbs states is in quantum Boltzmann machines
[3], which are a type of quantum machine learning model that are based on the
principles of classical Boltzmann machines [1], but incorporate quantum
effects such as entanglement and superposition. CBMs are variants of neural
networks embedded in an undirected graph, where the weights and biases of the
network represent the information encoded within it. Typically, the network
consists of two groups of nodes: the visible nodes, that determine the input
and output of the network; and the hidden nodes; which act as latent variables
[85]. The purpose of training a Boltzmann machine is to learn the sets of
weights and biases such that the resulting network is capable of approximating
a target probability distribution, whereby the model learns by feeding in
training data. The Boltzmann machine can then be used for both discriminative
and generative learning [3].

A Boltzmann machine is represented by an undirected graph, where $\bm{V}$ is
the set of vertices (nodes) and $\bm{E}$ is the set of edges. The states of
the CBM are defined by $\bm{Z}=\\{\bm{U},\bm{H}\\}$, where $\bm{U}$ represents
the state of the visible nodes, and $\bm{H}$ represents the state of the
hidden nodes. This defines an energy function

$E(\bm{Z})=-\sum_{i\in\bm{V}}b_{i}z_{i}-\sum_{\\{i,j\\}\in\bm{E}}w_{ij}z_{i}z_{j},$
(4)

with $w_{ij}\in\bm{w},b_{i}\in\bm{b}$, denoting the weights and biases
(parameters) of the model, respectively. $z_{i}$ is a binary unit, and to
remain consistent with quantum mechanics, we define $z_{i}\in\\{-1,1\\}$.
Thus, the probability to observe a particular configuration $\bm{U}$ of
visible nodes is defined as

$p_{\bm{U}}=\frac{1}{\mathcal{Z}(\beta,\bm{Z})}\sum_{\bm{H}}e^{-\beta
E(\bm{Z})},$ (5)

where in this case,

$\mathcal{Z}(\beta,\bm{Z})=\sum_{\bm{Z}}e^{-\beta E(\bm{Z})}.$ (6)

The Boltzmann distribution summed over the hidden variables in Eq. (5) is
called the marginal distribution. The goal of a CBM is to determine the
weights and biases of (4), such that $p_{\bm{U}}$ approximates well
$p^{\text{data}}_{\bm{U}}$ defined by the training data. Typically, the cost
function employed in the optimization procedure of a CBM is the Kullback-
Leibler divergence, defined as

$\mathcal{S}(p^{\text{data}}_{\bm{U}}\|p_{\bm{U}})=\sum_{\bm{U}}p^{\text{data}}_{\bm{U}}\log
p^{\text{data}}_{\bm{U}}-\sum_{\bm{U}}p^{\text{data}}_{\bm{U}}\log
p_{\bm{U}}.$ (7)

The main difference between a CBM and a QBM, is that the latter employs nodes
that are defined by Pauli matrices rather than binary units, and therefore
constructs a network defined by a parameterized Hamiltonian:

$\mathcal{H}(\bm{b},\bm{w})=-\sum_{i\in\bm{V}}b_{i}\sigma^{z}_{i}-\sum_{\\{i,j\\}\in\bm{E}}w_{ij}\sigma^{z}_{i}\sigma^{z}_{j}.$
(8)

In this case, the QBM is defined as the Gibbs state

$\rho(\beta,\mathcal{H})=\frac{e^{-\beta\mathcal{H}(\bm{b},\bm{w})}}{\mathcal{Z}(\beta,\mathcal{H})},$
(9)

where

$\mathcal{Z}(\beta,\mathcal{H})=\Tr{e^{-\beta\mathcal{H}(\bm{b},\bm{w})}}.$
(10)

The cost function is thus the generalized Kullback-Leibler divergence for
density matrices, which corresponds to the quantum relative entropy:

$\mathcal{S}(\sigma\|\rho(\beta,\mathcal{H}))=\Tr{\sigma\ln\sigma}-\Tr{\sigma\ln\rho(\beta,\mathcal{H})},$
(11)

where $\sigma$ is the target density matrix representing the data embedded
into a mixed state.

The goal of QBM training is to find a set of weights and biases of the
Hamiltonian, such that the Gibbs state approximates the target density matrix.
There are various methods of tackling the training of QBMs [23, 85, 18, 34,
31], and an analogous problem to QBM training is Hamiltonian learning [4].

##  3 Variational Gibbs State Preparation

In this section we will describe how one can prepare a Gibbs state using a VQA
with the free energy as a cost function, specifically discussing the algorithm
presented in Ref. [16]. Section 3.1 and 3.2 are intended to serve as a
qualitative description of the estimation of the von Neumann entropy and the
free energy, respectively. Consequently, in Section 3.3, a quantitative
description of the framework of the VQA is carried out, mathematically linking
the von Neumann entropy and energy expectation, with the estimation of the
free energy.

Suppose one chooses a Hamiltonian $\mathcal{H}$ and inverse temperature
$\beta$. For a general state $\rho$, one can define a generalized Helmholtz
free energy as

$\mathcal{F}(\rho)=\Tr{\mathcal{H}\rho}-\beta^{-1}\mathcal{S}(\rho),$ (12)

where the von Neumann entropy $\mathcal{S}(\rho)$ can be expressed in terms of
the eigenvalues, $p_{i}$, of $\rho$,

$\mathcal{S}(\rho)=-\sum_{i=0}^{d-1}p_{i}\ln p_{i}.$ (13)

Since the Gibbs state is the unique state that minimizes the free energy of
$\mathcal{H}$ [44], a variational form can be put forward that takes Eq. (12)
as an objective function, such that

$\rho(\beta,\mathcal{H})=\underset{\rho}{\arg\min}~{}\mathcal{F}(\rho).$ (14)

In this case, $p_{i}=\exp\left(-\beta
E_{i}\right)/\mathcal{Z}(\beta,\mathcal{H})$ is the probability of getting the
eigenstate $\ket{E_{i}}$ from the ensemble $\rho(\beta,\mathcal{H})$.

###  3.1 Computing the von Neumann Entropy

The difficulty in measuring the von Neumann entropy, defined by Eq. (13), of a
quantum state on a NISQ device is typically the challenging part of
variational Gibbs state preparation algorithms, since $\mathcal{S}(\rho)$ is
not an observable. To exactly compute the von Neumann entropy on a quantum
device, one would have to perform full state tomography on the Gibbs state,
which has a time complexity of $\mathcal{O}(3^{n})$ and a space complexity of
$\mathcal{O}(4^{n})$. As a consequence, different works have employed various
techniques to compute an approximation for the von Neumann entropy, such as
using: a Fourier series approximation to the von Neumann entropy [13]; a
thermal multi-scale entanglement renormalization ansatz [58]; sums of Rényi
entropies [21]; among others. On the other hand, Ref. [16] devised a method of
computing the von Neumann exactly, via post-processing of measurement results
acting on ancillary register.

When preparing an $n$-qubit state, the unitary gates used in quantum computers
ensure that the final quantum state of the entire register, starting from the
input state $\ket{0}^{\otimes n}$, remains pure. As a result, in order to
prepare an $n$-qubit Gibbs state on the system register, an $m\leq n$-qubit
ancillary register is required. For example, in the case of the infinite-
temperature Gibbs state, which is a fully mixed state, $m=n$ qubits are needed
in the ancillary register to achieve maximal von Neumann entropy. To
accurately evaluate the von Neumann entropy without any approximations, the
complete Boltzmann distribution is prepared on the ancillary register, thus,
$m=n$ is set regardless of the temperature.

###  3.2 Computing the Free Energy

Following the prescription of Ref. [16], we shall denote the ancillary
register as $A$, while the preparation of the Gibbs state will be carried out
on the system register $S$. The purpose of the VQA is to effectively create
the Boltzmann distribution on $A$, which is then imposed on $S$, via
intermediary CNOT gates, to generate a diagonal mixed state in the
computational basis. In the ancillary register we can choose a unitary ansatz
capable of preparing such a probability distribution. Thus, the ancillary
qubits are responsible for classically mixing in the probabilities of the
thermal state, while also being able to access these probabilities via
measurements in the computational basis. On the other hand, the system
register will host the preparation of the Gibbs state, as well as the
measurement of the expectation value of our desired Hamiltonian during the
optimization routine.

The specific design of the PQC allows classical post-processing of simple
measurement results, carried out on ancillary qubits in the computational
basis, to determine the von Neumann entropy. A diagrammatic representation of
the structure of the PQC is shown in Fig. 1. Note that while the PQC of the
algorithm has to have a particular structure — a unitary acting on the
ancillae and a unitary acting on the system, connected by intermediary CNOT
gates — it is not dependent on the choice of Hamiltonian $\mathcal{H}$,
inverse temperature $\beta$, or the variational ansätze, $U_{A}$ and $U_{S}$,
employed within. This is in addition to enjoying a sub-exponential scaling in
the number of shots needed to precisely compute the Boltzmann probabilities. A
qualitative analysis of a power-law scaling is presented in Section 4.2.

###  3.3 Modular Structure of the PQC

The PQC, as shown in Fig. 1 for the VQA, is composed of a unitary gate $U_{A}$
acting on the ancillary qubits, and a unitary gate $U_{S}$ acting on the
system qubits, with CNOT gates in between. Note that the circuit notation we
are using here means that there are $n$ qubits for both the system and the
ancillae, as well as $n$ CNOT gates that act in parallel, and are denoted as

$\textsc{CNOT}_{AS}\equiv\bigotimes\limits_{i=0}^{n-1}\textsc{CNOT}_{A_{i}S_{i}}.$
(15)

The parametrized unitary $U_{A}$ acting on the ancillae, followed by CNOT
gates between the ancillary and system qubits, is responsible for preparing a
probability distribution on the system. The parametrized unitary $U_{S}$ is
then applied on the system qubits to transform the computational basis states
into the eigenstates of the Hamiltonian.

Figure 1: PQC for Gibbs state preparation, with systems $A$ and $S$ each
carrying $n$ qubits. CNOT gates act between each qubit $A_{i}$ and
corresponding $S_{i}$.

Throughout the paper we will describe unitary operations and density matrices
in the computational basis: $\ket{0}\equiv\left(1~{}0\right)^{\top}$ and
$\ket{1}\equiv\left(0~{}1\right)^{\top}$, which is spanned by one qubit (and
products of the basis states in the case of multiple qubits). We will also
denote the $n$-fold tensor product of a state $\ket{\psi}$, of an $n$-qubit
register $R$, as $\ket{\psi}_{R}^{\otimes
n}\equiv\bigotimes_{i=0}^{n-1}\ket{\psi}_{i}$, where $i\in R$. Denote a
general unitary gate of dimension $d=2^{n}$, as $U_{A}=(u_{i,j})_{0\leq
i,j\leq d-1}$. Starting with the initial state of the $2n$-qubit register,
$\ket{0}_{AS}^{\otimes 2n}$, we apply the unitary gate $U_{A}$ on the ancillae
to get a quantum state $\ket{\psi}_{A}$, such that

$(U_{A}\otimes I_{S})\ket{0}_{AS}^{\otimes
2n}=\ket{\psi}_{A}\otimes\ket{0}_{S}^{\otimes n},$ (16)

where $\ket{\psi}_{A}=\sum_{i=0}^{d-1}u_{i,0}\ket{i}_{A}$ and $I_{S}$ is the
identity acting on the system. Since Eq. (16) is applied to the all-zero state
in the ancillary register, then the operation serves to extract the first
column of the unitary operator $U_{A}$. The next step is to prepare a
classical probability mixture on the system qubits, which can be done by
applying CNOT gates between each ancilla and system qubit. This results in a
state

$\displaystyle\textsc{CNOT}_{AS}\left(\ket{\psi}_{A}\otimes\ket{0}_{S}^{\otimes
n}\right)=\sum_{i=0}^{d-1}u_{i,0}\ket{i}_{A}\otimes\ket{i}_{S}.$ (17)

By then tracing out the ancillary qubits, we arrive at

$\displaystyle\text{Tr}_{A}\left\\{\left(\sum_{i=0}^{d-1}u_{i,0}\ket{i}_{A}\otimes\ket{i}_{S}\right)\left(\sum_{j=0}^{d-1}u_{j,0}^{*}\bra{j}_{A}\otimes\bra{j}_{S}\right)\right\\}=$
$\displaystyle\sum_{i,j=0}^{d-1}u_{i,0}u_{j,0}^{*}\innerproduct{i}{j}\outerproduct{i}{j}_{S}$
$\displaystyle=\ $
$\displaystyle\sum_{i=0}^{d-1}|u_{i,0}|^{2}\outerproduct{i}{i}_{S},$ (18)

ending up with a diagonal mixed state on the system, with probabilities given
directly by the absolute square of the entries of the first column of $U_{A}$,
that is, $p_{i}=|u_{i,0}|^{2}$. If the system qubits were traced out instead,
we would end up with the same diagonal mixed state,

$\displaystyle\text{Tr}_{S}\left\\{\left(\sum_{i=0}^{d-1}u_{i,0}\ket{i}_{A}\otimes\ket{i}_{S}\right)\left(\sum_{j=0}^{d-1}u_{j,0}^{*}\bra{j}_{A}\otimes\bra{j}_{S}\right)\right\\}=$
$\displaystyle\sum_{i,j=0}^{d-1}u_{i,0}u_{j,0}^{*}\innerproduct{i}{j}\outerproduct{i}{j}_{A}$
$\displaystyle=\ $
$\displaystyle\sum_{i=0}^{d-1}|u_{i,0}|^{2}\outerproduct{i}{i}_{A}.$ (19)

This implies that by measuring in the computational basis of the ancillary
qubits, we can determine the probabilities $p_{i}$, which can then be post-
processed to determine the von Neumann entropy $\mathcal{S}$ of the state
$\rho$ via Eq. (13) (since the entropy of $A$ is the same as that of $S$). As
a result, since $U_{A}$ only serves to create a probability distribution from
the entries of the first column, we can do away with a parametrized orthogonal
(real unitary) operator, thus requiring less gates and parameters for the
ancillary ansatz.

The unitary gate $U_{S}$ then serves to transform the computational basis
states of the system qubits to the eigenstates of the Gibbs state, once it is
optimized by the VQA, such that

$\displaystyle\rho=U_{S}\left(\sum_{i=0}^{d-1}|u_{i,0}|^{2}\outerproduct{i}{i}_{S}\right)U_{S}^{\dagger}=\sum_{i=0}^{d-1}p_{i}\outerproduct{\psi_{i}}{\psi_{i}},$
(20)

where the expectation value $\Tr{\mathcal{H}\rho}$ of the Hamiltonian can be
measured. Ideally, at the end of the optimization procedure,
$p_{i}=\exp\left(-\beta E_{i}\right)/\mathcal{Z}(\beta,\mathcal{H})$ and
$\ket{\psi_{i}}=\ket{E_{i}}$, so that we get

$\rho(\beta,\mathcal{H})=\sum_{i=0}^{d-1}\frac{e^{-\beta
E_{i}}}{\mathcal{Z}(\beta,\mathcal{H})}\outerproduct{E_{i}}{E_{i}}.$ (21)

The VQA therefore avoids the entire difficulty of measuring the von Neumann
entropy of a mixed state on a quantum computer, and instead transfers the task
of post-processing measurement results to the classical computer, which is
much more tractable.

###  3.4 Objective Function

Finally, we can define the objective function of our VQA to minimize the free
energy (12), via our constructed PQC, to obtain the Gibbs state

$\displaystyle\rho(\beta,\mathcal{H})$
$\displaystyle=\underset{\bm{\theta},\bm{\varphi}}{\arg\min}\
\mathcal{F}\left(\rho\left(\bm{\theta},\bm{\varphi}\right)\right)$
$\displaystyle=\underset{\bm{\theta},\bm{\varphi}}{\arg\min}\left(\Tr{\mathcal{H}\rho_{S}(\bm{\theta},\bm{\varphi})}-\beta^{-1}\mathcal{S}\left(\rho_{A}(\bm{\theta})\right)\right).$
(22)

It is noteworthy to mention that while the energy expectation depends on both
sets of angles $\bm{\theta}$ (as $U_{A}$ is responsible for parameterizing the
Boltzmann distribution) and $\bm{\varphi}$ (as $U_{S}$ is responsible for
parameterizing the eigenstates of the Gibbs state), the calculation of the von
Neumann entropy only depends on $\bm{\theta}$.

Furthermore, once we obtain the optimal parameters $\bm{\theta^{*}}$,
$\bm{\varphi^{*}}$, preparing the Gibbs state $\rho(\beta,\mathcal{H})$ on the
system qubits $S$, one can place the same unitary $U_{S}$ with optimal
parameters $\bm{\varphi^{*}}$ on the ancillary qubits to prepare the TFD state
on the quantum computer, as shown in Fig. 2. Notice that this is equivalent to
Eq. (17) followed by applying the optimized $U_{S}$ on both registers. Tracing
out either the ancilla or system register yields the same Gibbs state on the
other register.

Figure 2: Optimal PQC for TFD state preparation, with systems $A$ and $S$ each
carrying $n$ qubits. CNOT gates act between each qubit $A_{i}$ and
corresponding $S_{i}$.

##  4 Preparing Gibbs states of the XY model

In this Section we assess the performance of the VQA for Gibbs state
preparation of an $XY$ model. The $XY$ model [38] is defined as

$\mathcal{H}=-\sum_{i=1}^{n}\left(\frac{1+\gamma}{2}\sigma_{i}^{x}\sigma_{i+1}^{x}+\frac{1-\gamma}{2}\sigma_{i}^{y}\sigma_{i+1}^{y}\right)-h\sum_{i=1}^{n}\sigma_{i}^{z}.$
(23)

Here we only report one relevant property for implementing a problem-inspired
ansatz for $U_{S}$. The Hamiltonian in Eq. (23) commutes with the parity
operator $\mathcal{P}=\prod_{i=0}^{n-1}\sigma^{z}_{i}$. As a consequence, the
eigenstates of $\mathcal{H}$ have definite parity, and so will the eigenstates
of $\rho_{\beta}$.

To assess the performance of the VQA, we utilize the Uhlmann-Josza fidelity as
a figure of merit [72], defined as
$F\left(\rho,\sigma\right)=\left(\Tr{\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}}}\right)^{2}$.
This fidelity measure quantifies the “closeness” between the prepared state
and the target Gibbs state, making it a commonly-employed metric for
distinguishability. However, alternative measures exist, each with its own
interpretation. For instance, the trace distance [49] can be used, which
guarantees that if its value between two states is bounded by $\epsilon$, the
expectation values computed on the effectively prepared state will differ from
those of the true Gibbs state by at most an amount proportional to
$\epsilon$[30]. Another option is the relative entropy[49], which
characterizes the distinguishability between the two states as the surprise
that arises when an event occurs that is not possible with the true Gibbs
state [73].

We use an alternating, entangling brick-wall PQC for the unitary $U_{A}$,
composed of parametrized $R_{Y}(\theta_{i})$ gates, and CNOTs as the
entangling gates. This ansatz is hardware efficient and is sufficient to
produce real amplitudes for preparing the probability distribution. Note that
we require the use of entangling gates [16], as otherwise we will not be able
to prepare any arbitrary probability distribution, including the Boltzmann
distribution of the $XY$ model.

For the unitary $U_{S}$, We employ a brick-wall structure solely using parity-
preserving gates — $R_{XY}(\varphi_{i})$ followed by $R_{YX}(\varphi_{j})$
gates. If we combine the two gates, which we denote as
$R_{P}(\varphi_{i},\varphi_{j})$, we get

$\displaystyle R_{P}(\varphi_{i},\varphi_{j})$
$\displaystyle=R_{YX}(\varphi_{j})\cdot R_{XY}(\varphi_{i})$
$\displaystyle=\left(\begin{array}[]{cccc}\cos\left(\frac{\varphi_{i}+\varphi_{j}}{2}\right)&0&0&\sin\left(\frac{\varphi_{i}+\varphi_{j}}{2}\right)\\\
0&\cos\left(\frac{\varphi_{i}-\varphi_{j}}{2}\right)&-\sin\left(\frac{\varphi_{i}-\varphi_{j}}{2}\right)&0\\\
0&\sin\left(\frac{\varphi_{i}-\varphi_{j}}{2}\right)&\cos\left(\frac{\varphi_{i}-\varphi_{j}}{2}\right)&0\\\
-\sin\left(\frac{\varphi_{i}+\varphi_{j}}{2}\right)&0&0&\cos\left(\frac{\varphi_{i}+\varphi_{j}}{2}\right)\\\
\end{array}\right),$ (28)

which can be decomposed into two CNOT gates six $\sqrt{X}$ gates, and ten
$R_{Z}$ gates. The decomposed unitary is shown in Fig. 3. One layer of the
unitary acting on the system qubits consists of a brick-wall structure,
composed of an even-odd sublayer of $R_{P}$ gates, followed by an odd-even
sublayer of $R_{P}$ gates. Table 1 shows the scaling of the VQA assuming a
closed ladder connectivity, for $n>2$. An example of a PQC, split into a four-
qubit ancilla register, and a four-qubit system register, is shown in Fig. 4.

Figure 3: Decomposed $R_{P}$ gate in Eq. (28). Table 1: Scaling of the VQA assuming a closed ladder connectivity, for $n>2$, where $l_{A}$ and $l_{S}$ are the number of ancilla ansatz and system ansatz layers, respectively, and $P$ is 2 when $n$ is even and 3 when $n$ is odd. Circuit depth is calculated on the depth of CNOT gates. # of parameters | $n(l_{A}+1)+2nl_{S}$ | $\mathcal{O}(n(l_{A}+l_{S}))$
---|---|---
# of CNOT gates | $nl_{A}+2nl_{S}+n$ | $\mathcal{O}(n(l_{A}+l_{S}))$
# of $\sqrt{X}$ gates | $2n(l_{A}+1)+6nl_{S}$ | $\mathcal{O}(n(l_{A}+l_{S}))$
Circuit depth | $Pl_{A}+2Pl_{S}+1$ | $\mathcal{O}(l_{A}+l_{S})$
Figure 4: Example of an eight-qubit PQC, consisting of three ancilla layers
acting on a four-qubit register, and three system layers acting on another
four-qubit register. Each $R_{Y}$ gate is parametrized with one parameter
$\theta_{i}$, while each $R_{P}$ gate has two parameters $\varphi_{i}$ and
$\varphi_{j}$. The $R_{P}$ gate is defined in Eq. (28) with its decomposition
shown in Fig. 3. Note that the intermediary CNOT gates, as well as the CNOT
and $R_{P}$ gates acting on qubits two and three, and on qubits one and four
of the ancilla and system, respectively, can be carried out in parallel,
respectively.

###  4.1 Statevector Results

Fig. 5 shows the fidelity of the generated mixed state when compared with the
exact Gibbs state of the $XY$ model with $h=0.5$, and $\gamma$ ranging from
0.1 to 0.9 in steps of 0.1, across a range of temperatures for system sizes
between two to seven qubits. The VQA was carried out using statevector
simulations with the Broyden–Fletcher–Goldfarb–Shanno (BFGS) optimizer [50].
We used $n-1$ layers for both the ancilla ansatz and for the system ansatz,
with the scaling highlighted in Table 2. The number of layers was
heuristically chosen to satisfy, at most, a polynomial scaling in quantum
resources (but linear in depth), while achieving a fidelity higher than 99%.
Furthermore, in order to alleviate the issue of getting stuck in local minima,
the optimizer was embedded in a Monte Carlo framework, that is, taking ten
random initial positions and carrying out a local optimization from each
position — which we call a ‘run’ — and finally taking the global minimum to be
the minimum over all runs.

Figure 5: Fidelity $F$, of the obtained state via statevector simulations (using BFGS) with the exact Gibbs state, vs inverse temperature $\beta$, for two to seven qubits of the $XY$ model with $\gamma$ between 0.1 and 0.9, and $h=0.5$. A total of ten runs are made for each point, with the optimal state taken to be the one that maximizes the fidelity. Table 2: Scaling of the VQA assuming a closed ladder connectivity, for $n>2$, with $l_{A}=n-1$, and $l_{S}=n-1$, and $P$ is 2 when $n$ is even and 3 when $n$ is odd. Circuit depth is calculated on the depth of CNOT gates. # of parameters | $n(3n-2)$ | $\mathcal{O}(n^{2})$
---|---|---
# of CNOT gates | $3n^{2}-2n$ | $\mathcal{O}(n^{2})$
# of $\sqrt{X}$ gates | $2n(4n-3)$ | $\mathcal{O}(n^{2})$
Circuit depth | $3P(n-1)+1$ | $\mathcal{O}(n)$

A total of ten runs of BFGS per $\beta$ were carried out to verify the
reachability of the PQC, with Fig. 5 showcasing the maximal fidelity achieved
for each $\beta$ out of all runs. The results show that, indeed, our VQA, is
able to reach a very high fidelity $F>99\%$ for up to seven-qubit Gibbs states
of the $XY$ model. In the case of the extremal points, that is
$\beta\rightarrow 0$ and $\beta\rightarrow\infty$, the fidelity reaches unity,
for all investigated system sizes. For $\beta\rightarrow 0$, the problem
reduces to maximizing the von Neumann entropy, where the problem becomes
independent of the energy, and therefore, independent of the prepared
eigenstates by $U_{S}$. For $\beta\rightarrow\infty$, the problem reduces to
minimizing the energy, which is equivalent to finding the ground-state of the
Hamiltonian, and so, independent of $U_{A}$. This is in contrast with
intermediary temperatures $\beta\sim 1$, where the fidelity decreases with the
number of qubits. This could be attributed to the fact that at intermediate
temperatures, a substantial amount of Boltzmann probabilities are necessary
for evaluating the von Neumann entropy with a high precision, and similarly, a
large number of eigenstates contribute to preparing the Gibbs state, resulting
in a dip in the fidelity at those temperatures.

###  4.2 Error Analysis

In this Section, we perform an analysis of the error scaling with the number
of shots, in faithfully identifying the probability distribution,
$\bm{p}=\\{p_{0},p_{1},\dots,p_{d-1}\\}$, where $d=2^{n}$, prepared by
$U_{A}$. The outcome of one shot of a quantum circuit can be described by a
multinomial distribution, where $p_{i}$ is the probability of observing a bit
string $i$. The expected value of a multinomially distributed random bit
string $i$ is

$\mu=\mathbb{E}[i]=N_{s}p_{i},$ (29)

with the variance being

$\sigma^{2}=\mathrm{Var}[i]=N_{s}p_{i}(1-p_{i}),$ (30)

where $N_{s}$ is the number of shots. A quantity that can describe the
precision of measuring the bit string $i$, $p_{i}$ times, is the coefficient
of variation, or relative standard deviation,

$c_{\text{v}}=\frac{\sigma}{\mu}=\sqrt{\frac{1-p_{i}}{N_{s}p_{i}}}.$ (31)

Given that $p_{i}=\exp(-\beta E_{i})/\mathcal{Z}(\beta,\mathcal{H})$ for the
Boltzmann distribution, we get

$c_{\text{v}}=\sqrt{\frac{1-\frac{e^{-\beta
E_{i}}}{\mathcal{Z}(\beta,\mathcal{H})}}{N_{\text{s}}\frac{e^{-\beta
E_{i}}}{\mathcal{Z}(\beta,\mathcal{H})}}}=\sqrt{\frac{\mathcal{Z}(\beta,\mathcal{H})-e^{-\beta
E_{i}}}{N_{\text{s}}e^{-\beta
E_{i}}}}=\sqrt{\frac{1}{N_{\text{s}}}\left(\frac{\mathcal{Z}(\beta,\mathcal{H})}{e^{-\beta
E_{i}}}-1\right)}.$ (32)

From Eq. (32), if $\beta\rightarrow 0$, then
$\mathcal{Z}(\beta,\mathcal{H})/\exp(-\beta E_{i})\rightarrow d~{}\forall
i\in[d]$, which implies

$c_{\text{v}}(\beta\rightarrow 0)=\sqrt{\frac{d-1}{N_{\text{s}}}},$ (33)

hence, exhibiting an exponential scaling in the number of shots needed for
preparing the flat distribution. On the other hand, if
$\beta\rightarrow\infty$, then $\mathcal{Z}(\beta,\mathcal{H})/\exp(-\beta
E_{0})\rightarrow 1$, while all other probabilities tend to zero, and the
algorithm reduces to finding the ground-state of the Hamiltonian. As a result,

$c_{\text{v}}(\beta\rightarrow\infty)=0.$ (34)

For intermediary temperatures, we carry out a qualitative analysis. Fig. 6
reports the exponent $\alpha_{i}$, of a polynomial fit of the normalized
coefficient of variation, defined as

$c_{\text{v}}\sqrt{N_{\text{s}}}=Cn^{\alpha_{i}},$ (35)

for $n$ between 8 and 20, at different $\beta$, where the first 51 $p_{i}$ of
the $XY$ model with $h=0.5$, and $\gamma$ between 0.1 and 0.9 are considered.
We observe that as $\beta\rightarrow\infty$, the exponent $\alpha_{0}$
approaches zero only for the ground-state occupation probability $p_{0}$, in
accordance with Eq. (34). Simultaneously, as the occupation probability
$p_{i}$ of an excited state decreases to negligible values, the corresponding
exponent $\alpha_{i}$ becomes increasingly negative. However, for any finite
temperature, we find that a polynomial fit, with $\alpha_{i}\lesssim 6$,
provides a good approximation of the normalized coefficient of variation. This
suggests that achieving a faithful reconstruction of the probability
distribution, for the investigated system sizes of the $XY$ model in NISQ
algorithms, requires a sub-exponential scaling in the number of shots.

Figure 6: The exponent $\alpha_{i}$ in Eq. (35) is plotted vs $\beta$ for the
first 51 $p_{i}$ of the $XY$ model with $\gamma$ between 0.1 and 0.9, and
$h=0.5$, where each point is the polynomial fit for $n$ between 8 and 20.

##  5 Conclusion

In this manuscript, we provided a concise overview of various Gibbs state
preparation and Gibbs-sampling algorithms. However, our primary focus was on
investigating VQAs for Gibbs state preparation. To accomplish this, we
conducted a benchmark study using one of the latest variational Gibbs state
preparation algorithms on the $XY$ model, considering a wide range of
temperatures and $\gamma$ coefficients. Following the approach proposed by
Consiglio et al. [16], we leveraged the unique property of the Gibbs state as
the state that minimizes the Helmholtz free energy, which serves as an
appropriate objective function for the VQA. Through extensive statevector
simulations, we achieved fidelities of $F>99\%$ for system sizes up to seven
qubits.

Furthermore, we performed a qualitative analysis of the scalability of the VQA
for implementation on a NISQ device. We found that the number of gates,
iterations, and shots scale at most polynomially with the number of qubits in
the $XY$ model. The algorithm’s scalability not only makes it suitable for
near-term applications, but it holds significant potential for advancing
quantum thermodynamic experiments on quantum computers, and ensuring the
faithful preparation of Gibbs states for various computational tasks.

The Python code for running the statevector simulations, using Qulacs [67],
can be found on GitHub [17].

## Acknowledgments

MC would like to thank Jacopo Settino, Andrea Giordano, Carlo Mastroianni,
Francesco Plastina, Salvatore Lorenzo, Sabrina Maniscalco, John Goold, and
Tony J. G. Apollaro for fruitful and extensive discussions on Gibbs states,
and the challenging task of variationally preparing them on NISQ computers. MC
acknowledges funding by TESS (Tertiary Education Scholarships Scheme), and
project QVAQT (Quantum Variational Algorithms for Quantum Technologies)
REP-2022-003 financed by the Malta Council for Science & Technology, for and
on behalf of the Foundation for Science and Technology, through the FUSION:
R&I Research Excellence Programme.

## References

  * [1] Ackley, D.H., Hinton, G.E., Sejnowski, T.J.: A learning algorithm for boltzmann machines. Cognitive Science 9(1), 147–169 (1985). https://doi.org/10.1016/s0364-0213(85)80012-4
  * [2] Aharonov, D., Arad, I., Vidick, T.: Guest Column: The Quantum PCP Conjecture. SIGACT News 44(2), 47–79 (Jun 2013). https://doi.org/10.1145/2491533.2491549
  * [3] Amin, M.H., Andriyash, E., Rolfe, J., Kulchytskyy, B., Melko, R.: Quantum Boltzmann Machine. Phys. Rev. X 8, 021050 (May 2018). https://doi.org/10.1103/PhysRevX.8.021050
  * [4] Anshu, A., Arunachalam, S., Kuwahara, T., Soleimanifar, M.: Sample-efficient learning of interacting quantum systems. Nature Physics 17(8), 931–935 (Aug 2021). https://doi.org/10.1038/s41567-021-01232-0
  * [5] Bardet, I., Capel, A., Gao, L., Lucia, A., Pérez-García, D., Rouzé, C.: Rapid Thermalization of Spin Chain Commuting Hamiltonians. Phys. Rev. Lett. 130, 060401 (Feb 2023). https://doi.org/10.1103/PhysRevLett.130.060401
  * [6] Biamonte, J., Wittek, P., Pancotti, N., Rebentrost, P., Wiebe, N., Lloyd, S.: Quantum machine learning. Nature 549(7671), 195–202 (Sep 2017). https://doi.org/10.1038/nature23474
  * [7] Bilgin, E., Boixo, S.: Preparing thermal states of quantum systems by dimension reduction. Phys. Rev. Lett. 105(17), 1–4 (2010). https://doi.org/10.1103/PhysRevLett.105.170405
  * [8] Brandão, F.G.S.L., Kastoryano, M.J.: Finite Correlation Length Implies Efficient Preparation of Quantum Thermal States. Communications in Mathematical Physics 365(1), 1–16 (Jan 2019). https://doi.org/10.1007/s00220-018-3150-8
  * [9] Brandão, F.G.S.L., Svore, K.: Quantum Speed-ups for Semidefinite Programming (2016). https://doi.org/10.48550/arxiv.1609.05537
  * [10] Chen, C.F., Kastoryano, M.J., ao, F.G.S.L.B., Gilyén, A.: Quantum Thermal State Preparation (2023). https://doi.org/10.48550/arXiv.2303.18224
  * [11] Chiang, C.F., Wocjan, P.: Quantum Algorithm for Preparing Thermal Gibbs States \- Detailed Analysis (2010). https://doi.org/10.48550/arXiv.1001.1130
  * [12] Childs, A.M., Maslov, D., Nam, Y., Ross, N.J., Su, Y.: Toward the first quantum simulation with quantum speedup. Proceedings of the National Academy of Sciences 115(38), 9456–9461 (2018). https://doi.org/10.1073/pnas.1801723115
  * [13] Chowdhury, A.N., Low, G.H., Wiebe, N.: A Variational Quantum Algorithm for Preparing Quantum Gibbs States (2020). https://doi.org/10.48550/arXiv.2002.00055
  * [14] Chowdhury, A.N., Somma, R.D.: Quantum algorithms for gibbs sampling and hitting-time estimation. Quantum Inf. Comput. 17(1-2), 41–64 (2017). https://doi.org/10.26421/qic17.1-2-3
  * [15] Cohn, J., Yang, F., Najafi, K., Jones, B., Freericks, J.K.: Minimal effective Gibbs ansatz: A simple protocol for extracting an accurate thermal representation for quantum simulation. Phys. Rev. A 102, 022622 (Aug 2020). https://doi.org/10.1103/PhysRevA.102.022622
  * [16] Consiglio, M., Settino, J., Giordano, A., Mastroianni, C., Plastina, F., Lorenzo, S., Maniscalco, S., Goold, J., Apollaro, T.J.G.: Variational Gibbs State Preparation on NISQ devices (2023). https://doi.org/10.48550/arXiv.2303.11276
  * [17] Consiglio, Mirko: Variational Gibbs State Preparation. https://github.com/mirkoconsiglio/VariationalGibbsStatePreparation (2023)
  * [18] Coopmans, L., Kikuchi, Y., Benedetti, M.: Predicting Gibbs-State Expectation Values with Pure Thermal Shadows. PRX Quantum 4, 010305 (Jan 2023). https://doi.org/10.1103/PRXQuantum.4.010305
  * [19] Davies, E.B.: Markovian master equations. Communications in Mathematical Physics 39(2), 91–110 (Jun 1974). https://doi.org/10.1007/bf01608389
  * [20] Davies, E.B.: Markovian master equations. II. Mathematische Annalen 219(2), 147–158 (Jun 1976). https://doi.org/10.1007/bf01351898
  * [21] Foldager, J., Pesah, A., Hansen, L.K.: Noise-assisted variational quantum thermalization. Sci. Rep. 12(1), 3862 (Mar 2022). https://doi.org/10.1038/s41598-022-07296-z
  * [22] França, D.S.: Perfect sampling for quantum Gibbs states. Quantum Inf. Comput. 18(5-6), 361–388 (2018). https://doi.org/10.26421/qic18.5-6-1
  * [23] Gacon, J., Zoufal, C., Carleo, G., Woerner, S.: Simultaneous Perturbation Stochastic Approximation of the Quantum Fisher Information. Quantum 5, 567 (Oct 2021). https://doi.org/10.22331/q-2021-10-20-567
  * [24] Gao, P., Jafferis, D.L., Wall, A.C.: Traversable wormholes via a double trace deformation. Journal of High Energy Physics 2017(12), 151 (Dec 2017). https://doi.org/10.1007/JHEP12(2017)151, https://doi.org/10.1007/JHEP12(2017)151
  * [25] Gao, P., Liu, H.: Regenesis and quantum traversable wormholes. Journal of High Energy Physics 2019(10), 48 (Oct 2019). https://doi.org/10.1007/JHEP10(2019)048, https://doi.org/10.1007/JHEP10(2019)048
  * [26] Getelina, J.C., Gomes, N., Iadecola, T., Orth, P.P., Yao, Y.X.: Adaptive variational quantum minimally entangled typical thermal states for finite temperature simulations (2023). https://doi.org/10.48550/arXiv.2301.02592
  * [27] Guo, X.Y., Li, S.S., Xiao, X., Xiang, Z.C., Ge, Z.Y., Li, H.K., Song, P.T., Peng, Y., Wang, Z., Xu, K., Zhang, P., Wang, L., Zheng, D.N., Fan, H.: Variational quantum simulation of thermal statistical states on a superconducting quantum processer. Chinese Physics B 32(1), 010307 (Jan 2023). https://doi.org/10.1088/1674-1056/aca7f3
  * [28] Hastings, W.K.: Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57(1), 97–109 (1970). https://doi.org/10.1093/biomet/57.1.97
  * [29] Haug, T., Bharti, K.: Generalized quantum assisted simulator. Quantum Science and Technology 7(4), 045019 (Aug 2022). https://doi.org/10.1088/2058-9565/ac83e7
  * [30] Holmes, Z., Muraleedharan, G., Somma, R.D., Subasi, Y., Şahinoğlu, B.: Quantum algorithms from fluctuation theorems: Thermal-state preparation. Quantum 6, 825 (Oct 2022). https://doi.org/10.22331/q-2022-10-06-825
  * [31] Huijgen, O., Coopmans, L., Najafi, P., Benedetti, M., Kappen, H.J.: Training Quantum Boltzmann Machines with the $\beta$-Variational Quantum Eigensolver (2023). https://doi.org/10.48550/arXiv.2304.08631
  * [32] Israel, W.: Thermo-field dynamics of black holes. Physics Letters A 57(2), 107–110 (1976). https://doi.org/https://doi.org/10.1016/0375-9601(76)90178-X, https://www.sciencedirect.com/science/article/pii/037596017690178X
  * [33] Jaklič, J., Prelovšek, P.: Lanczos method for the calculation of finite-temperature quantities in correlated systems. Phys. Rev. B 49, 5065–5068 (Feb 1994). https://doi.org/10.1103/PhysRevB.49.5065
  * [34] Kālis, M., Locāns, A., Šikovs, R., Naseri, H., Ambainis, A.: A hybrid quantum-classical approach for inference on restricted Boltzmann machines (2023). https://doi.org/10.48550/arXiv.2304.12418
  * [35] Kastoryano, M.J., ao, F.G.S.L.B.: Quantum Gibbs Samplers: the commuting case (2016). https://doi.org/10.48550/arXiv.1409.3435
  * [36] Kieferová, M., Wiebe, N.: Tomography and generative training with quantum Boltzmann machines. Phys. Rev. A 96, 062327 (Dec 2017). https://doi.org/10.1103/PhysRevA.96.062327
  * [37] Lee, C.K., Zhang, S.X., Hsieh, C.Y., Zhang, S., Shi, L.: Variational Quantum Simulations of Finite-Temperature Dynamical Properties via Thermofield Dynamics (2022). https://doi.org/10.48550/arXiv.2206.05571
  * [38] Lieb, E., Schultz, T., Mattis, D.: Two soluble models of an antiferromagnetic chain. Annals of Physics 16(3), 407–466 (1961). https://doi.org/10.1016/0003-4916(61)90115-4
  * [39] Liu, J.G., Mao, L., Zhang, P., Wang, L.: Solving quantum statistical mechanics with variational autoregressive networks and quantum circuits. Machine Learning: Science and Technology 2(2), 025011 (Feb 2021). https://doi.org/10.1088/2632-2153/aba19d
  * [40] Lu, S., Bañuls, M.C., Cirac, J.I.: Algorithms for Quantum Simulation at Finite Energies. PRX Quantum 2, 020321 (May 2021). https://doi.org/10.1103/PRXQuantum.2.020321
  * [41] Maldacena, J.: Eternal black holes in anti-de sitter. Journal of High Energy Physics 2003(04), 021 (apr 2003). https://doi.org/10.1088/1126-6708/2003/04/021, https://dx.doi.org/10.1088/1126-6708/2003/04/021
  * [42] Maldacena, J., Stanford, D., Yang, Z.: Diving into traversable wormholes. Fortschritte der Physik 65(5), 1700034 (2017). https://doi.org/https://doi.org/10.1002/prop.201700034, https://onlinelibrary.wiley.com/doi/abs/10.1002/prop.201700034
  * [43] Martyn, J., Swingle, B.: Product spectrum ansatz and the simplicity of thermal states. Phys. Rev. A 100(3), 1–10 (2019). https://doi.org/10.1103/PhysRevA.100.032107
  * [44] Matsui, T.: Purification and uniqueness of quantum Gibbs states. Communications in Mathematical Physics 162(2), 321–332 (1994). https://doi.org/10.1007/bf02102020
  * [45] McArdle, S., Jones, T., Endo, S., Li, Y., Benjamin, S.C., Yuan, X.: Variational ansatz-based quantum simulation of imaginary time evolution. npj Quantum Information 5(1), 75 (Sep 2019). https://doi.org/10.1038/s41534-019-0187-2
  * [46] Metcalf, M., Moussa, J.E., de Jong, W.A., Sarovar, M.: Engineered thermalization and cooling of quantum many-body systems. Phys. Rev. Res. 2, 023214 (May 2020). https://doi.org/10.1103/PhysRevResearch.2.023214
  * [47] Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H., Teller, E.: Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics 21(6), 1087–1092 (1953). https://doi.org/10.1063/1.1699114
  * [48] Motta, M., Sun, C., Tan, A.T., O’Rourke, M.J., Ye, E., Minnich, A.J., Brandão, F.G., Chan, G.K.L.: Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution. Nat. Phys. 16(2), 205–210 (2020). https://doi.org/10.1038/s41567-019-0704-4
  * [49] Nielsen, M.A., Chuang, I.L.: Quantum Computation and Quantum Information: 10th Anniversary Edition. Cambridge University Press (2010). https://doi.org/10.1017/cbo9780511976667
  * [50] Nocedal, J., Wright, S.J.: Numerical Optimization. Springer, New York, NY, USA, 2e edn. (2006). https://doi.org/10.1007/978-0-387-40065-5
  * [51] Ozols, M., Roetteler, M., Roland, J.: Quantum Rejection Sampling. In: Proceedings of the 3rd Innovations in Theoretical Computer Science Conference. p. 290–308. Itcs ’12, Association for Computing Machinery, New York, NY, USA (2012). https://doi.org/10.1145/2090236.2090261
  * [52] Poulin, D., Wocjan, P.: Sampling from the thermal quantum Gibbs state and evaluating partition functions with a quantum computer. Phys. Rev. Lett. 103(22), 1–4 (2009). https://doi.org/10.1103/PhysRevLett.103.220502
  * [53] Premaratne, S.P., Matsuura, A.Y.: Engineering a Cost Function for Real-world Implementation of a Variational Quantum Algorithm. In: 2020 IEEE International Conference on Quantum Computing and Engineering (QCE). pp. 278–285 (2020). https://doi.org/10.1109/qce49297.2020.00042
  * [54] Propp, J.G., Wilson, D.B.: Exact sampling with coupled Markov chains and applications to statistical mechanics. Random Structures & Algorithms 9(1-2), 223–252 (1996). https://doi.org/10.1002/(sici)1098-2418(199608/09)9:1/2¡223::aid-rsa14¿3.0.co;2-o
  * [55] Rall, P., Wang, C., Wocjan, P.: Thermal State Preparation via Rounding Promises (2022). https://doi.org/10.48550/arXiv.2210.01670
  * [56] Riera, A., Gogolin, C., Eisert, J.: Thermalization in nature and on a quantum computer. Phys. Rev. Lett. 108(8) (2012). https://doi.org/10.1103/PhysRevLett.108.080402
  * [57] Sagastizabal, R., Premaratne, S.P., Klaver, B.A., Rol, M.A., Negîrneac, V., Moreira, M.S., Zou, X., Johri, S., Muthusubramanian, N., Beekman, M., Zachariadis, C., Ostroukh, V.P., Haider, N., Bruno, A., Matsuura, A.Y., DiCarlo, L.: Variational preparation of finite-temperature states on a quantum computer. npj Quantum Inf. 7(1), 1–7 (2021). https://doi.org/10.1038/s41534-021-00468-1
  * [58] Sewell, T.J., White, C.D., Swingle, B.: Thermal Multi-scale Entanglement Renormalization Ansatz for Variational Gibbs State Preparation (2022). https://doi.org/10.48550/arXiv.2210.16419
  * [59] Shabani, A., Neven, H.: Artificial quantum thermal bath: Engineering temperature for a many-body quantum system. Phys. Rev. A 94, 052301 (Nov 2016). https://doi.org/10.1103/PhysRevA.94.052301
  * [60] Shor, P.: Algorithms for quantum computation: discrete logarithms and factoring. In: Proceedings 35th Annual Symposium on Foundations of Computer Science. pp. 124–134 (1994). https://doi.org/10.1109/sfcs.1994.365700
  * [61] Shtanko, O., Movassagh, R.: Algorithms for Gibbs state preparation on noiseless and noisy random quantum circuits (2021). https://doi.org/10.48550/arXiv.2112.14688
  * [62] Silva, T.L., Taddei, M.M., Carrazza, S., Aolita, L.: Fragmented imaginary-time evolution for early-stage quantum signal processors (2022). https://doi.org/10.48550/arXiv.2110.13180
  * [63] Somma, R.D., Boixo, S., Barnum, H., Knill, E.: Quantum simulations of classical annealing processes. Phys. Rev. Lett. 101(13), 1–4 (2008). https://doi.org/10.1103/PhysRevLett.101.130504
  * [64] Sugiura, S., Shimizu, A.: Thermal Pure Quantum States at Finite Temperature. Phys. Rev. Lett. 108, 240401 (Jun 2012). https://doi.org/10.1103/PhysRevLett.108.240401
  * [65] Sugiura, S., Shimizu, A.: Canonical Thermal Pure Quantum State. Phys. Rev. Lett. 111, 010401 (Jul 2013). https://doi.org/10.1103/PhysRevLett.111.010401
  * [66] Sun, S.N., Motta, M., Tazhigulov, R.N., Tan, A.T., Chan, G.K.L., Minnich, A.J.: Quantum Computation of Finite-Temperature Static and Dynamical Properties of Spin Systems Using Quantum Imaginary Time Evolution. PRX Quantum 2, 010317 (Feb 2021). https://doi.org/10.1103/PRXQuantum.2.010317
  * [67] Suzuki, Y., Kawase, Y., Masumura, Y., Hiraga, Y., Nakadai, M., Chen, J., Nakanishi, K.M., Mitarai, K., Imai, R., Tamiya, S., Yamamoto, T., Yan, T., Kawakubo, T., Nakagawa, Y.O., Ibe, Y., Zhang, Y., Yamashita, H., Yoshimura, H., Hayashi, A., Fujii, K.: Qulacs: A fast and versatile quantum circuit simulator for research purpose. Quantum 5, 559 (Oct 2021). https://doi.org/10.22331/q-2021-10-06-559
  * [68] Tan, K.C.: Fast quantum imaginary time evolution (2020). https://doi.org/10.48550/arXiv.2009.1223
  * [69] Temme, K., Osborne, T.J., Vollbrecht, K.G., Poulin, D., Verstraete, F.: Quantum Metropolis sampling. Nature 471(7336), 87–90 (2011). https://doi.org/10.1038/nature09770
  * [70] Terhal, B.M., DiVincenzo, D.P.: Problem of equilibration and the computation of correlation functions on a quantum computer. Phys. Rev. A 61, 022301 (Jan 2000). https://doi.org/10.1103/PhysRevA.61.022301
  * [71] Trabesinger, A.: Quantum simulation. Nature Physics 8(4), 263–263 (Apr 2012). https://doi.org/10.1038/nphys2258
  * [72] Uhlmann, A.: Transition Probability (Fidelity) and Its Relatives. Foundations of Physics 41(3), 288–298 (2011). https://doi.org/10.1007/s10701-009-9381-y
  * [73] Vedral, V.: The role of relative entropy in quantum information theory. Rev. Mod. Phys. 74, 197–234 (Mar 2002). https://doi.org/10.1103/RevModPhys.74.197
  * [74] Verdon, G., Marks, J., Nanda, S., Leichenauer, S., Hidary, J.: Quantum Hamiltonian-Based Models and the Variational Quantum Thermalizer Algorithm (2019). https://doi.org/10.48550/arXiv.1910.02071
  * [75] Wang, X., Feng, X., Hartung, T., Jansen, K., Stornati, P.: Critical behavior of Ising model by preparing thermal state on quantum computer (2023). https://doi.org/10.48550/arXiv.2302.14279
  * [76] Wang, Y., Li, G., Wang, X.: Variational Quantum Gibbs State Preparation with a Truncated Taylor Series. Phys. Rev. Appl. 16(5), 1 (2021). https://doi.org/10.1103/PhysRevApplied.16.054035
  * [77] Warren, A., Zhu, L., Mayhall, N.J., Barnes, E., Economou, S.E.: Adaptive variational algorithms for quantum Gibbs state preparation (2022). https://doi.org/10.48550/arXiv.2203.12757
  * [78] Watrous, J.: Quantum Computational Complexity (2008). https://doi.org/10.48550/arxiv.0804.3401
  * [79] Wocjan, P., Temme, K.: Szegedy Walk Unitaries for Quantum Maps (2021). https://doi.org/10.48550/arXiv.2107.07365
  * [80] Wu, J., Hsieh, T.H.: Variational Thermal Quantum Simulation via Thermofield Double States. Phys. Rev. Lett. 123(22), 1–7 (2019). https://doi.org/10.1103/PhysRevLett.123.220502
  * [81] Yuan, X., Endo, S., Zhao, Q., Li, Y., Benjamin, S.C.: Theory of variational quantum simulation. Quantum 3, 1–41 (2019). https://doi.org/10.22331/q-2019-10-07-191
  * [82] Yung, M.H., Aspuru-Guzik, A.: A quantum-quantum Metropolis algorithm. Proc. Natl. Acad. Sci. U. S. A. 109(3), 754–759 (2012). https://doi.org/10.1073/pnas.1111758109
  * [83] Zhang, D., Bosse, J.L., Cubitt, T.: Dissipative Quantum Gibbs Sampling (2023). https://doi.org/10.48550/arXiv.2304.04526
  * [84] Zhu, D., Johri, S., Linke, N.M., Landsman, K.A., Huerta Alderete, C., Nguyen, N.H., Matsuura, A.Y., Hsieh, T.H., Monroe, C.: Generation of thermofield double states and critical ground states with a quantum computer. Proc. Natl. Acad. Sci. U. S. A. 117(41), 25402–25406 (2020). https://doi.org/10.1073/pnas.2006337117
  * [85] Zoufal, C., Lucchi, A., Woerner, S.: Variational quantum Boltzmann machines. Quantum Mach. Intell. 3(1), 1–15 (2021). https://doi.org/10.1007/s42484-020-00033-7

  *[NISQ]: noisy intermediate-scale quantum
  *[VQA]: variational quantum algorithm
  *[VQAs]: variational quantum algorithm
  *[PQC]: parametrized quantum circuit
  *[QAOA]: quantum approximate optimization algorithm
  *[ADAPT]: adaptive derivative assembled problem-tailored
  *[TFD]: thermofield double
  *[QITE]: quantum imaginary time evolution
  *[METTS]: minimally entangled typical thermal states
  *[CBMs]: classical Boltzmann machine
  *[CBM]: classical Boltzmann machine
  *[QBM]: quantum Boltzmann machine
  *[QBMs]: quantum Boltzmann machine
  *[BFGS]: Broyden–Fletcher–Goldfarb–Shanno