# Variational simulation of $d$-level systems on qubit-based quantum
simulators

Chufan Lyu  1  Institute of Fundamental and Frontier Sciences, University of
Electronic Science and Technology of China, Chengdu 611731, China  1  Key
Laboratory of Quantum Physics and Photonic Quantum Information, Ministry of
Education, University of Electronic Science and Technology of China, Chengdu
611731, China  Zuoheng Zou  1  Central Research Institute, Huawei
Technologies, Shenzhen, 518129, China  Xusheng Xu  1  Central Research
Institute, Huawei Technologies, Shenzhen, 518129, China  Man-Hong Yung  1
Central Research Institute, Huawei Technologies, Shenzhen, 518129, China  1
Shenzhen Institute for Quantum Science and Engineering, Southern University of
Science and Technology, Shenzhen, 518055, China  1  International Quantum
Academy, Shenzhen, 518048, China  1  Guangdong Provincial Key Laboratory of
Quantum Science and Engineering, Southern University of Science and
Technology, Shenzhen, 518055, China  1  Shenzhen Key Laboratory of Quantum
Science and Engineering, Southern University of Science and Technology,
Shenzhen, 518055, China  Abolfazl Bayat<EMAIL_ADDRESS>1
Institute of Fundamental and Frontier Sciences, University of Electronic
Science and Technology of China, Chengdu 611731, China  1  Key Laboratory of
Quantum Physics and Photonic Quantum Information, Ministry of Education,
University of Electronic Science and Technology of China, Chengdu 611731,
China

###### Abstract

Current quantum simulators are primarily qubit-based, making them naturally
suitable for simulating 2-level quantum systems. However, many systems in
nature are inherently $d$-level, including higher spins, bosons, vibrational
modes, and itinerant electrons. To simulate $d$-level systems on qubit-based
quantum simulators, an encoding method is required to map the $d$-level system
onto a qubit basis. Such mapping may introduce illegitimate states in the
Hilbert space which makes the simulation more sophisticated. In this paper, we
develop a systematic method to address the illegitimate states. In addition,
we compare two different mappings, namely binary and symmetry encoding
methods, and compare their performance through variational simulation of the
ground state and time evolution of various many-body systems. While binary
encoding is very efficient with respect to the number of qubits it cannot
easily incorporate the symmetries of the original Hamiltonian in its circuit
design. On the other hand, the symmetry encoding facilitates the
implementation of symmetries in the circuit design, though it comes with an
overhead for the number of qubits. Our analysis shows that the symmetry
encoding significantly outperforms the binary encoding, despite requiring
extra qubits. Their advantage is indicated by requiring fewer two-qubit gates,
converging faster, and being far more resilient to Barren plateaus. We have
performed variational ground state simulations of spin-1, spin-3/2, and
bosonic systems as well as variational time evolution of spin-1 systems. Our
proposal can be implemented on existing quantum simulators and its potential
is extendable to a broad class of physical models.

##  I Introduction

The importance of quantum computation lies in their capability for exponential
speedup in simulating quantum systems as well as addressing certain
computational problems. Quantum computers are programmable machines that can
implement any unitary operation on quantum systems and thus convert any
quantum state in the Hilbert space to another. Any unitary operation can be
decomposed into a combination of single-particle rotations as well as a two-
particle entangling gate, the so-called universal gates. This creates the
concept of digital quantum computation in which the design of a quantum
computer reduces to implementation of the universal gates. Near-term noisy
intermediate scale quantum (NISQ) devices, however, suffer from various
imperfection, including limited connectivity, short coherence time, noisy gate
operations, absence of quantum error correction, imperfect initialization, and
faulty readout [1]. These imperfections make universal quantum computation
inaccessible, at least, in a near future for NISQ devices. Nonetheless,
quantum simulators which are programmable machines with the capability of
realizing a limited class of unitary operations for emulating complex quantum
systems are feasible [2, 3]. Thanks to recent advancement in quantum
technologies quantum simulators are emerging in various physical platforms,
including optical lattices [4, 5, 6], Rydberg atoms [7, 8], optical chips [9,
10, 11], solid state systems [12], ion-traps [13, 14, 15, 16], and
superconducting devises [17, 18, 19, 20, 21, 22, 23, 24, 25]. These
implementations are mostly qubit-based, making them naturally suitable for
simulating two-level quantum systems.

Many systems in nature are inherently $d$-level (with $d>2$), including higher
spins [26], bosons [27], vibrational modes [28], and itinerant electrons [29].
In addition, there are classical optimization problems which can efficiently
mapped to qudit systems [30, 31]. In principle, there are two methods for
simulating such systems on quantum simulators. First, one can design and build
$d$-level quantum simulators which can directly emulate those systems [32, 33,
34, 35, 36, 37, 38, 39, 40]. Second, one can rely on existing qubit-based
quantum simulators but exploit proper encoding to effectively reproduce the
behavior of $d$-level systems. The former approach is still technologically
challenging and qudit quantum gates have only been considered for $d{=}3$ and
$d{=}4$ with several practical difficulties [32, 33, 34, 35, 36, 37, 38]. The
latter approach requires suitable encoding methods to map $d$-level systems to
qubits. In Refs. [41, 42, 43, 44, 45], several encoding methods have been
introduced to map a $d$-level system onto a qubit system. In these methods,
depending on the mapping algorithm, a single $d$-level system can be mapped to
a different number of qubits. While the most efficient mapping requires
$\lceil\log_{2}(d)\rceil$ qubits, the other may demand $d$ qubits. Note that,
the complexity does not solely depend on the mapping as the manipulation of
the qubits might be easier for a less efficient encoding. Simulating $d$-level
systems on qubit-based quantum simulators may face several challenges: (i)
mapping $d$-level systems to qubit systems can enlarge the corresponding
Hilbert space, implying that certain regions of the Hilbert space does not
represent the original system; (ii) the resulting qubit Hamiltonian may
include multi-body interactions, which complicates the quantum simulation;
(iii) symmetries of the original Hamiltonian may take complex forms after
mapping to the qubits.

While universal quantum computation is not achievable with NISQ simulators,
one may wonder whether such quantum devices can provide any computational
advantage over classical computers. Indeed, the advantage of NISQ simulators
has been demonstrated in sampling problems [10, 46, 47]. However, sampling
problems do not have practical applications. Hence, achieving practical
advantage with NISQ simulators is still an open problem. Variational Quantum
Algorithms (VQAs) are the most promising approach for achieving quantum
advantage in NISQ era. VQAs have been developed for a range of applications,
including quantum chemistry [48, 49, 50, 15, 51, 52], combinatorial
optimization [53, 54, 55, 30], quantum machine learning [56, 57, 58, 59, 60,
61, 62, 63], dynamical simulations in closed [64, 65, 66, 67, 68] and open
[69, 70, 71, 66, 72] systems, quantum sensing [73, 74, 75, 76, 77, 78, 79],
and condensed matter physics [16, 80, 81]. Variational Quantum Eigensolver
(VQE) is an important class of VQAs which tend to simulate the low-energy
eigenstates of a given Hamiltonian [50, 82, 83, 84, 85, 86, 87]. Exploiting
the inherent symmetries of the Hamiltonian in the VQE protocol, either through
tailored circuit design [88, 89, 81, 90, 91, 92, 93, 94] or by incorporating
them into the cost function as penalty terms [95, 96], is extremely beneficial
in simplifying the protocol. However, most of these symmetry exploiting
methods have been focused on qubit systems. Mapping symmetries of $d$-level
quantum systems into qubit quantum circuits remains unexplored.

In this paper, we focus on the variational simulation of $d$-level systems on
qubit-based quantum simulators. We first show that mapping from $d$-level
systems to qubits may result in introduction of illegitimate states in the
qubit Hilbert space which have no correspondence in the original system. We
provide a systematic penalty approach to prevent these states to spoil our
variational simulations. We then consider two different mappings, namely
binary and symmetry encoding methods. Binary encoding is efficient with
respect to the number of qubits, though, it cannot incorporate symmetries in
the circuit design. The symmetry encoding provides a natural way for
incorporating the symmetries of the original Hamiltonian into the circuit of
the qubit-based quantum simulators, although with an overhead of requiring
extra qubits. Through comprehensive variational quantum simulation of the
ground state and non-equilibrium dynamics of several spin and bosonic many-
body systems, we demonstrate that symmetry encoding outperforms the binary
mapping. This advantage reveals itself in fewer two-qubit gates, faster
convergence, and less sensitivity to training issues such as Barren plateaus.
Our proposal offers a promising and timely approach for simulating $d$-level
systems on qubit-based platforms, with potential applications extendable to
other types of quantum simulators.

##  II Variational Quantum Algorithms

Variational Quantum Algorithms (VQAs) have emerged as a promising class of
algorithms that leverage the hybrid quantum-classical approach to solve
optimization problems and simulate quantum systems, particularly beneficial in
the Noisy Intermediate-Scale Quantum (NISQ) era. In this section, we provide a
brief introduction to the VQAs, focusing on their applications in preparing
the ground state of a given Hamiltonian, or simulating non-equilibrium
dynamics using the shallowest possible quantum circuit.

The Variational Quantum Eigensolver (VQE) is a pivotal example of VQA
dedicated to finding the low energy eigenstates (e.g. the ground state) of a
Hamiltonian, which is central to understanding the physical properties of
quantum systems. In the VQE algorithm, a quantum state, denoted as
$|\psi(\vec{\theta})\rangle$, is prepared through a parameterized quantum
circuit acting on an initial state, expressed as
$|\psi(\vec{\theta})\rangle{=}U(\vec{\theta})|\psi_{0}\rangle$, for a given N
qubits. This parameterized quantum circuit, often referred as the ansatz,
incorporates a set of tunable parameters
$\vec{\theta}=(\theta_{1},\theta_{2},\dots,\theta_{L})$, and
$|\psi_{0}\rangle$ is the initial state. Varying $\vec{\theta}$ enables the
output of the quantum circuit, namely $|\psi(\vec{\theta})\rangle$, to explore
a segment of the Hilbert space, which is called reachable set. Notably, in
circuits of significant depth, and thus a large number of parameters $L$, the
reachable set can span the entire Hilbert space, allowing the preparation of
any quantum state of $N$ qubits. However, the goal is to keep the circuit as
shallow as possible, focusing on the most relevant part of the Hilbert space.
Specifically, VQE aims for a comparatively shallow ansatz that the reachable
set contains the ground state of the Hamiltonian of interest $H$. Various
ansatzes have been proposed, each differing in complexity [50, 97, 98, 99,
100, 101]. The expectation value of the Hamiltonian $\langle
H\rangle=\langle\psi(\vec{\theta})|H|\psi(\vec{\theta})\rangle$ is measured at
the output of the simulator. One then iteratively minimize $\langle H\rangle$
with respect to the parameters $\vec{\theta}$ using a classical optimization
algorithm, such as the gradient decent [102]. After a few iterations, $\langle
H\rangle$ reaches its minimum for a special set of parameters
$\vec{\theta}^{*}$. If the reachable set contains the ground state then the
variational circuit with optimal parameters $\vec{\theta}^{*}$ produces the
ground state of the Hamiltonian $H$.

The simulation of the ground state usually requires fewer resources than the
non-equilibrium dynamics of many-body systems. This is because the
entanglement in the ground state of local Hamiltonian is area-law while non-
equilibrium states usually lead to the emergence of far more complex volume-
law entangled states. Non-equilibrium dynamics of a quantum system is
described by Schrödinger equation:

$i\frac{d}{dt}|\psi(t)\rangle=H|\psi(t)\rangle.$ (1)

The solution for a time-independent Hamiltonian is given by
$|\psi(t)\rangle{=}e^{-iHt}|\psi_{0}\rangle$. There are two different methods
to simulate non-equilibrium dynamics of a many-body system on a digital
quantum simulator: (i) Trotterization of the unitary operation $e^{-iHt}$; and
(ii) variational simulation of the non-equilibrium dynamics. The former
requires a deep quantum circuit to keep the Trotterization error below an
acceptable range [81]. The latter is possibly more NISQ friendly as the use of
a classical optimization may help to simplify the quantum circuit. In this
paper, we focus on a variational algorithm, developed in Refs. [103, 66], for
simulating non-equilibrium dynamics of $d$-level systems on qubit-based
quantum computers. We employ a parameterized trial state
$|\psi(\vec{\theta}(t))\rangle$ with time-dependent parameters
$\vec{\theta}(t)$. For variational quantum simulation, the trial state is
prepared by a variational quantum circuit $U(\vec{\theta}(t))$ acting on an
initial state $|\psi_{0}\rangle$, e.g.
$|\psi(\vec{\theta}(t))\rangle=U(\vec{\theta}(t))|\psi_{0}\rangle$. Assuming
the state at time $t$ is obtained by
$|\psi(t)\rangle=|\psi(\vec{\theta}(t))\rangle$, the goal is to approximate
the state $|\psi(t+\delta t)\rangle$. According to the Schrödinger equation,
the state evolves from time $t$ to $t+\delta t$ as follows:

$|\psi(t+\delta t)\rangle\approx|\psi(\vec{\theta}(t))\rangle-i\delta
tH|\psi(\vec{\theta}(t))\rangle.$ (2)

Given that the variational trial state $|\psi(\vec{\theta}(t+\delta
t))\rangle$ is trying to find the optimal solution $\vec{\theta}(t+\delta t)$
that best approximates the state within its reachable set, we have:

$|\psi(\vec{\theta}(t+\delta
t))\rangle\approx|\psi(\vec{\theta}(t))\rangle-i\delta
tH|\psi(\vec{\theta}(t))\rangle.$ (3)

Note that, by focusing on the variation of the trial state with changing
parameters, we have:

$|\psi(\vec{\theta}(t+\delta
t))\rangle\approx|\psi(\vec{\theta}(t))\rangle+\sum_{j}\frac{\partial|\psi(\vec{\theta}(t))\rangle}{\partial\theta_{j}}\delta\theta_{j},$
(4)

where $\delta\theta_{j}=\theta_{j}(t+\delta t)-\theta_{j}(t)$. The objective
is then to find the optimal $\vec{\theta}(t+\delta t)$ that minimizes the
difference between the evolved state according to the Schrödinger equation and
the variational approximation:

$\sum_{j}\frac{\partial|\psi(\vec{\theta}(t))\rangle}{\partial\theta_{j}}\delta\theta_{j}\approx-i\delta
tH|\psi(\vec{\theta}(t))\rangle.$ (5)

Assuming the parameters are real, applying McLachlan’s variational principle
gives:

$\sum_{j}A_{i,j}^{R}\dot{\theta_{j}}=C_{i}^{I},$ (6)

where $A$ and $\vec{C}$ are defined as:

$\displaystyle
A_{i,j}=\frac{\partial\langle\psi(\vec{\theta}(t))|}{\partial\theta_{i}}\frac{\partial|\psi(\vec{\theta}(t))\rangle}{\partial\theta_{j}},$
(7) $\displaystyle
C_{i}=\frac{\partial\langle\psi(\vec{\theta}(t))|}{\partial\theta_{i}}H|\psi(\vec{\theta}(t))\rangle.$

Here, $A_{ij}^{R}$ and $C_{i}^{I}$ represent the real and imaginary parts of
$A_{ij}$ and $C_{i}$, respectively. The derivatives of the time-dependent
parameters, $\vec{\dot{\theta}}$, are thus given by
$\vec{\dot{\theta}}={A^{R}}^{-1}\vec{C}^{I}$. By updating the parameters as:

$\vec{\theta}(t+\delta t)=\vec{\theta}(t)+\delta t{A^{R}}^{-1}\vec{C}^{I},$
(8)

and repeating this process $T/\delta t$ times, the time evolution over a
duration $T$ is simulated. Fig. 1(e) illustrates the circuit implementation
for evaluating the real and imaginary parts of the matrix $A$ and vector $C$.

In the domain of VQAs, resource efficiency emerges as a critical property due
to the constraints imposed by NISQ devices. These devices are typically
characterized by a limited number of qubits, brief coherence times, and
elevated error rates. As a result, a quantification of resource efficiency is
essential for assessing the performance of VQAs. This efficiency can be
dissected into two primary categories: classical and quantum resources.
Classical resources primarily depend on the computational demands of the
optimization process, notably the number of iterations needed to tune the
parameters. Furthermore, each parameter requires individual measurements for
gradient evaluation if the optimizer in the VQAs adopted to gradient-based
optimizer. Thus, we define the measure of Classical Resources ($C_{R}$) as the
product of the number of parameters $L$ and the number of optimization
iterations $n_{I}$, expressed mathematically as:

$C_{R}=L\times n_{I}$ (9)

For quantum resources, the emphasis is on the circuit components that are most
susceptible to noises on NISQ devices. Since single-qubit operations are
relatively error-resilient and precise, the focus shifts to two-qubit
entangling operations, which are more prone to errors. Hence, we quantify
quantum resources by the number of two-qubit entangling gates used in the
circuit.

(a)

(b)

(c)

(d)

(e)

Figure 1: (a) The hardware efficient ansatz circuit for system size $N{=}4$.
The circuits highlighted in blue can be repeated to increase the
expressibility of the variational circuit. (b) Circuit implementing of
$\mathcal{N}(\theta_{x},\theta_{y},\theta_{z})$ as an entangling gate between
two qubits. (c) The $s_{z}$-conserving ansatz circuit for system size $N{=}4$.
(d) The total spin conserving ansatz circuit for system size $N{=}4$. (e)
Circuit implementation for evaluating the coefficients in a quantum simulator.
The qubit on the top is an ancillary qubit. The probability of measuring the
ancillary qubit in the state $|0\rangle$ is:
$P_{0}=\left(1+\Re(\langle\psi|U_{i}^{\dagger}U_{j}|\psi\rangle)\right)/2$.
One can utilize the same circuit with a phase gate
$\text{S}=\begin{pmatrix}1&0\\\ 0&i\end{pmatrix}$ applied to the ancillary
qubit before the last Hadamard gate to measure the imaginary part of the
coefficient.

Therefore, in consideration of developing a VQA, selecting an appropriate
ansatz is essential for achieving an optimal balance between minimizing the
use of error-prone two-qubit gates and controlling the number of variational
parameters. The former is crucial for quantum resource efficiency, and the
latter not only dictates the expansiveness of the reachable set but also
influences the efficiency of classical resources.

One of the most commonly used ansatzes in the literature is the hardware-
efficient ansatz, depicted in Fig. 1(a). This ansatz consists of alternating
layers of single-qubit rotations and two-qubit entangling gates. The single-
qubit rotations are defined by the formula:

$R_{\alpha}(\theta)=e^{-i\theta\sigma^{\alpha}/2},$ (10)

where $\sigma^{\alpha}$ is the Pauli operator, with $\alpha\in\\{x,y,z\\}$.
For the two-qubit operations, we employ the Controlled-Not (CNOT) gate, which
is specified by the matrix:

$\text{CNOT}=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&0&1\\\
0&0&1&0\end{pmatrix}.$ (11)

This structure is particularly favored for its simplicity and efficiency in
creating entanglement between qubits, which is critical for the success of
various variational quantum algorithms. Nonetheless, the hardware-efficient
ansatz suffers from Barren plateau [104, 105], which makes its training
challenging for large-scale quantum systems. Furthermore, this ansatz does not
inherently preserve any symmetries that may exist in the Hamiltonian, which
could enhance quantum simulation performance. As an alternative, one might
consider a more complex entangling gate with tunable parameters, such as
[106]:

$\mathcal{N}(\theta_{x},\theta_{y},\theta_{z})=e^{i(\theta_{x}\sigma^{x}\otimes\sigma^{x}+\theta_{y}\sigma^{y}\otimes\sigma^{y}+\theta_{z}\sigma^{z}\otimes\sigma^{z})},$
(12)

where $\sigma^{x}\otimes\sigma^{x}$, $\sigma^{y}\otimes\sigma^{y}$, and
$\sigma^{z}\otimes\sigma^{z}$ are the tensor products of the Pauli operators
acting on two qubits. The circuit implementation of this unitary operation is
depicted in Fig. 1(b). In the special case where
$\theta_{x}=\theta_{y}=\theta_{z}=\theta$, this unitary conserves the number
of excitations as well as the total spin. By combining this entangling gate
with single-qubit rotations, for example $R_{z}(\theta)$, one can construct an
excitation preserving ansatz, $U(\vec{\theta})$, satisfying
$[U(\vec{\theta}),S^{z}_{tot}]=0$, where
$S^{\alpha}_{tot}=\frac{1}{2}\sum_{i}\sigma^{\alpha}_{i}$ for
$\alpha\in\\{x,y,z\\}$, as illustrated in Figs. 1(c). In the absence of
single-qubit rotations on the $z$-axis, namely $R_{z}(\theta)$, the ansatz
preserves the total spin, i.e. $[U(\vec{\theta}),S_{\text{tot}}^{2}]=0$, where
$S_{\text{tot}}^{2}={(S_{tot}^{x})}^{2}+{(S_{tot}^{y})}^{2}+{(S_{tot}^{z})}^{2}$,
shown in Fig. 1(d).

##  III Qudit to Qubit Mapping

The first step towards the simulation of a $d$-level many-body quantum system
on a qubit-based quantum simulator is to exploit an encoding algorithm which
maps the $d$-level operators into qubit operators. Several encoding techniques
have been explored [41, 42, 43, 44]. The most straightforward method is the
binary encoding in which a $d$-level particle is mapped to
$M{=}\lceil\log_{2}(d)\rceil$ qubits, where $\lceil x\rceil$ represents the
smallest integer number that is greater or equal to the real number $x$. In
this encoding, each $d$-level basis state $|\widetilde{j}\rangle$
($j=0,1,\cdots,d-1$) is mapped to the binary representation of its index $j$,
which is converted into an $M$-qubit computational state. The binary
representation is padded with zeros on the left to match the number of qubits
$n$:

$|\widetilde{j}\rangle\mapsto|\mathbf{b}_{M}^{(j)}\rangle=|b_{M-1}b_{M-2}\ldots
b_{1}b_{0}\rangle$ (13)

where $b_{M-1}b_{M-2}\ldots b_{1}b_{0}$ is the binary representation of the
integer $j$ padded to $M$ bits. For instance, for a $9$-level qudit system, we
require $4$ qubits to encode $|\widetilde{0}\rangle=|0000\rangle$ all the way
to $|\widetilde{8}\rangle=|1000\rangle$. Note that throughout this paper the
tilde basis $|\widetilde{j}\rangle$ is used for qudits and the ordinary basis
$|j\rangle$ is used for qubits. It will be convenient to define the
transformation operator between the two bases as

$\mathcal{T}_{b}=\sum_{j=0}^{d-1}|\mathbf{b}_{M}^{(j)}\rangle\langle\widetilde{j}|.$
(14)

Note that while we always have
$\mathcal{T}_{b}^{\dagger}\mathcal{T}_{b}=\mathbb{I}_{\mathrm{qudits}}$, the
vice versa is not always valid, i.e.
$\mathcal{T}_{b}\mathcal{T}_{b}^{\dagger}\neq\mathbb{I}_{\mathrm{qubits}}$,
unless $M{=}\log_{2}(d)$. By using the transformation matrix
$\mathcal{T}_{b}$, one can map any qudit operator $\widetilde{H}$ into a qubit
basis through $H=\mathcal{T}_{b}\widetilde{H}\mathcal{T}_{b}^{\dagger}$ and
then write it in terms of Pauli strings with the help of the following maps:

$\displaystyle|0\rangle\langle 0|$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\mathbb{I}+\sigma^{z}),\;\;\;\;\;\;|1\rangle\langle
1|=\frac{1}{2}(\mathbb{I}-\sigma^{z}),$ $\displaystyle|0\rangle\langle 1|$
$\displaystyle=$
$\displaystyle\frac{1}{2}(\sigma^{x}+i\sigma^{y}),\;\;|1\rangle\langle
0|=\frac{1}{2}(\sigma^{x}-i\sigma^{y}),$ (15)

where $\mathbb{I}$ is the identity and $\sigma^{x,y,z}$ are the Pauli
operators.

If the qudit Hamiltonian $\widetilde{H}$ supports a symmetry $\widetilde{S}$,
namely $[\widetilde{H},\widetilde{S}]{=}0$, the corresponding qubit
Hamiltonian $H=\mathcal{T}_{b}\widetilde{H}\mathcal{T}_{b}^{\dagger}$ supports
an equivalent symmetry
$S=\mathcal{T}_{b}\widetilde{S}\mathcal{T}_{b}^{\dagger}$. However, even if
$\widetilde{S}$ might be a simple local operator the form of $S$ might be
complex consisting of several multi-qubit terms in the qubit basis. It has
been shown that in VQE simulations, the most resource efficient approach to
benefit from the symmetries is to incorporate them in the design of the
circuit such that $[U(\vec{\theta}),S]{=}0$ [97]. In this situation, the
circuit preserves the symmetries in the initial state and thus makes the VQE
simulation easier. The problem is that if the symmetry operator $S$ is complex
itself then designing such a circuit becomes challenging. This is often the
case for mapping qudits to qubits. Therefore, an alternative mapping might be
exploited, namely symmetry encoding. This approach maps each $d$-level basis
state $|\widetilde{j}\rangle$ into a superposition of $d{-}1$ computational
qubit states with exactly $j$ excitations:

$|\widetilde{j}\rangle\mapsto|\mathbf{f}^{(j)}_{d-1}\rangle=\frac{1}{\sqrt{\binom{d-1}{j}}}\sum_{\begin{subarray}{c}\forall
k\text{ with }j\\\
\text{excitations}\end{subarray}}|\mathbf{b}_{d-1}^{(k)}\rangle$ (16)

where $\binom{d-1}{i}$ is the binomial coefficient which guarantees the
normalization, $|\mathbf{b}_{d-1}^{(k)}\rangle$ is given in Eq. (13) and
summation runs over all values of $k=0,\cdots,d-1$ in which
$|\mathbf{b}_{d-1}^{(k)}\rangle$ has only $j$ excitations. One can define the
transformation operator between the two bases as

$\mathcal{T}_{s}=\sum_{j=0}^{d-1}|\mathbf{f}_{d-1}^{(j)}\rangle\langle\widetilde{j}|.$
(17)

The symmetry encoding scheme effectively preserves certain types of
symmetries, for example the conservation of components of the total spin in a
spin system, i.e.
$[\widetilde{H},\widetilde{S}^{\alpha}_{tot}]=[H,S^{\alpha}_{tot}]=0$, where
$\widetilde{S}^{\alpha}_{tot}=\frac{1}{2}\sum_{i}\widetilde{S}^{\alpha}_{i},(\alpha{=}x,y,z)$
(Similar definition is given for $S^{\alpha}_{tot}$ in the qubit basis). Here,
$\widetilde{S}^{\alpha}_{i}$ is the $\alpha$-component of the spin operator
acting on the $i$-th site of the $d$-level system. This conservation implies
that the symmetry-preserving encoding also preserves the conservation of total
spin from $d$-level to qubit systems, i.e.
$[\widetilde{H},{\widetilde{S}_{\text{tot}}}^{2}]=[H,S_{\text{tot}}^{2}]=0$,
where
$\widetilde{S}_{\text{tot}}^{2}=(\widetilde{S}_{tot}^{x})^{2}+(\widetilde{S}_{tot}^{y})^{2}+(\widetilde{S}_{tot}^{z})^{2}$
(again similar definition is valid for $S_{\text{tot}}^{2}$ in the qubit
basis).

Although, in general, using the symmetry encoding demands extra qubits than
the binary one, it allows one to use ansatzes which have been developed for
preserving certain symmetries such as total spin or number of excitations [91,
97]. Those ansatzes are way more resource efficient for variational simulation
than the conventional hardware efficient ansatz which does not preserve
symmetries. Therefore, the final performance of the symmetric encoding in
comparison with the binary encoding relies on two opposite factors, namely the
use of extra qubits for the encoding versus the access to more resource
efficient circuits. Determining which of these factors is going to play a more
important role in variational simulation of $d$-level systems is the subject
of this paper.

It is important to note that, for both of the encoding methods, the Hilbert
space of the qubit system is larger than (or at least equal to) the one for
the qudit system. For instance, in a many-body system of $N$ qudit particles
the size of the Hilbert space is $d^{N}$. While using the binary encoding
enlarges the size of the Hilbert space to $2^{NM}$ and using symmetry encoding
enlarges it even more to $2^{N(d-1)}$. As a result, certain states within the
encoded qubit system are “illegitimate”, i.e. having no counterpart in the
qudit system. Thus, for any quantum state $|\phi\rangle$ within the
illegitimate subspace one inherently gets $H|\phi\rangle=0$. The illegitimate
states can indeed complicate the simulation process, particularly in VQE
implementations, where such states may erroneously appear as the state with
minimum energy. This issue typically arises when the actual ground state of
the $d$-level system has an average energy greater than zero, namely
$\langle\widetilde{gs}|\widetilde{H}|\widetilde{gs}\rangle>0$. Under these
conditions, the illegitimate states, which inherently possess zero energy,
will be the ground state of the encoded qubit system. To address this issue,
one can incorporate a penalty term into the cost function of the VQE that
specifically pushes the energy of these states upward.

To exemplify the inclusion of penalty term, consider a symmetry encoding on a
spin-$1$ system with $3$-level basis states

$|\widetilde{0}\rangle=|s_{z}{=}+1\rangle,\quad|\widetilde{1}\rangle=|s_{z}{=}0\rangle,\quad|\widetilde{2}\rangle=|s_{z}{=}-1\rangle,$
(18)

where $|s_{z}\rangle$ represents the $z$-component of the spin-1 particles.
The symmetry encoding results in

$|\widetilde{0}\rangle\mapsto|00\rangle,\quad|\widetilde{1}\rangle\mapsto\frac{1}{\sqrt{2}}(|01\rangle+|10\rangle),\quad|\widetilde{2}\rangle\mapsto|11\rangle.$
(19)

In this mapping, the singlet state $\frac{1}{\sqrt{2}}(|01\rangle-|10\rangle)$
is not used and thus considered as illegitimate. To effectively penalize this
state, the following penalty term is naturally constructed for the neighboring
pairs of qubits which encode one spin-1 particle:

$\displaystyle\mathcal{P}_{i,i+1}=|00\rangle\langle 00|$
$\displaystyle+\frac{1}{2}(|01\rangle{+}|10\rangle)(\langle 01|{+}\langle
10|)$ (20) $\displaystyle+|11\rangle\langle 11|$
$\displaystyle-\frac{1}{2}(|01\rangle{-}|10\rangle)(\langle 01|{-}\langle
10|).$

By considering the Eq. (III) one can simplify this penalty term into:

$\displaystyle\mathcal{P}_{i,i+1}=-(\sigma^{x}_{i}\otimes\sigma^{x}_{i+1}+\sigma^{y}_{i}\otimes\sigma^{y}_{i+1}+\sigma^{z}_{i}\otimes\sigma^{z}_{i+1}).$
(21)

So, the cost function of the VQE will be updated to

$\text{cost}=\langle H\rangle+\beta\sum_{\text{odd
}i}\langle\mathcal{P}_{i,i+1}\rangle,$ (22)

where $\beta$ is a sufficiently large positive scalar. By minimizing this
updated cost function the VQE naturally penalizes the illegitimate states
energetically and ensures that the output state of the qubit-based simulator
accurately reflects the ground state of the $d$-level system.

##  IV Variational Ground State Simulation

(a)

(b)

Figure 2: (a) The phase diagram of the spin-$1$ bilinear biquadratic
Heisenberg model. (b) The fidelity between the true qudit ground state
$|\widetilde{E}_{1}\rangle$ of the spin-$1$ BBH model with a system size of
$N{=6}$ and the ground state of the encoded $12$-qubit Hamiltonian
$|E_{1}\rangle$ across the phase diagram, namely ${|\langle
E_{1}|\mathcal{T}|\widetilde{E}_{1}\rangle|}^{2}$. Note that the result is
valid for both binary and symmetry encoding methods.

In this section, we present the VQE algorithm for simulating the ground state
of $d$-level systems using qubit-based quantum simulators.

###  IV.1 VQE simulation of spin-1 models

Specifically, we first focus on the simulation of a spin-$1$ system (i.e.
$d{=}3$). We consider the Bilinear-Biquadratic Heisenberg (BBH) model of size
$N$, defined as:

$\widetilde{H}_{\text{BBH}}=J\sum_{i=1}^{N-1}\left[\cos(\theta)\vec{\widetilde{S}}_{i}\cdot\vec{\widetilde{S}}_{i+1}+\sin(\theta)(\vec{\widetilde{S}}_{i}\cdot\vec{\widetilde{S}}_{i+1})^{2}\right],$
(23)

where $J$ is the coupling strength. In our ground state simulation, we put
$J{=}1$.
$\vec{\widetilde{S}}_{i}=(\widetilde{S}^{x}_{i},\widetilde{S}^{y}_{i},\widetilde{S}^{z}_{i})$
represents the spin-$1$ operators at the $i$-th site, with matrices given by:

$\displaystyle\widetilde{S}^{x}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{pmatrix}0&1&0\\\ 1&0&1\\\
0&1&0\end{pmatrix}$ (24) $\displaystyle\widetilde{S}^{y}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{pmatrix}0&-i&0\\\ i&0&-i\\\
0&i&0\end{pmatrix}$ $\displaystyle\widetilde{S}^{z}$
$\displaystyle=\begin{pmatrix}1&0&0\\\ 0&0&0\\\ 0&0&-1\end{pmatrix}.$

The BBH model is a generalization of the Heisenberg model that includes both
bilinear and biquadratic interactions between neighboring spins. The parameter
$\theta$ controls the relative strength of these interactions. The spin-$1$
BBH model exhibits several symmetries, including the conservation of the total
spin, expressed as
$[\widetilde{H}_{\text{BBH}},\widetilde{S}_{\text{tot}}^{2}]=0$. This
conservation law implies that every eigenstate $|\widetilde{E}_{k}\rangle$ of
the BBH model possesses a well-defined total spin $s$, which is an integer for
even system sizes $N$ and a half-integer for odd $N$. Consequently, the
expectation value of the total spin operator for any eigenstate is given by
$\langle\widetilde{E}_{k}|\widetilde{S}_{\text{tot}}^{2}|\widetilde{E}_{k}\rangle=s(s+1)$.
Additionally, the conservation of total spin implies the conservation of its
components in all directions, i.e.,
$[\widetilde{H}_{\text{BBH}},\widetilde{S}^{\alpha}_{tot}]=0$, which ensures
that the $z$-component of the total spin, $\widetilde{S}_{tot}^{z}$, for each
eigenstate is also conserved such that
$\langle\widetilde{E}_{k}|\widetilde{S}_{tot}^{z}|\widetilde{E}_{k}\rangle=s_{z}$,
with $s_{z}$ ranging from $-s$ to $s$.

The spin-$1$ BBH model is well-known in condensed matter physics for
exhibiting a rich phase diagram as $\theta$ varies from $-\pi$ to $\pi$. The
schematic of the phase diagram is shown in Fig. 2(a). In the range
$\theta\in\left[-\frac{3}{4}\pi,\frac{1}{2}\pi\right]$ the system is in a
ferromagnetic phase, where all spins are aligned in the same direction, making
the ground state trivial to simulate. We therefore focus our VQE simulations
on the non-trivial phases of the BBH model. For
$\theta\in\left[-\frac{3}{4}\pi,-\frac{1}{4}\pi\right]$, the system is in the
dimerized phase with broken translational invariance and singlets of
neighboring spins. The ground state in this phase transition can be
characterized by the order parameter [107]:

$\mathcal{O}_{\text{dimer}}=|\langle H_{i}-H_{i+1}\rangle|,$

where
$H_{i}=\cos(\theta)\vec{\widetilde{S}}_{i}\cdot\vec{\widetilde{S}}_{i+1}+\sin(\theta)(\vec{\widetilde{S}}_{i}\cdot\vec{\widetilde{S}}_{i+1})^{2}$.
In large systems or periodic boundary conditions, the choice of index $i$ does
not matter. Nonetheless, in finite open systems we select $i=2$ to avoid being
at the boundary. Notably, for
$\theta\in\left[-\frac{1}{4}\pi,\frac{1}{4}\pi\right]$, the system is in the
Haldane phase, which is a topological phase. Particularly, at the
$\theta=\text{arctan}(1/3)$ point, known as the Affleck-Kennedy-Lieb-Tasaki
(AKLT) point, the ground state is four-fold degenerate, including a singlet
state with total spin $s=0$ and a triplet state with total spin $s=1$. The
Haldane topological phase can be characterized by a string order parameter
[107]:

$\mathcal{O}_{\text{Haldane}}=\lim_{r\to\infty}\left\langle
S^{z}_{i}\exp\left[i\pi\sum_{j=i+1}^{i+r-1}S^{z}_{j}\right]S^{z}_{i+r}\right\rangle,$
(25)

with $r$ being the size of the string. Although, ideally $r$ has to be large
in our finite system it will be a finite number. In the region where
$\theta\in\left[\frac{1}{4}\pi,\frac{1}{2}\pi\right]$, the ground state is in
the critical phase, which is a gapless phase. The existence of this phase can
be reflected by a nematic structure factor at $q=2\pi/3$, which is defined as
[107]

$S(q)=\frac{1}{N}\sum_{k,l}e^{iq(k-l)}\langle{(S^{z}_{k})}^{2}{(S^{z}_{l})}^{2}\rangle.$

The first step towards the VQE simulation of the Hamiltonian (23) is to
generate its qubit representation. Interestingly, for spin-1 systems both
binary and symmetry encodings demand exactly $2$ qubits for representing each
spin. This provides an opportunity to compare the performance of these two
encodings using the same number of qubits. In order to map the Hamiltonian
into a qubit basis one has to use the transformation $H=\mathcal{T}^{\otimes
N}\widetilde{H}\mathcal{T}^{\dagger\otimes N}$, where $\mathcal{T}$ stands for
$\mathcal{T}_{b}$ (for binary encoding) and $\mathcal{T}_{s}$ for symmetric
encoding. Note that one has to pay extra attention to this transformation, in
particular, when illegitimate states are also involved. This is because
$\mathcal{T}\mathcal{T}^{\dagger}\neq\mathbb{I}_{\text{qubit}}$. Consequently,
the transformation has to apply to the interaction terms
$\vec{\widetilde{S}}_{i}\cdot\vec{\widetilde{S}}_{i+1}$ and
$(\vec{\widetilde{S}}_{i}\cdot\vec{\widetilde{S}}_{i+1})^{2}$, rather than
directly on the spin-$1$ operators $\widetilde{S}^{\alpha}$ (with
$\alpha=x,y,z$) and then use tensor product. The exact form of the
transformation for the symmetric encoding is in the following (for the binary
transformation, please see appendix A)

$\displaystyle\vec{\widetilde{S}}_{i}\cdot\vec{\widetilde{S}}_{i+1}\mapsto\frac{1}{4}\Big{(}\sigma_{j+1}^{x}\sigma_{j+3}^{x}+\sigma_{j+1}^{x}\sigma_{j+2}^{x}$
$\displaystyle\quad+\sigma_{j}^{x}\sigma_{j+3}^{x}+\sigma_{j}^{x}\sigma_{j+2}^{x}+\sigma_{j+1}^{y}\sigma_{j+3}^{y}+\sigma_{j+1}^{y}\sigma_{j+2}^{y}$
$\displaystyle\quad+\sigma_{j}^{y}\sigma_{j+3}^{y}+\sigma_{j}^{y}\sigma_{j+2}^{y}+\sigma_{j}^{z}\sigma_{j+2}^{z}+\sigma_{j}^{z}\sigma_{j+3}^{z}$
$\displaystyle\quad+\sigma_{j+1}^{z}\sigma_{j+2}^{z}+\sigma_{j+1}^{z}\sigma_{j+3}^{z}\Big{)}$
(26)
$\displaystyle(\vec{\widetilde{S}}_{i}\cdot\vec{\widetilde{S}}_{i+1})^{2}\mapsto\frac{3}{4}I+\frac{1}{4}\Big{(}\sigma_{j+2}^{x}\sigma_{j+3}^{x}$
$\displaystyle\quad+\sigma_{j}^{x}\sigma_{j+1}^{x}+\sigma_{j}^{x}\sigma_{j+1}^{x}\sigma_{j+2}^{x}\sigma_{j+3}^{x}\Big{)}$
$\displaystyle\quad-\frac{1}{8}\Big{(}\sigma_{j+1}^{z}\sigma_{j+3}^{z}+\sigma_{j+1}^{y}\sigma_{j+3}^{y}+\sigma_{j+1}^{z}\sigma_{j+2}^{z}$
$\displaystyle\quad+\sigma_{j+1}^{y}\sigma_{j+2}^{y}+\sigma_{j}^{z}\sigma_{j+3}^{z}+\sigma_{j}^{y}\sigma_{j+3}^{y}$
$\displaystyle\quad+\sigma_{j}^{z}\sigma_{j+2}^{z}+\sigma_{j}^{y}\sigma_{j+2}^{y}+\sigma_{j+1}^{x}\sigma_{j+3}^{x}$
$\displaystyle\quad+\sigma_{j+1}^{x}\sigma_{j+2}^{x}+\sigma_{j}^{x}\sigma_{j+3}^{x}+\sigma_{j}^{x}\sigma_{j+2}^{x}\Big{)}$
$\displaystyle\quad+\frac{1}{8}\Big{(}\sigma_{j}^{y}\sigma_{j+1}^{x}\sigma_{j+2}^{y}\sigma_{j+3}^{x}+\sigma_{j}^{z}\sigma_{j+1}^{x}\sigma_{j+2}^{z}\sigma_{j+3}^{x}$
$\displaystyle\quad+\sigma_{j}^{y}\sigma_{j+1}^{x}\sigma_{j+2}^{x}\sigma_{j+3}^{y}+\sigma_{j}^{z}\sigma_{j+1}^{x}\sigma_{j+2}^{x}\sigma_{j+3}^{z}$
$\displaystyle\quad+\sigma_{j}^{x}\sigma_{j+1}^{y}\sigma_{j+2}^{y}\sigma_{j+3}^{x}+\sigma_{j}^{x}\sigma_{j+1}^{z}\sigma_{j+2}^{z}\sigma_{j+3}^{x}$
$\displaystyle\quad+\sigma_{j}^{x}\sigma_{j+1}^{y}\sigma_{j+2}^{x}\sigma_{j+3}^{y}+\sigma_{j}^{x}\sigma_{j+1}^{z}\sigma_{j+2}^{x}\sigma_{j+3}^{z}$
$\displaystyle\quad+\sigma_{j}^{z}\sigma_{j+1}^{y}\sigma_{j+2}^{z}\sigma_{j+3}^{y}+\sigma_{j}^{z}\sigma_{j+1}^{y}\sigma_{j+2}^{y}\sigma_{j+3}^{z}$
$\displaystyle\quad+\sigma_{j}^{y}\sigma_{j+1}^{z}\sigma_{j+2}^{z}\sigma_{j+3}^{y}+\sigma_{j}^{y}\sigma_{j+1}^{z}\sigma_{j+2}^{y}\sigma_{j+3}^{z}\Big{)}$
$\displaystyle\quad+\frac{1}{4}\Big{(}\sigma_{j+2}^{y}\sigma_{j+3}^{y}+\sigma_{j}^{y}\sigma_{j+1}^{y}+\sigma_{j}^{y}\sigma_{j+1}^{y}\sigma_{j+2}^{y}\sigma_{j+3}^{y}$
$\displaystyle\quad+\sigma_{j+2}^{z}\sigma_{j+3}^{z}+\sigma_{j}^{z}\sigma_{j+1}^{z}+\sigma_{j}^{z}\sigma_{j+1}^{z}\sigma_{j+2}^{z}\sigma_{j+3}^{z}\Big{)}$
(27)

where $j=2i-1$ represents the corresponding index of the spin-$1$ operator in
the qubit system.

It is insightful to check whether the ground state of the qubit Hamiltonian
$|E_{1}\rangle$ truly represents the true qudit ground state
$|\widetilde{E}_{1}\rangle$ through the fidelity $|\langle
E_{1}|\mathcal{T}|\widetilde{E}_{1}\rangle|^{2}$. In Fig. 2(b), we present the
fidelity result $|\langle E_{1}|\mathcal{T}|\widetilde{E}_{1}\rangle|^{2}$ for
the spin-$1$ BBH model with a system size of $N{=}6$, which corresponds to a
$12$-qubit encoded Hamiltonian, across the phase diagram. Apparently, in the
critical phase and part of the Haldane phase, this fidelity drops to zero
showing that illegitimate states are indeed taking the low energy states. Note
that, the result is valid for both binary and symmetry encoding methods.
Therefore, without including penalty terms, see Eq. (22), one cannot target
the true ground state.

(a)

(b)

(c)

(d)

Figure 3: The VQE simulation of the ground state of the spin-$1$ BBH model
with system size $N{=}6$ for different values of $\theta$. The corresponding
qubit Hamiltonian uses $12$ qubits, encoded using the binary encoding method.
The results, averaged over 50 independent runs, include error bars
representing the standard deviation. The panels show the average obtainable
fidelity as a function of optimization iteration for the following settings:
(a) $\theta=-0.68\pi$ (dimerized), (b) $\theta=-0.55\pi$ (dimerized), (c)
$\theta=-0.30\pi$ (dimerized), and (d) $\theta=-0.16\pi$ (Haldane).

In order to have a comparison between the binary and the symmetry encodings we
focus on the VQE simulation of the ground state. In both cases to simulate a
many-body system of size $N$, we need $2N$ qubits. We start with the
simulation using the binary encoding method in which we use the hardware
efficient circuit as the variational ansatz. The results of the VQE
simulations for the BBH model with a system size of $N{=}6$, namely a $12$
qubit system after applying binary encoding, are presented in Fig. 3. The
average obtainable fidelity of the ground state with respect to the true
ground state of the BBH model is plotted as a function of the optimization
iterations for different values of $\theta$. All the simulations in this paper
use the gradient-based L-BFGS-B optimizer [102] for minimizing their cost
function. The results are averaged over 50 independent runs. It is clear from
the simulations that none of them has converged to the ground state of the BBH
model with a fidelity $\mathcal{F}\geq 0.85$. Additionally, we see a notable
performance drop as circuit depth increases beyond a certain point. This is
due to the appearance of the Barren plateaus problem in the optimization
landscape.

The VQE simulations with symmetry encoding utilize the total spin $s$
preserving ansatz, as illustrated in Fig. 1(d). In this case, it is crucial to
select an initial state that aligns with the symmetry preserving nature of the
ansatz. Specifically, the ground state of the BBH model exhibits a total spin
of $s{=}0$ across the dimerized, Haldane, and critical phases. Since the qubit
Hamiltonian consists of $2N$ qubits, the initial state is taken to be
$|\psi_{0}\rangle{=}\otimes^{N}|\psi^{-}\rangle$, where
$|\psi^{-}\rangle=(|01\rangle-|10\rangle)/\sqrt{2}$. This initial state
naturally has a total spin of $s{=}0$. As we will see, unlike the binary
encoding, the symmetry encoding method can faithfully simulate the ground
state in all phases. It is worth emphasizing that, the $|\psi^{-}\rangle$ is
an illegitimate state for the qubits, however, as
$|\psi_{0}\rangle{=}\otimes^{N}|\psi^{-}\rangle$ conserves the total spin it
eventually reaches the true ground state.

(a)

(b)

(c)

(d)

Figure 4: The VQE simulation of the ground state of the spin-$1$ BBH model
with system size $N{=}6$ for different values of $\theta$ in the dimerized
phase. The corresponding qubit Hamiltonian uses $12$ qubits, encoded using the
symmetry encoding method. The results, averaged over 50 independent runs,
include error bars representing the standard deviation. The upper panels
(a)-(b) show the average energy as a function of optimization iteration for
$\theta{=}{-}0.71\pi$ and $\theta{=}{-}0.30\pi$, respectively. The lower
panels (c)-(d) show the corresponding average fidelity for
$\theta{=}{-}0.71\pi$ and $\theta{=}{-}0.30\pi$.

Let’s consider the BBH model with system size $N=6$ at $\theta=-0.71\pi$ and
$\theta=-0.3\pi$, both within the dimerized phase. In this region, the ground
state of the qubit Hamiltonian represents the true ground state of the
original system. Therefore, the illegitimate states do not influence the
simulation outcomes. Consequently, the cost function was simply defined as the
expectation of the Hamiltonian, i.e. $\langle H\rangle$. The results of these
VQE simulations are depicted in Figs. 4(a)-(d), where the average energy and
obtainable fidelity are plotted as a function of optimization iterations for
different ansatz layers. These results, whose error bars are computed over 50
independent runs, demonstrate that in both cases, the VQE simulations
successfully converged to the ground state of the BBH model with a fidelity
$\mathcal{F}\geq 0.95$ after $\sim 1400$ optimization iterations, employing
$6$ and $5$ layers of the total spin-preserving ansatz, for $\theta=-0.71\pi$
and $\theta=-0.3\pi$ respectively.

(a)

(b)

(c)

(d)

Figure 5: The VQE simulation of the ground state of the spin-$1$ BBH model
with system size $N{=}6$ for different values of $\theta$ in the Haldane
phase. The corresponding qubit Hamiltonian uses $12$ qubits, encoded using the
symmetry encoding method. The results, averaged over 50 independent runs,
include error bars representing the standard deviation. The upper panels
(a)-(b) show the average energy as a function of optimization iteration for
$\theta{=}{-}0.16\pi$ and $\theta{=}\arctan(1/3)$ (the AKLT point),
respectively. The lower panels (c)-(d) show the corresponding average fidelity
for $\theta{=}{-}0.16\pi$ and $\theta{=}\arctan(1/3)$.

Subsequently, we explore the Haldane phase of the BBH model, particularly at
$\theta=-0.16\pi$ and $\theta=\text{arctan}(1/3)$. Similar to the dimerized
phase, in this phase also the ground state of the qubit Hamiltonian faithfully
represents the true original Hamiltonian of the system. In other words, the
illegitimate states become high-energy eigenstates. At the
$\theta=\text{arctan}(1/3)$ point, known as the Affleck-Kennedy-Lieb-Tasaki
(AKLT) point, the ground state is four-fold degenerate, including a singlet
state with total spin $s=0$ and a triplet state with total spin $s=1$. The
AKLT point has a matrix product representation which makes its simulation
accessible to classical computers. By using a total spin preserving ansatz,
one can target any of these degenerate states by a proper choice of the
initial state with a given total spin $s$. For instance, for simulating the
ground state with total spin $s{=}0$ one can use the initial state
$|\psi_{0}\rangle=\otimes^{N}|\psi^{-}\rangle$. On the other hand, for
simulating the other ground states with $s{=}1$ one can replace one of the
singlet states $|\psi^{-}\rangle$ by
$|\psi^{+}\rangle{=}(|01\rangle+|10\rangle)/\sqrt{2}$ (for $s_{z}=0$) or
$|\phi^{\pm}\rangle{=}(|00\rangle\pm|11\rangle)/\sqrt{2}$ (for $s_{z}=\pm 1$).
The VQE simulations with different numbers of ansatz layers are presented in
Figs. 5(a)-(d). The error bars are computed for 50 different runs. Both
simulations require 5 layers of the total spin-preserving ansatz to converge
to the ground state of the BBH model with fidelity $\mathcal{F}\geq 0.95$. The
convergence requires approximately $1750$ and $1400$ optimization iterations
for $\theta=-0.16\pi$ and $\theta=\text{arctan}(1/3)$, respectively. These
results demonstrate the capability of the VQE algorithm to accurately simulate
the Haldane phase.

(a)

(b)

Figure 6: The VQE simulation of the ground state of the spin-$1$ BBH model
with system size $N{=}6$ for the critical phase. The corresponding qubit
Hamiltonian uses $12$ qubits, encoded using the symmetry encoding method. The
results, averaged over 50 independent runs, include error bars representing
the standard deviation. The average energy for $\theta{=}0.32\pi$, shown in
panel (a), and the corresponding obtainable fidelity, shown in panel (b), is
plotted as a function of optimization iteration.

We further explore the critical phase of the BBH model at $\theta=0.32\pi$, a
regime where illegitimate states significantly influence the simulation
outcomes as they are becoming the low energy eigenstates of the system. In
fact, for a system of $N=6$, the real ground state of the original Hamiltonian
corresponds to the eigenstate number $1228$ in the qubit Hamiltonian. In order
to target the true ground state, we integrate a penalty term, as specified in
Eqs. (21)–(22), into the VQE cost function with $\beta=10$. The VQE results
for various depths of the circuit, averaged over 50 independent trials, are
depicted in Figs. 6(a)-(b). These results demonstrate that the simulations
consistently converged to the ground state of the BBH model with a fidelity
$\mathcal{F}\geq 0.95$ using a quantum circuit with 6 layers. Inevitably, the
penalty term increases the computational cost of the VQE simulations, as
indicated by the higher number of optimization iterations compared to other
phases. Nonetheless, the simulations successfully capture the ground state of
the BBH model in the critical phase, even though the true ground state of the
original Hamiltonian is a high-energy eigenstate of the qubit Hamiltonian.

(a)

(b)

(c)

Figure 7: The value of order parameters for the dimerized, Haldane, and
critical phases of the BBH model, calculated using the output state of the VQE
simulation, are shown in panels (a), (b), and (c), respectively. The results
are averaged over 50 independent runs. As the number of total spin preserving
ansatz layers increases, the obtained order parameters more closely align with
the corresponding theoretical values.

As discussed before, each phase of the BBH model is characterized by an order
parameter. One can measure the order parameters at the output of the
simulator. The results of these order parameters are presented in Fig. 7,
which confirm that the VQE simulations accurately reproduce the ground state
properties of the BBH model across the dimerized, Haldane, and critical
phases.

Binary encoding, $N=12$
---
| $\theta=-0.75\pi$ | $\theta=-0.36\pi$ | $\theta=-0.15\pi$
#CNOT | 132 | 132 | 132
#Para. | 312 | 312 | 312
$\mathcal{F}$ |  0.432 $\pm$ 0.417 |  0.806 $\pm$ 0.208 |  0.513 $\pm$ 0.256
Iter. |  3076.4 $\pm$ 1082.7 |  2785.3 $\pm$ 1201.8 |  3244.2 $\pm$ 1211.9
$C_{R}$ ($\times 10^{3}$) |  959.8 $\pm$ 337.8 |  869.0 $\pm$ 375.0 |  1012.2 $\pm$ 378.1
Symmetry encoding, $N=18$
| $\theta=-0.75\pi$ | $\theta=-0.36\pi$ | $\theta=-0.15\pi$
#CNOT | 612 | 459 | 612
#Para. | 204 | 153 | 204
$\mathcal{F}$ |  0.970 $\pm$ 0.008 |  0.951 $\pm$ 0.010 |  0.973 $\pm$ 0.011
Iter. |  8806 $\pm$ 1839.7 |  4888.7 $\pm$ 1479.5 |  6958.2 $\pm$ 1712.9
$C_{R}$ ($\times 10^{3}$) |  1796.4 $\pm$ 375.3 |  748.0 $\pm$ 226.4 |  1419.5 $\pm$ 350.4
Table 1: Comparison between binary and symmetry encoding methods at ($\theta$)
for $N=12$ and $N=18$. The table shows the number of CNOT gates, the number of
parameters, fidelity ($\mathcal{F}$), iteration count, and classical resources
($C_{R}$), with standard deviations for fidelity, iteration and classical
resource metrics.

###  IV.2 VQE simulation of spin-3/2 models

For spin-$1$ systems, where both binary and symmetry encodings require the
same number of qubits, the result for both cases shows a clear advantage on
the symmetry encoding. This advantage is expected as exploiting symmetries in
VQE simulations reduces the search subspace in the Hilbert subspace to only
the relevant part in which all the states have the desired symmetries. This
matches with the previous results for VQE simulation of qubit systems [97]. By
increasing $d$, however, the binary and symmetry encoding methods require
different numbers of qubits. In this situation, it is not clear whether the
symmetry encoding still performs better as it demands more qubits for the
mapping. To explore this situation, we focus on a spin-$\frac{3}{2}$ model. We
consider the same Hamiltonian as the spin-$1$ BBH model, shown in Eq. (23),
but with the spin-1 operators replaced by spin-$\frac{3}{2}$ operators. The
corresponding spin-$3/2$ operators are represented by matrices:

$\displaystyle\widetilde{S}^{x}$
$\displaystyle=\frac{1}{2}\begin{pmatrix}0&\sqrt{3}&0&0\\\ \sqrt{3}&0&2&0\\\
0&2&0&\sqrt{3}\\\ 0&0&\sqrt{3}&0\end{pmatrix},$
$\displaystyle\widetilde{S}^{y}$
$\displaystyle=\frac{1}{2i}\begin{pmatrix}0&\sqrt{3}&0&0\\\ -\sqrt{3}&0&2&0\\\
0&-2&0&\sqrt{3}\\\ 0&0&-\sqrt{3}&0\end{pmatrix},$
$\displaystyle\widetilde{S}^{z}$
$\displaystyle=\frac{1}{2}\begin{pmatrix}3&0&0&0\\\ 0&1&0&0\\\ 0&0&-1&0\\\
0&0&0&-3\end{pmatrix}.$

The spin-$3/2$ BBH model is a 4-level system with the basis determined by the
eigenvalues of the $S_{z}$ operator. The basis is then given by

$\displaystyle|\widetilde{0}\rangle$ $\displaystyle=$
$\displaystyle|s_{z}{=}{+}\frac{3}{2}\rangle,\quad|\widetilde{1}\rangle=|s_{z}{=}{+}\frac{1}{2}\rangle,$
(28) $\displaystyle|\widetilde{2}\rangle$ $\displaystyle=$
$\displaystyle|s_{z}{=}{-}\frac{1}{2}\rangle,\quad|\widetilde{3}\rangle=|s_{z}{=}{-}\frac{3}{2}\rangle.$
(29)

To encode these states into qubits using binary encoding one requires $2$
qubits as:

$\displaystyle|\widetilde{0}\rangle\mapsto|00\rangle,\quad|\widetilde{1}\rangle\mapsto|01\rangle,$
(30)
$\displaystyle|\widetilde{2}\rangle\mapsto|10\rangle,\quad|\widetilde{3}\rangle\mapsto|11\rangle.$

On the other hand, the symmetry encoding method of the spin-$3/2$ BBH model,
requires $3$ qubits for each spin. The mapping for the symmetry encoding
method is as follows:

$\displaystyle|\widetilde{0}\rangle\mapsto|000\rangle,\quad|\widetilde{1}\rangle\mapsto\frac{1}{\sqrt{3}}(|001\rangle+|010\rangle+|100\rangle),$
(31)
$\displaystyle|\widetilde{2}\rangle\mapsto\frac{1}{\sqrt{3}}(|110\rangle+|101\rangle+|011\rangle),\quad|\widetilde{3}\rangle\mapsto|111\rangle.$

Thus, exploiting the symmetry in the spin-$3/2$ BBH model comes at the cost of
using more qubits with the symmetry encoding method. Here, we compare the
performance of binary encoding and symmetry encoding methods for VQE
simulation of the spin-$3/2$ BBH model. In Table 1, we present VQE simulation
results for this model with $N=6$ spins, which is mapped to qubit systems of
size $2N$ and $3N$ for binary and symmetry encoding methods, respectively. The
results are averaged over 50 independent runs. The table shows the number of
CNOT gates (i.e. quantum resources), the number of parameters to optimize,
obtainable fidelity, optimization iterations and classical resources, defined
in Eq. (9), for three different values of $\theta$. The results using binary
encoding with a hardware-efficient ansatz are shown in the top half of the
table, which fails to converge to fidelity $\mathcal{F}\geq 0.95$ for all
values of $\theta$. More importantly, adding more layers to the ansatz
decreases the fidelity, indicating the presence of a Barren plateau problem in
the VQE simulations. In contrast, the results using symmetry encoding with a
total spin preserving ansatz, shown in the bottom half of the table,
successfully converge to fidelity $\mathcal{F}\geq 0.95$ for all tested
$\theta$ values. The results demonstrate that the symmetry encoding
significantly outperforms the binary encoding despite the overhead of extra
qubits. Note that these results are presented without implementing any Barren
plateau mitigation techniques, which could potentially enhance both the
performance of the binary and the symmetry encoding methods.

| $U=-2$ | $U=2$
---|---|---
#CNOT | 330 | 297
#Para. | 230 | 207
$\mathcal{F}$ |  0.975$\pm$0.004 |  0.999$\pm$0.001
Iter. |  11802.9$\pm$2669.6 |  7538.8$\pm$2183.5
$C_{R}$ ($\times 10^{3}$) |  2714.7$\pm$614.0 |  1560.5$\pm$452.0
Table 2: Simulation results for the Bose-Hubbard model with different
interaction strengths ($U$). This table shows the number of CNOT gates, number
of parameters, fidelity ($\mathcal{F}$), iteration count, and classical
resources ($C_{R}$), with respective standard deviations, across two
interaction strengths: $U=-2$, and $U=2$.

###  IV.3 VQE simulation of bosonic systems

Finally, we turn our attention to the simulation of bosonic systems in which
$d\rightarrow\infty$. Specifically, we focus on the Bose-Hubbard model defined
by the Hamiltonian:

$\widetilde{H}_{\text{BH}}=-t\sum_{i=1}^{N-1}(a_{i}^{\dagger}a_{i+1}+a_{i+1}^{\dagger}a_{i})+\frac{U}{2}\sum
n_{i}(n_{i}-1),$ (32)

where $N$ is the size of the lattice, $a_{i}^{\dagger}$ and $a_{i}$ are the
bosonic creation and annihilation operators at the $i$-th site, respectively,
and $n_{i}=a_{i}^{\dagger}a_{i}$ is the number operator. The parameter $t$
denotes the hopping strength, and $U$ represents the on-site interaction. Note
that the Bose-Hubbard Hamiltonian has a symmetry of preserving the total
number of bosons as $[\widetilde{H}_{\text{BH}},\mathcal{W}]=0$, where
$\mathcal{W}=\sum_{i}n_{i}$ is the total number operator. Theoretically, the
Bose-Hubbard model is infinite-dimensional, posing significant challenges for
direct simulation on a qubit-based quantum simulator. However, by truncating
the bosonic Hilbert space to a finite dimension, effective simulation becomes
feasible. In fact, for a given number of bosons $n_{bos}$, namely
$\langle\mathcal{W}\rangle=n_{bos}$, the Hilbert space of each site is
naturally limited to $d{=}n_{bos}+1$. Therefore, if the truncation number is
below this value it may result in some errors. Nonetheless, the truncation
error may still be negligible in the limit of $U\gg 0$, where repulsive
interaction between bosonic excitations makes the occupation number of each
site as low as possible. In contrast, negative $U$, i.e. attractive
interaction, tends to accumulate all the excitations on one site which
naturally increases the Hilbert space of each site and thus enhances the
truncation error. For our simulations, we adopt a truncation limit of four
levels per site (i.e. $d{=}4$) represented by the bosonic Fock states
$\\{|\widetilde{0}\rangle,|\widetilde{1}\rangle,|\widetilde{2}\rangle,|\widetilde{3}\rangle\\}$.
Similar to the case of spin-$3/2$, the binary encoding with hardware efficient
ansatz cannot converge to good results. Therefore, we only focus on the
symmetry encoding approach. The encoding is as follows:

$\displaystyle|\widetilde{0}\rangle$ $\displaystyle\mapsto|000\rangle,$ (33)
$\displaystyle|\widetilde{1}\rangle$
$\displaystyle\mapsto\frac{1}{\sqrt{3}}(|100\rangle+|010\rangle+|001\rangle),$
$\displaystyle|\widetilde{2}\rangle$
$\displaystyle\mapsto\frac{1}{\sqrt{3}}(|110\rangle+|101\rangle+|011\rangle),$
$\displaystyle|\widetilde{3}\rangle$ $\displaystyle\mapsto|111\rangle.$

Let’s consider a lattice of $N{=}4$ sites which is mapped to a system of $3N$
qubits through
$H_{\text{BH}}=\mathcal{T}_{s}\widetilde{H}_{\text{BH}}\mathcal{T}_{s}^{\dagger}$.
In other words, the 4-site Bose-Hubbard model is mapped to a 12-qubit system.
We also assume that there are only $n_{bos}{=}3$ bosons in the systems. The
symmetry encoding method effectively preserves the conservation of the total
number of bosons. This conservation translates to the conservation of
excitations in the qubit system, which also implies the conservation of $z$
component of the total spin, denoted by $[H_{\text{BH}},S^{z}_{tot}]=0$. Thus,
we utilize the $s_{z}$-conserving ansatz, as illustrated in Fig. 1(c).

Regarding the on-site interaction strength $U$, we investigate two distinct
regimes, namely $U=\pm 2$. Ground state simulations with positive $U$ values
are relatively trivial as bosons naturally avoid clustering, resulting in the
vacuum state as the ground state. Therefore, our simulations concentrate on
the subspace with three bosons, i.e. $n_{bos}{=}3$, across four sites, i.e.
$N{=}4$. The initial state of the VQE circuit is taken to be
$|\psi_{0}\rangle=|\widetilde{3},\widetilde{0},\widetilde{0},\widetilde{0}\rangle$
whose qubit representation is $|\psi_{0}\rangle=|111000000000\rangle$. With
negative $U$, bosons tend to gather, challenging the limit of our truncated
Hilbert space and complicating the VQE simulations. The results of these
simulations are summarized in Table 2, showing the resources required for the
VQE to achieve a fidelity of $\mathcal{F}\geq 0.95$, averaged over 50
independent runs. The data clearly indicates that as $U$ transitions from
negative to positive, the computational demands for the VQE simulations
decrease, implying easier simulation. These results show that the VQE
simulations robustly converge to the ground state of the Bose-Hubbard model
across all explored regimes, demonstrating the effectiveness of the VQE
algorithm in accurately simulating bosonic systems.

(a)

(b)

Figure 8: The variational time evolution simulation of the spin-$1$ BBH model
with system size $N{=}2$ over time $\tau{=}6/J$ for: (a)
$\theta{=}\frac{\pi}{4}$, where $\frac{\cos(\theta)}{\sin(\theta)}=1$; and (b)
$\theta{=}\arctan(\frac{1}{2})$, where
$\frac{\cos(\theta)}{\sin(\theta)}=0.5$. The corresponding qubit Hamiltonian
uses $4$ qubits, encoded using the symmetry encoding method. The average
fidelity of the variational time evolution results states with two states
$|\widetilde{0}{,}\widetilde{2}\rangle$, and
$|\widetilde{2}{,}\widetilde{0}\rangle$ is plotted as a function of the
evolution time, labeled VTE in the figure. The theoretical result of the time
evolution of the system is shown in dash lines.

##  V Variational Time Evolution Simulation

In this section, we extend our study to include the variational time evolution
of $d$-level systems using qubit-based quantum simulators. We specifically
demonstrate the time evolution simulation of a spin-$1$ BBH model (i.e.
$d{=}3$) shown in Eq. (23). Notably, this Hamiltonian conserves the
$z$-component of the total spin, namely
($[\widetilde{H}_{\text{BBH}},S^{z}_{tot}]=0$). Hence, for example, if the
initial state of the system is chosen to be
$|\widetilde{0},\widetilde{2}\rangle$, as the time evolution proceeds, the
state of the system can only be a superposition of states with the same total
spin $s_{z}$, namely $|\widetilde{0},\widetilde{2}\rangle$,
$|\widetilde{2},\widetilde{0}\rangle$, and
$|\widetilde{1},\widetilde{1}\rangle$. Furthermore, for the case of
$\theta{=}\frac{\pi}{4}$, the Hamiltonian also conserve
$\frac{1}{2}\sum_{i}(S^{z}_{i})^{2}$. This implies that for $\theta=\pi/4$,
the quantum state $|\widetilde{0},\widetilde{2}\rangle$ remains in the
subspace spanned by $|\widetilde{0},\widetilde{2}\rangle$, and
$|\widetilde{2},\widetilde{0}\rangle$. In fact, the time evolution will act
similar to a swap operation between the two states
$|\widetilde{0},\widetilde{2}\rangle$ and
$|\widetilde{2},\widetilde{0}\rangle$.

In order to perform variational time evolution, we consider a system of $N=2$
spin-1 particles interacting under the BBH Hamiltonian with two different
values of $\theta$. Especially, $\theta=\pi/4$
($\frac{\cos(\theta)}{\sin(\theta)}=1$), where the time evolution acts similar
to a state swapping process, and $\theta=\arctan(1/2)$
($\frac{\cos(\theta)}{\sin(\theta)}=0.5$), to demonstrate a typical non-
equilibrium time evolution. We employ the symmetry encoding (19) to evolve the
system up to time $\tau{=}6/J$, which requires $4$ qubits. The initial state
of the system is chosen to be
$|\Psi(0)\rangle=|\widetilde{0},\widetilde{2}\rangle$ which will be mapped to
$|1100\rangle$ for the qubit simulation. The time evolution is given by
$|\Psi(\tau)\rangle=\exp{\left(-iH_{\text{BBH}}\tau\right)}|\Psi(0)\rangle$.
We utilize the total spin preserving ansatz, as illustrated in Fig. 1(d). The
variational time evolution simulation results are presented in Figs. 8(a)-(b),
where two fidelities, namely
$|\langle\widetilde{0},\widetilde{2}|\Psi(\tau)\rangle|^{2}$ and
$|\langle\widetilde{2},\widetilde{0}|\Psi(\tau)\rangle|^{2}$, are plotted as a
function of time $J\tau$. One can clearly see a good match between the exact
results and the results from the variational time evolution. This demonstrates
the capability of symmetry encoding in capturing the properties of the non-
equilibrium dynamics of $d$-level systems through a variational time evolution
algorithm.

##  VI Conclusion

In this paper, we investigate variational quantum simulation of $d$-level
systems on qubit-based quantum simulators. In particular, we consider higher
spins and bosonic systems. Both ground state and non-equilibrium dynamics have
been explored. We have compared two different types of encoding methods for
such simulations, namely binary and symmetry encodings. While binary encoding
is efficient with respect to the number of qubits, requiring only
$\lceil\log_{2}(d)\rceil$ qubits per particle, it may not be able to
incorporate the symmetries in the design of the quantum circuit. On the other
hand, symmetry encoding requires $d{-}1$ qubits per particle to ensure
preserving of typical symmetries, such as $S_{tot}^{2}$ and $S_{tot}^{z}$. In
both of these methods, the Hilbert space of the qubit system is typically
larger than the Hilbert space of the original $d$-level system. This naturally
introduces illegitimate states in the Hilbert space of the qubit system which
has no correspondence in the real $d$-level system. By mapping the $d$-level
Hamiltonian into qubit systems, sometimes the ground state of the qubit system
lies in the illegitimate subspace. We have systematically introduced a penalty
term which pushes the illegitimate states into high-energy states to guarantee
the success of the VQE simulations. Interestingly, through simulation of a
wide range of phases in spin-$1$, spin-$3/2$, as well as the Bose-Hubbard
model, we observe that the binary encoding is not scalable as the hardware
efficient ansatz cannot be trained due to the presence of Barren plateau. On
the other hand, symmetry encoding, despite requiring extra qubits, does not
face such trainability issue thanks to exploiting symmetries. Our proposal is
very timely as it can be implemented on existing quantum simulators and can
easily be generalized to cover a wide range of qudit models.

##  VII Acknowledgement

A.B. acknowledges support from the National Natural Science Foundation of
China (Grants No. 12050410253, No. 92065115, and No. 12274059), and the
Ministry of Science and Technology of China (Grant No. QNJ2021167001L)

## References

  * Preskill [2018] J. Preskill, Quantum 2, 79 (2018).
  * Feynman [1982] R. P. Feynman, International Journal of Theoretical Physics 21, 467 (1982).
  * Lloyd [2018] S. Lloyd, Quantum approximate optimization is computationally universal (2018), arXiv:1812.11075 [quant-ph] .
  * Bordia  _et al._ [2017] P. Bordia, H. Lüschen, S. Scherg, _et al._ , Phys. Rev. X 7, 041047 (2017).
  * Schreiber  _et al._ [2015] M. Schreiber, S. S. Hodgman, P. Bordia, _et al._ , Science 349, 842 (2015).
  * Gross and Bloch [2017] C. Gross and I. Bloch, Science 357, 995 (2017).
  * Omran  _et al._ [2019] A. Omran, H. Levine, A. Keesling, _et al._ , Science 365, 570 (2019).
  * Keesling  _et al._ [2019] A. Keesling, A. Omran, H. Levine, _et al._ , Nature 568, 207 (2019).
  * Wang  _et al._ [2017] J. Wang, S. Paesani, R. Santagati, _et al._ , Nat. Phys. 13, 551 (2017).
  * Zhong  _et al._ [2020] H.-S. Zhong, H. Wang, Y.-H. Deng, _et al._ , Science 370, 1460 (2020).
  * Carolan  _et al._ [2020] J. Carolan, M. Mohseni, J. P. Olson, _et al._ , Nat. Phys. 16, 322 (2020).
  * Li  _et al._ [2017] J. Li, R. Fan, H. Wang, _et al._ , Phys. Rev. X 7, 031011 (2017).
  * Lanyon  _et al._ [2011] B. P. Lanyon, C. Hempel, D. Nigg, _et al._ , Science 334, 57 (2011).
  * Zhang  _et al._ [2017] J. Zhang, G. Pagano, P. W. Hess, _et al._ , Nature 551, 601 (2017).
  * Hempel  _et al._ [2018] C. Hempel, C. Maier, J. Romero, _et al._ , Phys. Rev. X. 8, 031022 (2018).
  * Kokail  _et al._ [2019] C. Kokail, C. Maier, R. van Bijnen, _et al._ , Nature 569, 355 (2019).
  * Han  _et al._ [2024] Z. Han, C. Lyu, Y. Zhou, _et al._ , Phys. Rev. Res. 6, 013015 (2024).
  * Salathé  _et al._ [2015] Y. Salathé, M. Mondal, M. Oppliger, _et al._ , Phys. Rev. X. 5, 021027 (2015).
  * Wang  _et al._ [2020] C. S. Wang, J. C. Curtis, B. J. Lester, _et al._ , Phys. Rev. X 10, 021060 (2020).
  * Karamlou  _et al._ [2021] A. H. Karamlou, W. A. Simon, A. Katabarwa, _et al._ , npj Quantum Inf. 7, 156 (2021).
  * Neill  _et al._ [2021] C. Neill, T. McCourt, X. Mi, _et al._ , Nature 594, 508 (2021).
  * Han  _et al._ [2021] J. Han, W. Cai, L. Hu, _et al._ , Phys. Rev. Lett. 127, 020504 (2021).
  * Braumüller  _et al._ [2022] J. Braumüller, A. H. Karamlou, Y. Yanay, _et al._ , Nat. Phys. 18, 172 (2022).
  * Zhang  _et al._ [2022] X. Zhang, W. Jiang, J. Deng, _et al._ , Nature 607, 468–473 (2022).
  * Shi  _et al._ [2023] Y.-H. Shi, Y. Liu, Y.-R. Zhang, _et al._ , Phys. Rev. Lett. 131, 080401 (2023).
  * Levitt [2013] M. H. Levitt,  _Spin dynamics: basics of nuclear magnetic resonance_ (John Wiley & Sons, 2013).
  * Fisher  _et al._ [1989] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
  * Wilson  _et al._ [1955] E. B. Wilson, J. C. Decius, P. C. Cross, and B. R. Sundheim, J. Electrochem. Soc. 102, 235Ca (1955).
  * Shimizu [1981] M. Shimizu, Rep. Prog. Phys. 44, 329 (1981).
  * Glos  _et al._ [2022] A. Glos, A. Krawiec, and Z. Zimborás, npj Quantum Inf. 8, 39 (2022).
  * Tabi  _et al._ [2020] Z. Tabi, K. H. El-Safty, Z. Kallus, _et al._ , in  _2020 IEEE International Conference on Quantum Computing and Engineering (QCE)_ (2020) pp. 56–62.
  * Liu  _et al._ [2023] P. Liu, R. Wang, J.-N. Zhang, _et al._ , Phys. Rev. X 13, 021028 (2023).
  * Meth  _et al._ [2024] M. Meth, J. F. Haase, J. Zhang, _et al._ , Simulating 2d lattice gauge theories on a qudit quantum computer (2024), arXiv:2310.12110 [quant-ph] .
  * Zache  _et al._ [2023] T. V. Zache, D. González-Cuadra, and P. Zoller, Quantum 7, 1140 (2023).
  * González-Cuadra  _et al._ [2022] D. González-Cuadra, T. V. Zache, J. Carrasco, _et al._ , Phys. Rev. Lett. 129, 160501 (2022).
  * Hrmo  _et al._ [2023] P. Hrmo, B. Wilhelm, L. Gerster, _et al._ , Nat. Commun. 14, 2242 (2023).
  * Ringbauer  _et al._ [2022] M. Ringbauer, M. Meth, L. Postler, _et al._ , Nat. Phys. 18, 1053 (2022).
  * Luo and Wang [2014] M. Luo and X. Wang, Sci. China Phys. Mech. Astron. 57, 1712 (2014).
  * Meth  _et al._ [2023] M. Meth, J. F. Haase, J. Zhang, _et al._ , arXiv:2310.12110 (2023).
  * González-Cuadra  _et al._ [2022] D. González-Cuadra, T. V. Zache, J. Carrasco, _et al._ , Physical Review Letters 129, 160501 (2022).
  * Wu and Lidar [2002] L.-A. Wu and D. A. Lidar, J. Math. Phys. 43, 4506 (2002).
  * Batista and Ortiz [2004] C. D. Batista and G. Ortiz, Adv. Phys. 53, 1 (2004).
  * Sawaya  _et al._ [2020] N. P. D. Sawaya, T. Menke, T. H. Kyaw, _et al._ , npj Quantum Inf. 6, 49 (2020).
  * Dicke [1954] R. H. Dicke, Phys. Rev. 93, 99 (1954).
  * Yang  _et al._ [2024] K. Yang, Y. Zhu, X. Zeng, _et al._ , Cost of locally approximating high-dimensional ground states of contextual quantum models (2024), arXiv:2405.04884 [quant-ph] .
  * Zhong  _et al._ [2021] H.-S. Zhong, Y.-H. Deng, J. Qin, _et al._ , Phys. Rev. Lett. 127, 180502 (2021).
  * Deng  _et al._ [2023] Y.-H. Deng, Y.-C. Gu, H.-L. Liu, _et al._ , Phys. Rev. Lett. 131, 150601 (2023).
  * Arute  _et al._ [2020] F. Arute, K. Arya, R. Babbush, _et al._ , Science 369, 1084–1089 (2020).
  * O’Malley  _et al._ [2016] P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, _et al._ , Phys. Rev. X. 6, 031007 (2016).
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, _et al._ , Nature 549, 242 (2017).
  * Colless  _et al._ [2018] J. I. Colless, V. V. Ramasesh, D. Dahlen, _et al._ , Phys. Rev. X. 8, 011021 (2018).
  * Nam  _et al._ [2020] Y. Nam, J.-S. Chen, N. C. Pisenti, _et al._ , npj Quantum Inf. 6, 33 (2020).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm (2014), arXiv:1411.4028 [quant-ph] .
  * Bravyi  _et al._ [2020] S. Bravyi, A. Kliesch, R. Koenig, and E. Tang, Phys. Rev. Lett. 125, 260505 (2020).
  * Karácsony  _et al._ [2024] M. Karácsony, L. Oroszlány, and Z. Zimboras, SciPost Phys. Core 7, 032 (2024).
  * Biamonte  _et al._ [2017] J. Biamonte, P. Wittek, N. Pancotti, _et al._ , Nature 549, 195 (2017).
  * Arunachalam and de Wolf [2017] S. Arunachalam and R. de Wolf, A survey of quantum learning theory (2017), arXiv:1701.06806 [quant-ph] .
  * Ciliberto  _et al._ [2018] C. Ciliberto, M. Herbster, A. D. Ialongo, _et al._ , Proc. R. Soc. A: Math. Phys. Eng. Sci. 474, 20170551 (2018).
  * Dunjko and Briegel [2018] V. Dunjko and H. J. Briegel, Reports on Progress in Physics 81, 074001 (2018).
  * Farhi and Neven [2018] E. Farhi and H. Neven, Classification with quantum neural networks on near term processors (2018), arXiv:1802.06002 [quant-ph] .
  * Schuld and Killoran [2019] M. Schuld and N. Killoran, Phys. Rev. Lett. 122, 040504 (2019).
  * Wilkinson and Hartmann [2022] S. A. Wilkinson and M. J. Hartmann, Evaluating the performance of sigmoid quantum perceptrons in quantum neural networks (2022), arXiv:2208.06198 [quant-ph] .
  * Zapletal  _et al._ [2023] P. Zapletal, N. A. McMahon, and M. J. Hartmann, Error-tolerant quantum convolutional neural networks for symmetry-protected topological phases (2023), arXiv:2307.03711 [quant-ph] .
  * Cirstoiu  _et al._ [2020] C. Cirstoiu, Z. Holmes, J. Iosue, _et al._ , npj Quantum Inf. 6, 1 (2020).
  * Gibbs  _et al._ [2021] J. Gibbs, K. Gili, Z. Holmes, _et al._ , Long-time simulations with high fidelity on quantum hardware (2021), arXiv:2102.04313 [quant-ph] .
  * Yuan  _et al._ [2019] X. Yuan, S. Endo, Q. Zhao, _et al._ , Quantum 3, 191 (2019).
  * McArdle  _et al._ [2019] S. McArdle, T. Jones, S. Endo, _et al._ , npj Quantum Inf. 5, 1 (2019).
  * Heya  _et al._ [2023] K. Heya, K. M. Nakanishi, K. Mitarai, _et al._ , Subspace variational quantum simulator (2023), arXiv:1904.08566 [quant-ph] .
  * Huh  _et al._ [2014] J. Huh, S. Mostame, T. Fujita, _et al._ , New J. Phys. 16, 123008 (2014).
  * Hu  _et al._ [2020] Z. Hu, R. Xia, and S. Kais, Sci. Rep. 10, 1 (2020).
  * Endo  _et al._ [2020] S. Endo, J. Sun, Y. Li, _et al._ , Phys. Rev. Lett. 125, 010501 (2020).
  * Haug and Bharti [2022] T. Haug and K. Bharti, Quantum Sci. Technol. 7, 045019 (2022).
  * Meyer  _et al._ [2021] J. J. Meyer, J. Borregaard, and J. Eisert, npj Quantum Inf. 7, 1 (2021).
  * Meyer [2021] J. J. Meyer, Quantum 5, 539 (2021).
  * Beckey  _et al._ [2022] J. L. Beckey, M. Cerezo, A. Sone, and P. J. Coles, Phys. Rev. Res. 4, 013083 (2022).
  * Kaubruegger  _et al._ [2019] R. Kaubruegger, P. Silvi, C. Kokail, _et al._ , Phys. Rev. Lett. 123, 260505 (2019).
  * Koczor  _et al._ [2020] B. Koczor, S. Endo, T. Jones, _et al._ , New J. Phys. 22, 083038 (2020).
  * Ma  _et al._ [2021] Z. Ma, P. Gokhale, T.-X. Zheng, _et al._ , in  _2021 IEEE International Conference on Quantum Computing and Engineering (QCE)_ (IEEE, 2021).
  * Haug and Kim [2022] T. Haug and M. S. Kim, Phys. Rev. A 106, 052611 (2022).
  * Sagastizabal  _et al._ [2021] R. Sagastizabal, S. P. Premaratne, B. A. Klaver, _et al._ , npj Quantum Inf. 7, 130 (2021).
  * Lyu  _et al._ [2020] C. Lyu, V. Montenegro, and A. Bayat, Quantum 4, 324 (2020).
  * Nakanishi  _et al._ [2019] K. M. Nakanishi, K. Mitarai, and K. Fujii, Phys. Rev. Res. 1, 033062 (2019).
  * Higgott  _et al._ [2019] O. Higgott, D. Wang, and S. Brierley, Quantum 3, 156 (2019).
  * McClean  _et al._ [2017] J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. De Jong, Phys. Rev. A 95, 042308 (2017).
  * Santagati  _et al._ [2018] R. Santagati, J. Wang, A. A. Gentile, _et al._ , Sci. Adv. 4, eaap9646 (2018).
  * Hong  _et al._ [2023] C.-L. Hong, L. Colmenarez, L. Ding, _et al._ , Quantum parallelized variational quantum eigensolvers for excited states (2023), arXiv:2306.11844 [quant-ph] .
  * Feulner and Hartmann [2022] V. Feulner and M. J. Hartmann, Phys. Rev. B 106, 144426 (2022).
  * Barkoutsos  _et al._ [2018] P. K. Barkoutsos, J. F. Gonthier, I. Sokolov, _et al._ , Phys. Rev. A 98, 022322 (2018).
  * Wang  _et al._ [2009] H. Wang, S. Ashhab, and F. Nori, Phys. Rev. A 79, 042335 (2009).
  * Seki  _et al._ [2020] K. Seki, T. Shirakawa, and S. Yunoki, Phys. Rev. A 101, 052340 (2020).
  * Gard  _et al._ [2020] B. T. Gard, L. Zhu, G. S. Barron, _et al._ , npj Quantum Inf. 6, 10 (2020).
  * Barron  _et al._ [2021] G. S. Barron, B. T. Gard, O. J. Altman, _et al._ , Phys. Rev. Appl. 16, 034003 (2021).
  * Zhang  _et al._ [2021] F. Zhang, N. Gomes, N. F. Berthusen, _et al._ , Phys. Rev. Res. 3, 013039 (2021).
  * Zheng  _et al._ [2023] H. Zheng, Z. Li, J. Liu, _et al._ , PRX Quantum 4, 020327 (2023).
  * McClean  _et al._ [2016] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New J. Phys. 18, 023023 (2016).
  * Ryabinkin  _et al._ [2018] I. G. Ryabinkin, S. N. Genin, and A. F. Izmaylov, J. Chem. Theory Comput. 15, 249 (2018).
  * Lyu  _et al._ [2023] C. Lyu, X. Xu, M.-H. Yung, and A. Bayat, Quantum 7, 899 (2023).
  * Wiersema  _et al._ [2020] R. Wiersema, C. Zhou, Y. de Sereville, _et al._ , PRX Quantum 1, 020319 (2020).
  * Haegeman  _et al._ [2012] J. Haegeman, B. Pirvu, D. J. Weir, _et al._ , Phys. Rev. B 85, 100408 (2012).
  * Vicentini  _et al._ [2019] F. Vicentini, A. Biella, N. Regnault, and C. Ciuti, Phys. Rev. Lett. 122, 250503 (2019).
  * Huang  _et al._ [2022] Y. Huang, Q. Li, X. Hou, _et al._ , Phys. Rev. A 105, 052414 (2022).
  * Liu and Nocedal [1989] D. C. Liu and J. Nocedal, Math. Program. 45, 503 (1989).
  * Li and Benjamin [2017] Y. Li and S. C. Benjamin, Phys. Rev. X 7, 021050 (2017).
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, _et al._ , Nat. Commun. 9, 4812 (2018).
  * Holmes  _et al._ [2022] Z. Holmes, K. Sharma, M. Cerezo, and P. J. Coles, PRX Quantum 3, 010313 (2022).
  * Vatan and Williams [2004] F. Vatan and C. Williams, Phys. Rev. A 69, 032315 (2004).
  * De Chiara  _et al._ [2011] G. De Chiara, M. Lewenstein, and A. Sanpera, Phys. Rev. B 84, 054451 (2011).

##  Appendix A Qubit Hamiltonian of Spin-1 BBH model with Binary Encoding

Here, we provide the exact form of result obtained by applying the binary
encoding method to the spin-1 Bilinear-Biquadratic Heisenberg (BBH)
Hamiltonian, as shown in Eqs. (13)-(14). The qubit Hamiltonian for the spin-1
BBH model, after implementing the transformation described in Eq. (III), is as
follows:

$\displaystyle\vec{\widetilde{S}}_{i}\cdot\vec{\widetilde{S}}_{i+1}\mapsto\frac{1}{8}\big{(}\sigma_{j+1}^{x}\sigma_{j+3}^{x}+\sigma_{j+1}^{x}\sigma_{j+2}^{z}\sigma_{j+3}^{x}+\sigma_{j+1}^{x}\sigma_{j+2}^{x}$
(34)
$\displaystyle-\sigma_{j+1}^{x}\sigma_{j+2}^{x}\sigma_{j+3}^{z}+\sigma_{j}^{z}\sigma_{j+1}^{x}\sigma_{j+3}^{x}+\sigma_{j}^{z}\sigma_{j+1}^{x}\sigma_{j+2}^{z}\sigma_{j+3}^{x}$
$\displaystyle+\sigma_{j}^{z}\sigma_{j+1}^{x}\sigma_{j+2}^{x}-\sigma_{j}^{z}\sigma_{j+1}^{x}\sigma_{j+2}^{x}\sigma_{j+3}^{z}+\sigma_{j}^{x}\sigma_{j+3}^{x}$
$\displaystyle+\sigma_{j}^{x}\sigma_{j+2}^{z}\sigma_{j+3}^{x}+\sigma_{j}^{x}\sigma_{j+2}^{x}-\sigma_{j}^{x}\sigma_{j+2}^{x}\sigma_{j+3}^{z}$
$\displaystyle-\sigma_{j}^{x}\sigma_{j+1}^{z}\sigma_{j+3}^{x}-\sigma_{j}^{x}\sigma_{j+1}^{z}\sigma_{j+2}^{z}\sigma_{j+3}^{x}-\sigma_{j}^{x}\sigma_{j+1}^{z}\sigma_{j+2}^{x}$
$\displaystyle+\sigma_{j}^{x}\sigma_{j+1}^{z}\sigma_{j+2}^{x}\sigma_{j+3}^{z}+\sigma_{j+1}^{y}\sigma_{j+3}^{y}+\sigma_{j+1}^{y}\sigma_{j+2}^{z}\sigma_{j+3}^{y}$
$\displaystyle+\sigma_{j+1}^{y}\sigma_{j+2}^{y}-\sigma_{j+1}^{y}\sigma_{j+2}^{y}\sigma_{j+3}^{z}+\sigma_{j}^{z}\sigma_{j+1}^{y}\sigma_{j+3}^{y}$
$\displaystyle+\sigma_{j}^{z}\sigma_{j+1}^{y}\sigma_{j+2}^{z}\sigma_{j+3}^{y}+\sigma_{j}^{z}\sigma_{j+1}^{y}\sigma_{j+2}^{y}-\sigma_{j}^{z}\sigma_{j+1}^{y}\sigma_{j+2}^{y}\sigma_{j+3}^{z}$
$\displaystyle+\sigma_{j}^{y}\sigma_{j+3}^{y}+\sigma_{j}^{y}\sigma_{j+2}^{z}\sigma_{j+3}^{y}+\sigma_{j}^{y}\sigma_{j+2}^{y}-\sigma_{j}^{y}\sigma_{j+2}^{y}\sigma_{j+3}^{z}$
$\displaystyle-\sigma_{j}^{y}\sigma_{j+1}^{z}\sigma_{j+3}^{y}-\sigma_{j}^{y}\sigma_{j+1}^{z}\sigma_{j+2}^{z}\sigma_{j+3}^{y}-\sigma_{j}^{y}\sigma_{j+1}^{z}\sigma_{j+2}^{y}$
$\displaystyle+\sigma_{j}^{y}\sigma_{j+1}^{z}\sigma_{j+2}^{y}\sigma_{j+3}^{z}\big{)}+\frac{1}{4}\big{(}\sigma_{j}^{z}\sigma_{j+2}^{z}+\sigma_{j}^{z}\sigma_{j+3}^{z}$
$\displaystyle+\sigma_{j+1}^{z}\sigma_{j+2}^{z}+\sigma_{j+1}^{z}\sigma_{j+3}^{z}\big{)}$
$\displaystyle(\vec{\widetilde{S}}_{i}\cdot\vec{\widetilde{S}}_{i+1})^{2}\mapsto\frac{3}{4}I+\frac{1}{4}\Big{(}\sigma_{j+2}^{z}+\sigma_{j}^{z}-\sigma_{j+3}^{z}-\sigma_{j}^{z}\sigma_{j+3}^{z}$
(35)
$\displaystyle\quad-\sigma_{j+1}^{z}\sigma_{j+2}^{z}-\sigma_{j+1}^{z}\Big{)}+\frac{1}{8}\Big{(}\sigma_{j}^{x}\sigma_{j+1}^{x}\sigma_{j+2}^{x}\sigma_{j+3}^{x}+\sigma_{j}^{y}\sigma_{j+1}^{x}\sigma_{j+2}^{y}\sigma_{j+3}^{x}$
$\displaystyle\quad-\sigma_{j}^{x}\sigma_{j+1}^{x}\sigma_{j+2}^{y}\sigma_{j+3}^{y}+\sigma_{j}^{y}\sigma_{j+1}^{x}\sigma_{j+2}^{x}\sigma_{j+3}^{y}+\sigma_{j}^{z}\sigma_{j+1}^{x}\sigma_{j+2}^{x}\sigma_{j+3}^{z}$
$\displaystyle\quad+\sigma_{j}^{x}\sigma_{j+1}^{y}\sigma_{j+2}^{y}\sigma_{j+3}^{x}+\sigma_{j}^{y}\sigma_{j+1}^{y}\sigma_{j+2}^{y}\sigma_{j+3}^{y}+\sigma_{j}^{x}\sigma_{j+1}^{y}\sigma_{j+2}^{x}\sigma_{j+3}^{y}$
$\displaystyle\quad+\sigma_{j+1}^{x}\sigma_{j+2}^{x}\sigma_{j+3}^{z}-\sigma_{j}^{z}\sigma_{j+1}^{y}\sigma_{j+2}^{y}-\sigma_{j+1}^{x}\sigma_{j+2}^{x}+\sigma_{j}^{z}\sigma_{j+1}^{y}\sigma_{j+2}^{y}\sigma_{j+3}^{z}$
$\displaystyle\quad-\sigma_{j}^{y}\sigma_{j+3}^{y}+\sigma_{j}^{x}\sigma_{j+1}^{z}\sigma_{j+2}^{z}\sigma_{j+3}^{x}-\sigma_{j}^{y}\sigma_{j+2}^{z}\sigma_{j+3}^{y}+\sigma_{j}^{x}\sigma_{j+1}^{z}\sigma_{j+3}^{x}$
$\displaystyle\quad+\sigma_{j}^{y}\sigma_{j+1}^{z}\sigma_{j+3}^{y}-\sigma_{j}^{x}\sigma_{j+2}^{z}\sigma_{j+3}^{x}+\sigma_{j}^{y}\sigma_{j+1}^{z}\sigma_{j+2}^{z}\sigma_{j+3}^{y}$
$\displaystyle\quad-\sigma_{j}^{x}\sigma_{j+3}^{x}\Big{)}+\frac{1}{4}\Big{(}\sigma_{j+2}^{z}\sigma_{j+3}^{z}+\sigma_{j}^{z}\sigma_{j+1}^{z}+\sigma_{j}^{z}\sigma_{j+1}^{z}\sigma_{j+2}^{z}\sigma_{j+3}^{z}\Big{)}$