  {Many-body quantum states at thermal equilibrium are ubiquitous in nature.
Investigating their dynamical properties is a formidable task due to the
complexity of the Hilbert space they live in. Quantum computers may have the
potential to effectively simulate quantum systems, provided that the many-body
state under scrutiny can be faithfully prepared via an efficient algorithm.
With this aim, we present a hybrid quantum--classical variational quantum
algorithm for the preparation of the Gibbs state of the quantum $XY$ model. Our
algorithm is based on the Grover and Rudolph parametrized quantum circuit for
the preparation of the Boltzmann weights of the Gibbs state, and on a
parity-preserving ansatz for the allocation of the eigenenergy basis to their
respective Boltzmann weight. We explicitly show, with a paradigmatic few-body
case instance, how the symmetries of a many-body system can be exploited to
significantly reduce the exponentially increasing number of variational
parameters needed in the Grover and Rudolph algorithm. Finally, we show that
the density matrix, of the Gibbs state of the $XY$ model, obtained by
statevector simulations for different parameters, exhibits a fidelity close to
unity with the exact Gibbs state; this highlights the potential use of our
protocol on current quantum computers.
