  The simulation of strongly correlated many-electron systems is one of the
most promising applications for near-term quantum devices. Here we use a class
of eigenvalue solvers (presented in Phys. Rev. Lett. 126, 070504 (2021)) in
which a contraction of the Schr\"odinger equation is solved for the
two-electron reduced density matrix (2-RDM) to resolve the energy splittings of
ortho-, meta-, and para-isomers of benzyne ${\textrm C_6} {\textrm H_4}$. In
contrast to the traditional variational quantum eigensolver, the contracted
quantum eigensolver solves an integration (or contraction) of the many-electron
Schr\"odinger equation onto the two-electron space. The quantum solution of the
anti-Hermitian part of the contracted Schr\"odinger equation (qACSE) provides a
scalable approach with variational parameters that has its foundations in 2-RDM
theory. Experimentally, a variety of error mitigation strategies enable the
calculation, including a linear shift in the 2-RDM targeting the iterative
nature of the algorithm as well as a projection of the 2-RDM onto the convex
set of approximately $N$-representable 2-RDMs defined by the 2-positive (DQG)
$N$-representability conditions. The relative energies exhibit single-digit
millihartree errors, capturing a large part of the electron correlation energy,
and the computed natural orbital occupations reflect the significant
differences in the electron correlation of the isomers.
