Introduction

Quantum computers are anticipated to enable simulations of quantum systems more efficiently and accurately than classical computers1,2. A promising algorithm to perform this task on emerging noisy intermediate-scale quantum (NISQ)3,4,5 computers is the variational quantum eigensolver (VQE)6,7,8,9,10,11,12. The VQE is a hybrid quantum-classical algorithm that estimates the lowest eigenvalue of a Hamiltonian H by minimizing the energy expectation value E(θ) = 〈ψ(θ)∣H∣ψ(θ)〉 with respect to a parametrized state \(\left|\psi (\boldsymbol\theta )\right\rangle =U(\boldsymbol\theta )\left|{\psi }_{0}\right\rangle\). Here, θ is a set of variational parameters, and the unitary U(θ) is an ansatz. Compared to other purely quantum algorithms for eigenvalue estimation, like the quantum-phase estimation algorithm13,14, the VQE requires shallower quantum circuits. This makes the VQE more noise resistant, at the expense of requiring a higher number of quantum measurements and additional classical postprocessing.

The VQE can solve the electronic structure problem6,15 by estimating the lowest eigenvalue of an electronic Hamiltonian. A major challenge for the practical realization of a molecular VQE simulation on NISQ computers is to construct a variationally flexible ansatz U(θ) that: (1) accurately approximates the ground state of H; (2) is easy to optimize; and (3) can be implemented by a shallow circuit.

These desired qualities are satisfied, to various levels, by several types of ansätze. The unitary coupled-cluster (UCC) type, was the first to be used for molecular VQE simulations16. The UCC is motivated by the classical coupled-cluster theory15, and corresponds to a series of unitary evolutions of fermionic-excitation operators, which we refer to as “fermionic excitation evolutions” (see the section “Ansatz elements”). A prominent example of a UCC ansatz is the UCC Singles and Doubles (UCCSD)17,18,19,20,21,22, which corresponds to a series of single and double-fermionic-excitation evolutions. The UCCSD has been used successfully to implement the VQE for small molecules16,17,23. Due to their physically motivated fermionic structure, UCC ansätze respects the symmetries of electronic wavefunctions, which makes these ansätze accurate and easy to optimize. Even a relatively simple UCC ansatz, like the UCCSD, is highly accurate for weakly correlated systems, such as molecules near their equilibrium configuration16,17,24,25. However, UCC ansätze are general-purpose built and do not take into account details of the system of interest. They contain redundant excitation terms, resulting in unnecessarily high numbers of variational parameters as well as long ansatz circuits. Moreover, to simulate strongly correlated systems, UCC ansätze requires higher-order excitations and/or multiple-step Trotterization24, which creates additional overhead for the quantum hardware.

Another type of “hardware-efficient” ansätze26,27,28,29,30 corresponds to universal unitary transformations implemented as periodic sequences of parametrized one- and two-qubit gates. These ansätze are implemented by shallow circuits, and can be highly variationally flexible. However, as they lack a physically motivated structure, these ansätze require a large number of variational parameters and may suffer by vanishing energy gradients along their variational parameters, making classical optimization intractable for large molecules31,32. In some scenarios, this is known as the barren plateau problem32,33,34,35.

Recently, a number of works36,37,38,39,40,41,42,43,44 suggested new “iterative” VQE protocols, which instead of using general-purpose, fixed ansätze, construct problem-tailored ansätze on the go. These algorithms can construct arbitrarily accurate ansätze that are optimized in the number of variational parameters and the ansatz-circuit depth, at the expense of requiring a larger number of quantum computer measurements. The ADAPT-VQE protocols36,37 are perhaps the most prominent family of iterative VQE protocols. The fermionic-ADAPT-VQE36, which was the first iterative VQE protocol, constructs its ansatz by iteratively appending parametrized unitary operators, which we refer to as “ansatz elements”. The ansatz element at each iteration is sampled from a pool of spin-complement single- and double-fermionic-excitation evolutions, based on an energy-gradient hierarchy. The fermionic-ADAPT-VQE was demonstrated to achieve chemical accuracy (10−3 Hartree), using an ansatz with several times fewer variational parameters, and a correspondingly shallower circuit, than the UCCSD. In the follow-up work37, the qubit-ADAPT-VQE utilizes an ansatz-element pool of more variationally flexible and rudimentary Pauli string exponentials. Due to this, the qubit-ADAPT-VQE constructs even shallower ansatz circuits than the fermionic-ADAPT-VQE, thus being, to the best of our knowledge, the currently most circuit-efficient, physically motivated VQE algorithm. However, the use of more rudimentary unitary operations comes at the expense of requiring additional variational parameters and iterations to construct an ansatz for a given accuracy.

In this work, we utilize unitary operations that, despite the lack of some of the physical features captured by fermionic-excitation evolutions, achieve the accuracy of fermionic excitations evolutions as well as the hardware efficiency of Pauli string exponentials. These operations can be used to construct circuit-efficient molecular ansätze without incurring as many additional variational parameters and iterations, as the qubit-ADAPT-VQE. We call these unitary operations “qubit excitation evolutions”. Qubit-excitation evolutions23,45,46,47 are unitary evolutions of “qubit excitation operators”, which satisfy “qubit commutation relations”46,47. Qubit-excitation evolutions can be implemented by circuits that act on fixed numbers of qubits, as opposed to fermionic-excitation evolutions, which act on a number of qubits that scale at least as \(O({{{{{{{\mathrm{log}}}}}}}\,}_{2}{N}_{{{{{{{{\rm{MO}}}}}}}}})\) with the number of considered molecular spin orbitals NMO. We show numerically that qubit-excitation evolutions can approximate an electronic wavefunction almost as accurately as fermionic-excitation evolutions can. On the other hand, qubit-excitation evolutions enjoy higher complexity than Pauli string exponentials, thus allowing for more rapid construction of the ansatz. We utilize qubit-excitation evolutions to introduce the qubit-excitation-based adaptive variational quantum eigensolver (QEB-ADAPT-VQE) protocol. As the name suggests, the QEB-ADAPT-VQE is an ADAPT-VQE protocol for molecular simulations that grows a problem-tailored ansatz from an ansatz-element pool of qubit-excitation evolutions. The QEB-ADAPT-VQE also features a modified ansatz-growing strategy, which allows for a more efficient ansatz construction at the expense of a constant-factor increase of quantum computer measurement. We benchmark the performance of the QEB-ADAPT-VQE with classical numerical simulations for small molecules: LiH, H6, and BeH2. In the section “Energy dissociation curves”, we compare the QEB-ADAPT-VQE to the standard UCCSD-VQE by presenting energy-dissociation curves obtained with each of the two methods. In the section “Energy convergence”, we compare the QEB-ADAPT-VQE to the fermionic-ADAPT-VQE and to the qubit-ADAPT-VQE by presenting energy convergence plots, obtained with each of the three ADAPT-VQE protocols.

Results

Theoretical background and notation

We begin with a theoretical introduction (required for the self-completeness of the paper) and by defining our notation. Finding the ground-state electron wavefunction \(\left|{E}_{0}\right\rangle\) and corresponding energy E0 of a molecule (or an atom) is known as the “electronic structure problem”15. This problem can be solved by solving the time-independent Schrödinger equation \(H\left|{{{\Phi }}}_{0}\right\rangle ={E}_{0}\left|{{{\Phi }}}_{0}\right\rangle\), where H is the electronic Hamiltonian of the molecule. Within the Born–Oppenheimer approximation, where the nuclei of the molecule are treated as motionless, H can be secondly quantized as

$$H=\mathop{\sum }\limits_{i,k}^{{N}_{{{{{{{{\rm{MO}}}}}}}}}}{h}_{i,k}{a}_{i}^{{{{\dagger}}} }{a}_{k}+\mathop{\sum }\limits_{i,j,k,l}^{{N}_{{{{{{{{\rm{MO}}}}}}}}}}{h}_{i,j,k,l}{a}_{i}^{{{{\dagger}}} }{a}_{j}^{{{{\dagger}}} }{a}_{k}{a}_{l}.$$
(1)

As already mentioned, NMO is the number of molecular spin orbitals, \({a}_{i}^{{{{\dagger}}} }\) and ai are fermionic creation and annihilation operators, corresponding to the ith molecular spin orbital, and the factors hij and hijkl are one- and two-electron integrals, written in a spin-orbital basis15. The Hamiltonian expression in Eq. (1) can be mapped to quantum-gate operators using an encoding method, e.g., the Jordan–Wigner48 or the Bravyi–Kitaev49 methods. Throughout this work, we assume the more straightforward Jordan–Wigner encoding, where the occupancy of the ith molecular spin orbital is represented by the state of the ith qubit.

The fermionic operators \({a}_{i}^{{{{\dagger}}} }\) and ai satisfy anticommutation relations

$$\{{a}_{i},{a}_{j}^{{{{\dagger}}} }\}={\delta }_{i,j},\,\,\{{a}_{i},{a}_{j}\}=\{{a}_{i}^{{{{\dagger}}} },{a}_{j}^{{{{\dagger}}} }\}=0.$$
(2)

Within the Jordan–Wigner encoding, \({a}_{i}^{{{{\dagger}}} }\) and ai can be written in terms of quantum-gate operators as

$${a}_{i}^{{{{\dagger}}} }={Q}_{i}^{{{{\dagger}}} }\mathop{\prod }\limits_{r=0}^{i-1}{Z}_{r}=\frac{1}{2}({X}_{i}-{{{{{{{\rm{i}}}}}}}}{Y}_{i})\mathop{\prod }\limits_{r=0}^{i-1}{Z}_{r}\ \ {{{{{{{\rm{and}}}}}}}}$$
(3)
$${a}_{i}={Q}_{i}\mathop{\prod }\limits_{r=0}^{i-1}{Z}_{r}=\frac{1}{2}({X}_{i}+{{{{{{{\rm{i}}}}}}}}{Y}_{i})\mathop{\prod }\limits_{r=0}^{i-1}{Z}_{r},$$
(4)

where

$${Q}_{i}^{{{{\dagger}}} }\equiv \frac{1}{2}({X}_{i}-{{{{{{{\rm{i}}}}}}}}{Y}_{i})\ \ {{{{{{{\rm{and}}}}}}}}$$
(5)
$$\,{Q}_{i}\equiv \frac{1}{2}({X}_{i}+{{{{{{{\rm{i}}}}}}}}{Y}_{i}).$$
(6)

We refer to \({Q}_{i}^{{{{\dagger}}} }\) and Qi as qubit creation and annihilation operators, respectively. They act to change the occupancy of spin-orbital i. The Pauli-z strings, in Eqs. (3) and (4), compute the parity of the state and act as exchange phase factors that account for the fermionic anticommutation of a† and a. Substituting Eqs. (3) and (4) into Eq. (1), H can be written as

$$H=\mathop{\sum}\limits_{r}{h}_{r}\mathop{\prod }\limits_{s=0}^{{N}_{MO}-1}{\sigma }_{s}^{r},$$
(7)

where σs is a Pauli operator (Xs, Ys, Zs, or Is) acting on qubit s, and hr (not to be confused with hik and hijkl) is a real scalar coefficient. The number of terms in EQ. (7) scales as \(O({N}_{{{{{{{{\rm{MO}}}}}}}}}^{4})\).

Once H is mapped to a Pauli operator representation, the VQE can be used to minimize the expectation value E(θ) = 〈ψ(θ)∣H∣ψ(θ)〉. The VQE relies upon the Rayleigh–Ritz variational principle

$$\langle \psi (\boldsymbol\theta )| H| \psi (\boldsymbol\theta )\rangle \ge {E}_{0},$$
(8)

to find an estimate for E0. The VQE is a hybrid quantum-classical algorithm that uses a quantum computer to prepare the trial state \(\left|\psi (\boldsymbol\theta )\right\rangle\) and evaluate E(θ), and a classical computer to process the measurement data and update θ at each iteration. The trial state \(\left|\psi (\boldsymbol\theta )\right\rangle =U(\boldsymbol\theta )\left|{\psi }_{0}\right\rangle\) is generated by an ansatz, U(θ), applied to an initial reference state \(\left|{\psi }_{0}\right\rangle\).

The ADAPT-VQE protocols

The ADAPT-VQE protocols iteratively construct problem-tailored ansätze on the go. At the mth iteration one or several unitary operators, \(\{{U}_{r}^{(m)}({\theta }_{r}^{(m)})\}\), which we refer to as ansatz elements, are appended to the left of the already existing ansatz, U(θ(m−1)):

$$U({\boldsymbol\theta }^{(m)})=\mathop{\prod}\limits_{r}{U}_{r}^{(m)}({\theta }_{r}^{(m)})U({\boldsymbol\theta }^{(m-1)})=\mathop{\prod }\limits_{p=m}^{1}\mathop{\prod}\limits_{r}{U}_{r}^{(p)}({\theta }_{r}^{(p)}).$$
(9)

The ansatz elements, \(\left\{{U}_{r}^{(m)}({\theta }_{r}^{(m)})\right\}\), at each iteration, are chosen from a finite ansatz-element pool \({\mathbb{P}}\), based on an ansatz-growing strategy that aims to achieve the lowest estimate of E(θ(m)). After a new ansatz U(θ(m)) is constructed, the new set of variational parameters \({\boldsymbol\theta }^{(m)}={\boldsymbol\theta }^{(m-1)}\cup \left\{{\theta }_{r}^{(m)}\right\}\) is optimized by the VQE, and a new estimate for E(θ(m)) is obtained. This iterative greedy strategy results in an ansatz that is tuned specifically to the system being simulated, and can approximate the ground eigenstate of the system with considerably fewer variational parameters and a shallower ansatz circuit, than general-purpose fixed ansätze, like the UCCSD.

In the fermionic-ADAPT-VQE, the ansatz-element pool \({\mathbb{P}}\) is a set of spin-complement pairs of single and double-fermionic-excitation evolutions. In the qubit-ADAPT-VQE, \({\mathbb{P}}\) is a set of parametrized exponentials of XY-Pauli strings. The growth strategy of the fermionic-ADAPT-VQE and the qubit-ADAPT-VQE is to add, at each iteration, the ansatz element with the largest energy-gradient magnitude

$${\left|\frac{\partial }{\partial {\theta }^{(m)}}\langle {\psi }^{(m-1)}| {U}^{(m){{{\dagger}}} }({\theta }^{(m)})H{U}^{(m)}({\theta }^{(m)})| {\psi }^{(m-1)}\rangle \right|}_{\theta = 0},$$

where \(\left|{\psi }^{(m-1)}\right\rangle\) is the trial state at the end of the (m−1)th iteration. For detailed descriptions of the fermionic-ADAPT-VQE and the qubit-ADAPT-VQE, we refer the reader to refs. 36 and 37, respectively.

Ansatz elements

Single- and double-fermionic-excitation evolutions can construct an ansatz that approximates an electronic wavefuction to an arbitrary accuracy50,51. Single and double-fermionic-excitation operators are defined, respectively, by the skew-Hermitian operators

$${T}_{ik}\equiv {a}_{i}^{{{{\dagger}}} }{a}_{k}-{a}_{k}^{{{{\dagger}}} }{a}_{i}\ \ {{{{{{{\rm{and}}}}}}}}$$
(10)
$${T}_{ijkl}\equiv {a}_{i}^{{{{\dagger}}} }{a}_{j}^{{{{\dagger}}} }{a}_{k}{a}_{l}-{a}_{k}^{{{{\dagger}}} }{a}_{l}^{{{{\dagger}}} }{a}_{i}{a}_{j}.$$
(11)

Single and double-fermionic-excitation evolutions are thus given, respectively, by the unitaries

$${A}_{ik}(\theta )={e}^{\theta {T}_{ik}}=\exp \left[\theta ({a}_{i}^{{{{\dagger}}} }{a}_{k}-{a}_{k}^{{{{\dagger}}} }{a}_{i})\right]\ \ {{{{{{{\rm{and}}}}}}}}$$
(12)
$${A}_{ijkl}(\theta )={e}^{\theta {T}_{ijkl}}=\exp \left[\theta ({a}_{i}^{{{{\dagger}}} }{a}_{j}^{{{{\dagger}}} }{a}_{k}{a}_{l}-{a}_{k}^{{{{\dagger}}} }{a}_{l}^{{{{\dagger}}} }{a}_{i}{a}_{j})\right].$$
(13)

Using Eqs. (3) and (4), for i < j < k < l, Aik and Aijkl can be expressed in terms of quantum-gate operators as

$${A}_{ik}(\theta )=\exp \left[{{{{{{{\rm{i}}}}}}}}\frac{\theta }{2}({X}_{i}{Y}_{k}-{Y}_{i}{X}_{k})\mathop{\prod }\limits_{r=i+1}^{k-1}{Z}_{r}\right]\ \ {{{{{{{\rm{and}}}}}}}}$$
(14)
$${A}_{ijkl}(\theta )= \, \exp \left[{{{{{{{\rm{i}}}}}}}}\frac{\theta }{8}\left({X}_{i}{Y}_{j}{X}_{k}{X}_{l}+{Y}_{i}{X}_{j}{X}_{k}{X}_{l}+{Y}_{i}{Y}_{j}{Y}_{k}{X}_{l}+{Y}_{i}{Y}_{j}{X}_{k}{Y}_{l}\right.\right.\\ \, \left.\left. \qquad-\,{X}_{i}{X}_{j}{Y}_{k}{X}_{l}-{X}_{i}{X}_{j}{X}_{k}{Y}_{l}-{Y}_{i}{X}_{j}{Y}_{k}{Y}_{l}-{X}_{i}{Y}_{j}{Y}_{k}{Y}_{l}\right)\mathop{\prod }\limits_{r=i+1}^{j-1}{Z}_{r}\mathop{\prod }\limits_{{r}^{\prime}=k+1}^{l-1}{Z}_{{r}^{\prime}}\right].$$
(15)

As seen from Eqs. (14) and (15), fermionic-excitation evolutions act on a number of qubits that scales as O(NMO). Therefore, they are implemented by circuits whose size (in terms of number of CNOTs) also scales as O(NMO). We derived a CNOT-efficient method to construct circuits for fermionic excitations evolutions in ref.  47. The circuits for a single and double-fermionic-excitation evolution have CNOT counts of 2(k − i) + 1 and 2(l + j − i − k) + 9, respectively.

Qubit-excitation operators are defined by the qubit annihilation and creation operators, Qi and \({Q}_{i}^{{{{\dagger}}} }\) (Eqs. (5) and (6)), which satisfy the qubit-commutation relations

$$\{{Q}_{i},{Q}_{i}^{{{{\dagger}}} }\}={\delta }_{i,j},[{Q}_{i},{Q}_{j}^{{{{\dagger}}} }]=0\,{{\mbox{if}}}\,\,i\,\ne\, j,\ \ {{{{{{{\rm{and}}}}}}}}\ \ [{Q}_{i},{Q}_{j}]=[{Q}_{i}^{{{{\dagger}}} },{Q}_{j}^{{{{\dagger}}} }]=0\,{{{{{\mathrm{for}}}}}}\; {{{{{\mathrm{all}}}}}}\,i,j.$$
(16)

Some authors have referred to these commutation relations as parafermionic46. Single- and double-qubit-excitation operators are given, respectively, by the skew-Hermitian operators

$${\tilde{T}}_{ik}\equiv {Q}_{i}^{{{{\dagger}}} }{Q}_{k}-{Q}_{k}^{{{{\dagger}}} }{Q}_{i}\ \ {{{{{{{\rm{and}}}}}}}}$$
(17)
$${\tilde{T}}_{ijkl}\equiv {Q}_{i}^{{{{\dagger}}} }{Q}_{j}^{{{{\dagger}}} }{Q}_{k}{Q}_{l}-{Q}_{k}^{{{{\dagger}}} }{Q}_{l}^{{{{\dagger}}} }{Q}_{i}{Q}_{j}.$$
(18)

Thus, single and double-qubit-excitation evolutions are given, respectively, by the unitary operators

$${\tilde{A}}_{ik}(\theta )={e}^{\theta {\tilde{T}}_{ik}}=\exp \left[\theta ({Q}_{i}^{{{{\dagger}}} }{Q}_{k}-{Q}_{k}^{{{{\dagger}}} }{Q}_{i})\right]\ \ {{{{{{{\rm{and}}}}}}}}$$
(19)
$${\tilde{A}}_{ijkl}(\theta )={e}^{\theta {\tilde{T}}_{ijkl}}=\exp \left[\theta ({Q}_{i}^{{{{\dagger}}} }{Q}_{j}^{{{{\dagger}}} }{Q}_{k}{Q}_{l}-{Q}_{k}^{{{{\dagger}}} }{Q}_{l}^{{{{\dagger}}} }{Q}_{i}{Q}_{j})\right].$$
(20)

Using Eqs. (5) and (6), \({\tilde{A}}_{ik}\) and \({\tilde{A}}_{ijkl}\) can be re-expressed in terms of quantum-gate operators as

$${\tilde{A}}_{ik}(\theta )=\exp \left[{{{{{{{\rm{i}}}}}}}}\frac{\theta }{2}({X}_{i}{Y}_{k}-{Y}_{i}{X}_{k})\right]\ \ {{{{{{{\rm{and}}}}}}}}$$
(21)
$${\tilde{A}}_{ijkl}(\theta )= \, \exp \left[{{{{{{{\rm{i}}}}}}}}\frac{\theta }{8}\left({X}_{i}{Y}_{j}{X}_{k}{X}_{l}+{Y}_{i}{X}_{j}{X}_{k}{X}_{l}+{Y}_{i}{Y}_{j}{Y}_{k}{X}_{l}+{Y}_{i}{Y}_{j}{X}_{k}{Y}_{l}\right.\right.\\ \, \left.\left.-{X}_{i}{X}_{j}{Y}_{k}{X}_{l}-{X}_{i}{X}_{j}{X}_{k}{Y}_{l}-{Y}_{i}{X}_{j}{Y}_{k}{Y}_{l}-{X}_{i}{Y}_{j}{Y}_{k}{Y}_{l}\right)\right].$$
(22)

As seen from Eqs. (21) and (22), unlike fermionic-excitation evolutions, qubit-excitation evolutions act on a fixed number of qubits, and can be implemented by circuits that have a fixed number of CNOTs. Single-qubit-excitation evolutions can be performed by the circuit in Fig. 1, with a CNOT count of 2. Double-qubit-excitation evolutions can be performed by the circuit in Fig. 2, which was introduced in ref.  47, with a CNOT count of 13.

Fig. 1: A circuit to implement a single-qubit-excitation evolution.
figure 1

A single-qubit-excitation evolution is defined by the unitary operator \({\tilde{A}}_{ik}(\theta )=\exp \left[{{{{{{{\rm{i}}}}}}}}\frac{\theta }{2}({X}_{i}{Y}_{k}-{Y}_{i}{X}_{k})\right]\), where X and Y are the Pauli x and y operators (the subscript denotes the qubit on which these operators act). qi denote the state of qubit i. Rx(θ) and Rz(θ) denote single-qubit rotation gates around the x and z axes, respectively, by θ.

Full size image
Fig. 2: A circuit to implement a double-qubit-excitation evolution.
figure 2

A double-qubit-excitation evolution is defined by the unitary operator \({\tilde{A}}_{ijkl}(\theta )=\exp \left[{{{{{{{\rm{i}}}}}}}}\frac{\theta }{8}({X}_{i}{Y}_{j}{X}_{k}{X}_{l}+{Y}_{i}{X}_{j}{X}_{k}{X}_{l}+{Y}_{i}{Y}_{j}{Y}_{k}{X}_{l}+{Y}_{i}{Y}_{j}{X}_{k}{Y}_{l}-{X}_{i}{X}_{j}{Y}_{k}{X}_{l}-{X}_{i}{X}_{j}{X}_{k}{Y}_{l}-{Y}_{i}{X}_{j}{Y}_{k}{Y}_{l}-{X}_{i}{Y}_{j}{Y}_{k}{Y}_{l})\right]\), where X and Y are the Pauli x and y operators. qi denote the state of qubit i. H denotes the Hadamard gate (not to be confused with the molecular Hamiltonian), and Ry(θ) and Rz(θ) denote single-qubit rotation gates around the y and z axes, respectively, by θ.

Full size image

For larger systems, qubit-excitation evolutions are increasingly more CNOT-efficient compared to fermionic-excitation evolutions, whose CNOT count scales as O(NMO) in the Jordan–Wigner encoding and as \(O({{{{{{\mathrm{log}}}}}}}\,{N}_{{{{{{{{\rm{MO}}}}}}}}})\) in the Bravyi–Kitaev encoding. On the other hand, single- and double-qubit-excitation evolutions, as seen from Eqs. (21) and (22), correspond to combinations of 2 and 8, mutually commuting Pauli string exponentials, respectively. Hence, by constructing ansätze with qubit-excitation evolutions instead of Pauli string exponentials, we decrease the number of variational parameters. A further advantage of qubit-excitation evolutions is that they allow for the local circuit optimizations of ref.  47, which Pauli string exponentials do not.

When comparing the QEB-ADAPT-VQE with the fermionic-ADAPT-VQE (see the section “Energy convergence”), we assume the use of the qubit- and fermionic-excitation evolutions circuits derived in ref.  47. To our knowledge, these are the most CNOT-efficient circuits for these two types of unitary operations. For the qubit-ADAPT-VQE, we assume that an exponential of a Pauli string of length l is implemented by a standard CNOT staircase construction6,47,52, with a CNOT count of 2(l − 1). Global circuit optimization is beyond the scope of this paper.

The QEB-ADAPT-VQE protocol

In the previous section, we formally introduced qubit-excitation evolutions and presented the circuits that implement such unitary evolutions. Here, we describe the three preparation components, and the fourth iterative component, of the QEB-ADAPT-VQE protocol.

First, we transform the molecular Hamiltonian H to a quantum-gate-operator representation as described earlier. This transformation is a standard step in every VQE algorithm. It involves the calculation of the one- and two-electron integrals hik and hijkl (Eq. (1)), which can be done efficiently (in time polynomial in NMO) on a classical computer6.

Second, we define an ansatz-element pool \({\mathbb{P}}(\tilde{A},{N}_{{{{{{{{\rm{MO}}}}}}}}})\) of all unique single and double-qubit-excitation evolutions, \({\tilde{A}}_{ik}(\theta )\) and \({\tilde{A}}_{ijkl}(\theta )\), respectively, for i, j, k, l ∈ {0, NMO − 1}. The size of this pool is \(| | {\mathbb{P}}(\tilde{A},{N}_{{{{{{{{\rm{MO}}}}}}}}})| | =\left({{N}_{{{{{{{{\rm{MO}}}}}}}}}} \atop {2}\right)+3\left({{N}_{{{{{{{{\rm{MO}}}}}}}}}} \atop {4}\right)\). Here, ∣∣ ⋅ ∣∣ denotes a set’s cardinality.

Third, we choose an initial reference state \(\left|{\psi }_{0}\right\rangle\). For faster convergence, \(\left|{\psi }_{0}\right\rangle\) should have a significant overlap with the unknown ground state, \(\left|{E}_{0}\right\rangle\). In the classical numerical simulations presented in this paper, we use the conventional choice of the Hartree–Fock state53.

Fourth, we initialize the iteration number to m = 1, and the ansatz to the identity U → U(0) = I. Then, we initiate the QEB-ADAPT-VQE iterative loop. We start by describing the six steps of the mth iteration of the QEB-ADAPT-VQE. We then comment on these steps.

  1. 1.

    Prepare state \(\left|{\psi }^{(m-1)}\right\rangle =U({\boldsymbol\theta }^{(m-1)})\left|{\psi }_{0}\right\rangle\), with θ(m−1) as determined in the previous iteration.

  2. 2.

    For each qubit-excitation evolution \({\tilde{A}}_{p}({\theta }_{p})={e}^{{\theta }_{p}{\tilde{T}}_{p}}\in {\mathbb{P}}(\tilde{A},{N}_{{{{{{{{\rm{MO}}}}}}}}})\), calculate the energy gradient:

    $$\, {\left.\frac{\partial }{\partial {\theta }_{p}}{E}^{(m-1)}({\theta }_{p})\right|}_{{\theta }_{p} = 0}={\left.\frac{\partial }{\partial {\theta }_{p}}\left\langle {\psi }^{(m-1)}\right|{\tilde{A}}_{p}^{{{{\dagger}}} }({\theta }_{p})H{\tilde{A}}_{p}({\theta }_{p})\left|{\psi }^{(m-1)}\right\rangle \right|}_{{\theta }_{p} = 0}\\ \, \qquad\qquad\quad\qquad\;\; =\left\langle {\psi }^{(m-1)}\right|[H,\tilde{{T}_{p}}]\left|{\psi }^{(m-1)}\right\rangle .$$
    (23)
  3. 3.

    Identify the set of n qubit-excitation evolutions, \({\tilde{{\mathbb{A}}}}^{(m)}(n)\), with largest energy gradient magnitudes. For \({\tilde{A}}_{p}({\theta }_{p})\in {\tilde{{\mathbb{A}}}}^{(m)}(n)\):

    • (a) Run the VQE to find \(\mathop{\min }\nolimits_{{\boldsymbol\theta }^{(m-1)},{\theta }_{p}}E({\boldsymbol\theta }^{(m-1)},{\theta }_{p})=\mathop{\min }\nolimits_{{\boldsymbol\theta }^{(m-1)},{\theta }_{p}}\left\langle {\psi }_{0}\right|{U}^{{{{\dagger}}} }({\boldsymbol\theta }^{(m-1)}){\tilde{A}}_{p}^{{{{\dagger}}} }({\theta }_{p})H{\tilde{A}}_{p}({\theta }_{p})U({\boldsymbol\theta }^{(m-1)})\left|{\psi }_{0}\right\rangle .\)

    • (b) Calculate the energy reduction \({{\Delta }}{E}_{p}^{(m)}={E}^{(m-1)}-\mathop{\min }\nolimits_{{\boldsymbol\theta }^{(m-1)},{\theta }_{p}}E({\boldsymbol\theta }^{(m-1)},{\theta }_{p})\) for each p.

    • (c) Save the (re)optimized values of θ(m−1) ∪ {θp} as \({\boldsymbol\theta }_{p}^{(m)}\) for each p.

  4. 4.

    Identify the largest energy reduction \({{\Delta }}{E}^{(m)}\equiv {{\Delta }}{E}_{p^{\prime} }^{(m)}=\mathop{\max }\nolimits_{p}\{{{\Delta }}{E}_{p}^{(m)}\}\), and the corresponding qubit-excitation evolution \({\tilde{A}}^{(m)}({\theta }^{(m)})\equiv {\tilde{A}}_{p^{\prime} }({\theta }_{p^{\prime} })\).

    If ΔE(m) < ϵ, where ϵ > 0 is an energy threshold:

    1. (a)

      Exit

    Else:

    1. (a)

      Append \({\tilde{A}}^{(m)}({\theta }^{(m)})\) to the ansatz: \(U({\boldsymbol\theta }^{(m)})={\tilde{A}}^{(m)}({\theta }^{(m)})U({\boldsymbol\theta }^{(m-1)})\)

    2. (b)

      Set \({E}^{(m)}={E}^{(m-1)}-{{\Delta }}{E}_{p^{\prime} }^{(m)}\)

    3. (c)

      Set the values of the new set of variational parameters, \({\boldsymbol\theta }^{(m)}={\boldsymbol\theta }^{(m-1)}\cup \{{\theta }_{p^{\prime} }\}\), to \({\boldsymbol\theta }_{p^{\prime} }^{(m)}\)

  5. 5.

    (Optional) If the ground state of the system of interest is known, a priori, to have the same spin as \(\left|{\psi }_{0}\right\rangle\), append to the ansatz the spin-complementary of \({\tilde{A}}^{(m)}({\theta }^{(m)})\), \({\tilde{A}}^{^{\prime} (m)}({\theta }^{^{\prime} (m)})\), unless \({\tilde{A}}^{(m)}({\theta }^{(m)})\equiv {\tilde{A}}^{^{\prime} (m)}\big({\theta }^{^{\prime} (m)}\big)\):

    $$U({\boldsymbol\theta }^{(m)})={\tilde{A}}^{^{\prime} (m)}({\theta }^{^{\prime} (m)}){\tilde{A}}^{(m)}({\theta }^{(m)})U({\boldsymbol\theta }^{(m-1)}).$$
    (24)
  6. 6.

    Enter the m + 1 iteration by returning to step 1

We now provide some more information about the steps of the protocol. The QEB-ADAPT-VQE loop starts by preparing the trial state \(\left|{\psi }^{(m-1)}\right\rangle\) obtained in the (m − 1)th iteration. To identify a suitable qubit-excitation evolution to append to the ansatz, first we calculate (step 2) the gradient of the energy expectation value, with respect to the variational parameter of each qubit-excitation evolution in \({\mathbb{P}}(\tilde{A},{N}_{{{{{{{{\rm{MO}}}}}}}}})\).The gradients are evaluated at θp = 0 because of the presumption that \(\left|{\psi }_{0}\right\rangle\) is close to the ground state, which suggests that the optimized value of θp is close to 0. The gradients (Eq. (23)) are calculated by measuring, on a quantum computer, the expectation value of the commutator of H and the corresponding qubit-excitation operator \({\tilde{T}}_{p}\), with respect to \(\left|{\psi }^{(m-1)}\right\rangle\). The expression for the gradient in Eq. (23) is derived explicitly in Supplementary Note 1. Note that, steps 1 and 2 are identical to those of the original fermionic-ADAPT-VQE.

The gradients calculated in step 2, indicate how much each qubit excitation can decrease E(m − 1). However, the largest gradient does not necessarily correspond to the largest energy reduction, optimized over all variational parameters. In step 3, we identify the set of n qubit-excitation evolutions with the largest energy-gradient magnitudes: \({\tilde{{\mathbb{A}}}}^{(m)}(n)\in {\mathbb{P}}(\tilde{A},{N}_{{{{{{{{\rm{MO}}}}}}}}})\). We assume that \({\tilde{{\mathbb{A}}}}^{(m)}(n)\) likely contains the qubit-excitation evolution that reduces E(m − 1) the most. For each of the n qubit-excitation evolutions in \({\tilde{{\mathbb{A}}}}^{(m)}(n)\), we run the VQE with the ansatz from the previous iteration to calculate how much it contributes to the energy reduction. Step 3 is not present in the original fermionic-ADAPT-VQE, which directly grows its ansatz by the ansatz element with largest energy-gradient magnitude (equivalent to n = 1). Performing step 3 for n > 1 further reduces the ansatz circuit at the expense of more quantum computer measurements. A study of the performance of the QEB-ADAPT-VQE for different values of n is presented in Supplementary Note 5. The study shows that for the three molecules considered in this paper, LiH, H6, and BeH2, a CNOT reduction between 15 and 25% is achieved for n = 10.

In step 4, we pick the qubit excitation, \({\tilde{A}}^{(m)}({\theta }^{(m)})\), with the largest contribution to the energy reduction, ΔE(m). If ΔE(m) is below some threshold ϵ > 0, we exit the iterative loop. If instead the ∣ΔE(m)∣ > ϵ, we add \({\tilde{A}}^{(m)}({\theta }^{(m)})\) to the ansatz.

If it is known, a priori, that the ground state of the simulated system has spin zero as the Hartree–Fock state does, we assume that qubit-excitation evolutions come in spin-complement pairs. Hence, we append the spin-complement of \({\tilde{A}}^{(m)}({\theta }^{(m)})\), \(\tilde{A}{^{\prime} }^{(m)}(\theta {^{\prime} }^{(m)})\) (step 5) to the ansatz. However, unlike the fermionic-ADAPT-VQE, the QEB-ADAPT-VQE assigns independent variational parameters to the two spin-complement excitation evolutions. The reason for this is that qubit-excitation evolutions do not account for the parity of the state. Hence, additional variational flexibility is required to obtain the correct relative sign between the two spin-complement qubit-excitation evolutions. Performing step 5 roughly halves the number of iterations required to construct an ansatz for a particular accuracy.

In Supplementary Note 4, we discuss the computational complexity of the QEB-ADAPT-VQE. As a worst-case estimate, the QEB-ADAPT-VQE might require as many as \(O(n{{N}_{MO}}^{16})\) quantum computer measurements.

Classical numerical simulations

We perform classical numerical VQE simulations for LiH, H6, and BeH2 to compare the use of qubit and fermionic excitations in the construction of molecular ansätze and to benchmark the performance of the QEB-ADAPT-VQE. LiH and BeH2 have been simulated with VQE-based protocols on real quantum computers and are often used in the field of quantum-computational chemistry to classically benchmark various VQE protocols17,18,36,39,40. Similar to refs.  36,37, we use H6 as a prototype of a molecule with a strongly correlated ground state. Our numerical results are based on a custom code, designed to implement ADAPT-VQE protocols for arbitrary ansatz-element pools and ansatz-growing strategies. The code is optimized to analytically calculate excitation-based state vectors (see Supplementary Note 2). The code uses the openfermion-psi454 package to second-quantize the Hamiltonian, and subsequently to transform it to quantum-gate-operator representation. For all simulations presented in this paper, we represent the molecular Hamiltonians in the Slater type orbital-3 Gaussians (STO-3G) spin-orbital basis set55,56, without assuming frozen orbitals. In this basis set, LiH, H6, and BeH2, have 12, 12, and 14 spin orbitals, respectively, which are represented by 12, 12, and 14 qubits. For the optimization of variational parameters, we use the gradient-descend Broyden Fletcher Goldfarb Shannon (BFGS) minimization method57 from Scipy58. We also supply to the BFGS an analytically calculated energy-gradient vector (see Supplementary Note 3), for faster optimization. We note that in the presence of high noise levels, gradient-descend minimizers are likely to struggle to find the global energy minimum59,60, while direct search minimizers, like the Nelder–Mead61, are likely to perform better62,63.

Qubit versus fermionic excitations

In this section, we compare qubit and fermionic-excitation evolutions in their ability to construct ansätze to approximate electronic wavefunctions. Directly comparing the QEB-ADAPT-VQE and the fermionic-ADAPT-VQE (as we do in the section “Energy convergence”) does not constitute a fair comparison of the two types of excitation evolutions: the QEB-ADAPT-VQE assigns one variational parameter per qubit-excitation evolution in its ansatz, whereas the fermionic-ADAPT-VQE assigns one variational parameter per spin-complement pair of fermionic-excitation evolutions. Consequently, here we compare the QEB-ADAPT-VQE for n = 1 and step 5 not implemented, to the fermionic-ADAPT-VQE when it grows its ansatz by appending individual fermionic-excitation evolutions (instead of spin-complement pairs of fermionic-excitation evolutions). In this way, the two protocols differ only in using a pool of qubit-excitation evolutions, and a pool of fermionic-excitation evolutions, respectively.

Figure 3 shows energy-convergence plots, obtained with the two protocols as explained above, for the ground states of LiH (Fig. 3a), H6 (Fig. 3b), and BeH2 (Fig. 3c) at bond distances of rLi-H = 1.546 Å, rH-H = 1.5 Å, and rBe-H = 1.316 Å, respectively. All plots are terminated for ϵ = 10−12 Hartree. The two protocols converge similarly, with the fermionic-ADAPT-VQE converging slightly faster for more than ~50 ansatz elements. This difference is most evident for the more strongly correlated H6 (Fig. 3b), where the fermionic-ADAPT-VQE requires up to 20% fewer excitation evolutions than the QEB-ADAPT-VQE to achieve a given accuracy. These observations suggest that fermionic-excitation-based ansätze might be able to approximate strongly correlated states a bit better than qubit-excitation-based ansätze. To further investigate this observation, in Fig. 4 we include energy-convergence plots, similar to those in Fig. 3, but for bond distances of rLi-H = 3 Å (Fig. 4a), rH-H = 3 Å (Fig. 4b), and rBe-H = 3 Å (Fig. 4c). At larger bond distances the ground states of the LiH, and BeH2 are more strongly correlated, so we expect to see a larger difference in the convergence rates of the two protocols.

Fig. 3: Comparison of the qubit and fermionic-excitation evolutions at equilibrium bond distances.
figure 3

The three subfigures present energy-convergence plots for the ground states of: a LiH, b H6, and c BeH2, in the STO-3G orbital basis set, at bond distances of rLi−H = 1.546 Å, rH−H = 1.5 Å, and rBe−H = 1.316 Å, respectively. The blue plots are obtained with the QEB-ADAPT-VQE for n = 1 and step 5 not implemented. The red plots are obtained with the fermionic-ADAPT-VQE using an ansatz-element pool of non-spin-complement fermionic-excitation evolutions. The plots are terminated at ϵ = 10−12 Hartree.

Full size image
Fig. 4: Comparison of the qubit and fermionic-excitation evolutions at large bond distances.
figure 4

The three subfigures present energy-convergence plots for the ground states of: a LiH, b H6, and c BeH2, in the STO-3G orbital basis set, at bond distances of rLi−H = 3 Å, rH−H = 3 Å, and rBe−H = 3 Å, respectively. The blue plots are obtained with the QEB-ADAPT-VQE for n = 1 and step 5 not implemented. The red plots are obtained with the fermionic-ADAPT-VQE using an ansatz-element pool of non-spin-complement fermionic-excitation evolutions. The plots are terminated at ϵ = 10−12 Hartree.

Full size image

In Fig. 4a, c, we see that for LiH and BeH2, at rLi–H = 3 Å and rBe–H = 3 Å, respectively, indeed there is a larger difference in the convergence rates of the two protocols, in favor of the fermionic-ADAPT-VQE. This is more evident for BeH2 where the fermionic-ADAPT-VQE requires ~20% fewer ansatz elements, on average, than the QEB-ADAPT-VQE, to achieve a given accuracy. These results further indicate that fermionic-excitation-based ansätze can approximate strongly correlated states better than qubit-excitation-based ansätze.

Energy-dissociation curves

Figure 5 shows energy-dissociation curves for LiH, H6, and BeH2, obtained with the QEB-ADAPT-VQE for n = 10 and energy-reduction thresholds ϵ4 = 10−4 Hartree, ϵ6 = 10−6 Hartree and ϵ8 = 10−8 Hartree. Dissociation curves obtained with the Hartree–Fock (HF) method, the full configuration interaction (FCI) method, and the VQE, using an untrotterized UCCSD ansatz (UCCSD-VQE) are also included for comparison. The UCCSD includes spin-conserving single and double-fermionic evolutions only, for a fairer comparison to the QEB-ADAPT-VQE.

Fig. 5: Energy-dissociation curves for LiH, H6, and BeH2 molecules in the STO-3G orbital basis set.
figure 5

a–c Absolute energy as a function of bond distance. d–f Energy error with respect to the exact FCI energy as a function of bond distance. g–i Number of ansatz variational parameters required to reach the energy accuracies in (d–f). The QEB-ADAPT-VQE is performed for n = 10 and step 5 implemented. The number of variational parameters for the UCCSD is 92, 117, and 204 for LiH, H6, and BeH2, respectively. The number of variational parameters is also equivalent to the number of ansatz elements of each ansatz.

Full size image

Figure 5a–c shows the absolute values for the ground-state energy estimates. All methods except the HF, produce close energy estimates that cannot be clearly distinguished. In Fig. 5d–f, the exact FCI energy is subtracted in order to differentiate better the different methods and their corresponding errors.

The UCCSD-VQE achieves chemical accuracy over all bond distances for LiH (Fig. 5d) and over bond distances close to equilibrium configuration for H6 (Fig. 5e) and BeH2 (Fig. 5f). However, the UCCSD-VQE fails to achieve chemical accuracy for bond distances away from equilibrium configuration for H6 and BeH2, where the ground states become more strongly correlated.

The QEB-ADAPT-VQE for ϵ4, similarly to the UCCSD-VQE, struggles to achieve chemical accuracy for strongly correlated ground states. However, for ϵ6 and ϵ8 the QEB-ADAPT-VQE achieves chemical accuracy over all investigated bond distances, for all three molecules. This indicates that the QEB-ADAPT-VQE can successfully construct ansätze to accurately approximate strongly correlated states.

However, the real strength of the QEB-ADAPT-VQE, similarly to other ADAPT-VQE protocols, is not just in constructing accurate ansätze, but in constructing accurate problem-tailored ansätze with few variational parameters, and corresponding shallow ansatz circuits. Figure 5g–i shows plots of the number of variational parameters used by the ansatz of each method as a function of bond distance. In the cases of LiH (Fig. 5g) and BeH2 (Fig. 5i), the ansätze constructed by the QEB-ADAPT-VQE for ϵ6 and ϵ8 are not only more accurate than the UCCSD but also have significantly fewer parameters. However, in the case of H6, the QEB-ADAPT-VQE on average requires more parameters than the UCCSD. The reason for this is that H6 is more strongly correlated than LiH and BeH2, so even an optimally constructed ansatz would require more variational parameters than the UCCSD, to accurately approximate the ground state of H6.

An interesting observation is the abrupt changes in the number of variational parameters used by the QEB-ADAPT-VQE for H6 at bond distances of around 1, 2, and 2.75 Å. The reason for these changes are molecular structure transformations, where different eigenstates of H become lowest in energy (energy-level crossings).

Energy convergence

In this section, we compare the QEB-ADAPT-VQE against the fermionic-ADAPT-VQE and the qubit-ADAPT-VQE using energy-convergence plots (see Fig. 6). To ensure a fair comparison we choose the following settings for the three protocols: We perform the QEB-ADAPT-VQE for n = 1, using an ansatz-element pool of all unique single- and double-qubit-excitation evolutions. The fermionic-ADAPT-VQE is performed as in ref. 36, using an ansatz-element pool of all unique single and double spin-complement fermionic-excitation evolutions. For the qubit-ADAPT-VQE, we use an ansatz element of all evolutions of XY-Pauli strings of length 2 and 4 that have an odd number of Ys. This pool consists of O(N4) Pauli string evolutions that can be combined to obtain all qubit-excitation evolutions in the ansatz element of the QEB-ADAPT-VQE (see the section “Ansatz elements”). Because of this, the comparison between the QEB-ADAPT-VQE and qubit-ADAPT-VQE, in terms of ansatz-circuit efficiency, can be considered fair. We note that the authors of ref. 37 proved that the qubit-ADAPT-VQE actually can construct an ansatz that exactly recovers the FCI wavefunction, using a reduced ansatz-element pool of only 2NMO − 2 Pauli string evolutions. This reduced pool can decrease the number of quantum computer measurements required to evaluate the energy gradients at each iteration (see step 2 of the QEB-ADAPT-VQE) from \(O({N}_{MO}^{8})\) to \(O({N}_{MO}^{5})\). However, the reduced ansatz-element pool will also result in a slower and less circuit-efficient ansatz construction, so using this reduced pool in the comparison with the QEB-ADAPT-VQE would not be fair.

Fig. 6: Comparison of the QEB-ADAPT-VQE, the fermionic-ADAPT-VQE, and the qubit-ADAPT-VQE.
figure 6

The subfigures above present energy-convergence plots for the ground states of LiH, H6, and BeH2, in the STO-3G orbital basis set, at bond distances rLi-H = 1.546 Å, rH–H = 1.5 Å, and rBe–H = 1.316 Å. The plots compare the QEB-ADAPT-VQE (blue), the fermionic-ADAPT-VQE (red), and the qubit-ADAPT-VQE (green) protocols in terms of the number of iterations (a–c), the number of parameters (d–f), and the number of CNOTs (g–i). The QEB-ADAPT-VQE is performed for n = 1. All convergence plots are terminated for an energy-reduction threshold of ϵ = 10−12 Hartree. The CNOT counts in (g–i) are obtained assuming the use of the quantum circuits discussed in the section “Ansatz elements''.

Full size image

We compare the three protocols in terms of three cost metrics, required to construct an ansatz to achieve a specific accuracy: (1) the number of iterations; (2) the number of variational parameters; and (3) the number of CNOTs. The number of iterations and the number of variational parameters (the number of iterations is the same as the number of variational parameters for the fermionic-ADAPT-VQE and the qubit-ADAPT-VQE, but not for the QEB-ADAPT-VQE) determine the total number of quantum computer measurements (see Supplementary Note 4). The CNOT count of the ansatz circuit is approximately proportional to its depth. Hence, the CNOT count can be used as a measure of the run time of the quantum subroutine of the VQE, which also reflects the error accumulated by the quantum hardware. Due to the limited coherence times of NISQ computers, the CNOT count is considered as a primary cost metric.

Figure 6 shows energy-convergence plots, obtained with the three ADAPT-VQE protocols, for LiH, H6, and BeH2 at bond distances of rLi−H = 1.546 Å, rH−H = 1.5 Å, and rBe−H = 1.316 Å, respectively. All energy-convergence plots are terminated at ϵ = 10−12 Hartree.

In Fig. 6a–c, we notice that the QEB-ADAPT-VQE and the fermionic-ADAPT-VQE perform similarly in terms of the number of iterations. This implies that the QEB-ADAPT-VQE and the fermionic-ADAPT-VQE use approximately the same number of the qubit and fermionic-excitation evolutions, respectively, when constructing their respective ansätze. This result is expected because the two types of excitation evolutions perform similarly in constructing electronic wavefunction ansätze. Since qubit-excitation evolutions are implemented by simpler circuits than fermionic-excitation evolutions, the QEB-ADAPT-VQE systematically outperforms the fermionic-ADAPT-VQE in terms of CNOT count in Fig. 6g–i.

While the QEB-ADAPT-VQE and the fermionic-ADAPT-VQE require similar numbers of iterations (Fig. 6a–c), the QEB-ADAPT-VQE requires up to twice as many variational parameters (Fig. 6d–f). This difference is due to the fact that the QEB-ADAPT-VQE assigns one parameter to each qubit-excitation evolutions in its ansatz, whereas the fermionic-ADAPT-VQE assigns one parameter to a pair of spin-complement fermionic-excitation evolutions.

Figure 6a–d shows that the QEB-ADAPT-VQE converges faster, requiring systematically fewer iterations and variational parameters than the qubit-ADAPT-VQE. As suggested in the section “Ansatz elements”, this result is due to the fact that single- and double-qubit-excitation evolutions correspond to combinations of 2 and 8 Pauli string exponentials.

In terms of CNOT count (Fig. 6g–i), the qubit-ADAPT-VQE is more efficient than the QEB-ADAPT-VQE at low accuracies. However, for higher accuracies, and correspondingly larger ansätze, the QEB-VQE-ADAPT starts to systematically outperform the qubit-ADAPT-VQE in terms of CNOT efficiency. This result can be attributed to the fact that qubit evolutions allow for the local circuit optimizations introduced in ref.  47, whereas Pauli string evolutions, albeit more variationally flexible, do not allow for any local circuit optimizations.

As a side point, it is interesting to note that when the fermionic-ADAPT-VQE is performed with a pool of independent single- and double-fermionic evolutions (Figs. 3 and 4) it is able to converge, albeit more slowly, to higher final accuracies than when it is performed with a pool of spin-complement pairs of single and double-fermionic evolutions (Fig. 6). This is owing to the fact that the pool of independent fermionic-excitation is more variationally flexible.

Discussion

In this work, we investigated the use of qubit excitations to construct electronic VQE ansätze. We demonstrated numerically that in general an ansatz of qubit-excitation evolutions can approximate a molecular electronic wavefunction almost as accurately as an ansatz of fermionic-excitation evolutions. However, fermionic-excitation-based ansätze were found to be a slightly more accurate per number of excitation evolutions when approximating strongly correlated states. These results suggest that, on their own, the Pauli-z strings, which measure the parity of the state and account for the anticommutation of the fermionic-excitation operators, play little role in the variational flexibility of an electronic wavefunction ansatz. These results agree with previous findings in refs.  37,45. Another advantage of fermionic-excitation evolutions is that they can form spin-complement pairs of fermionic-excitation evolutions. Such spin-complement pairs can then be used to enforce parity conservation and thus reduce the number of variational parameters of an ansatz by up to a factor of 2. Nonetheless, fermionic-excitation evolutions are implemented by circuits whose size, in terms of CNOT count, scales linearly (logarithmically) in the Jordan–Wigner (Bravyi–Kitaev) encoding with the system size, as opposed to qubit-excitation evolutions, which enjoy the quantum-computational benefit of being implemented by fixed-size circuits. Therefore, for NISQ devices, where the number of CNOTs is a primary cost factor, qubit-excitation evolutions are more suitable for constructing electronic ansätze.

Motivated by the accuracy and circuit efficiency of qubit-excitation-based ansätze, we introduce the qubit-excitation-based adaptive variational quantum eigensolver (QEB-ADAPT-VQE). The QEB-ADAPT-VQE simulates molecular electronic ground states with a problem-tailored ansatz, grown iteratively by appending single and double-qubit-excitation evolutions. We benchmarked the performance of the QEB-ADAPT-VQE with classical numerical simulations for LiH, H6, and BeH2. In particular, we compared the QEB-ADAPT-VQE to the original fermionic-ADAPT-VQE, and its more slowly converging, but a more circuit-efficient cousin, the qubit-ADAPT-VQE. Compared to the fermionic-ADAPT-VQE, the QEB-ADAPT-VQE requires up to twice as many variational parameters. However, the QEB-ADAPT-VQE requires asymptotically fewer CNOTs, owing to its use of qubit-excitation evolutions.

The simulations also showed that the qubit-ADAPT-VQE is more CNOT-efficient than the QEB-ADAPT-VQE in achieving low accuracies that correspond to small ansatz circuits. However, for higher accuracies and correspondingly larger ansatz circuits, the QEB-ADAPT-VQE systematically outperformed the qubit-ADAPT-VQE in terms of CNOT efficiency. The primary reason for this is that qubit evolutions allow for local circuit optimizations, while the more rudimentary Pauli string evolutions, utilized by the qubit-ADAPT-VQE, do not. In practice, we are often just interested in reaching chemical accuracy. Therefore, one might question what is the usefulness of constructing more CNOT-efficient ansätze with the QEB-ADAPT-VQE for accuracy higher than chemical accuracy. Although the numerical results presented here are not sufficient to draw a general conclusion, they indicate that the CNOT efficiency of the QEB-ADAPT-VQE becomes more evident for larger ansatz circuits. Therefore, for larger molecules, the QEB-ADAPT-VQE will likely be able to reach chemical accuracy using fewer CNOTs than the qubit-ADAPT-VQE. Our simulation results also demonstrated that in terms of convergence speed, the QEB-ADAPT-VQE requires fewer variational parameters, and correspondingly fewer ansatz-constructing iterations, than the qubit-ADAPT-VQE.

These results imply that the QEB-ADAPT-VQE is more circuit-efficient and converges faster than the qubit-ADAPT-VQE, which to our knowledge was the previously most circuit-efficient, scalable VQE protocol for molecular modeling. We do remark though, that in our comparison of the QEB-ADAPT-VQE and the qubit-ADAPT-VQE, we ignored the fact that the latter protocol can use a reduced ansatz element of O(NMO) Pauli string evolutions, as shown in ref. 37. Using a reduced ansatz-element pool would decrease the number of required quantum computer measurements, but will also result in a slower and less efficient ansatz construction. Moreover, the complexity of a single iteration of both the QEB-VQE-ADAPT and the qubit-ADAPT-VQE, might actually be dominated by running the VQE (see Supplementary Note 4). Therefore, reducing the size of the ansatz-element pool might not affect the overall complexity of the protocol. We also note that, in theory, hardware-efficient ansätze and the ansätze of the IQCC protocol suggested in refs. 39,40 can be implemented by shallower circuits than the ansätze constructed by the QEB-ADAPT-VQE. However, hardware-efficient ansätze and the IQCC are unlikely to be scalable for large systems: the optimization of hardware-efficient ansätze is likely to become intractable for large systems; and the IQCC requires evaluating a number of expectation values, exponential in the number of variational parameters.

As further work, three potential upgrades to the QEB-VQE-ADAPT can be considered. First, the ansatz-element pool of the QEB-VQE-ADAPT can be expanded to include non-symmetry-preserving terms as suggested in ref. 64. Potentially, this expanded pool could further improve the speed of convergence and boost the resilience to symmetry-breaking errors of the QEB-VQE-ADAPT. Second, methods from ref. 41 can be used to “prune”, from the already constructed ansatz, qubit-excitation evolutions that have little contribution to the energy reduction. This could potentially optimize further the constructed ansatz. Third, the QEB-VQE-ADAPT functionality can be expanded to enable estimations of energies of low-lying excited states. This will be the topic of another work (see ref. 65 for a preprint).