Quantum state preparation of multiconfigurational states for quantum chemistry

Abstract

The ability to prepare states for quantum chemistry is a promising feature of quantum computers, and efficient techniques for chemical state preparation is an active area of research. In this paper, we implement and investigate two methods of quantum circuit preparation for multiconfigurational states for quantum chemical applications. It has previously been shown that controlled Givens rotations are universal for quantum chemistry. To prepare a selected linear combination of Slater determinants (represented as occupation number configurations) using Givens rotations, the gates that rotate between the reference and excited determinants need to be controlled on qubits outside the excitation (external controls), in general. We implement a method to automatically find the external controls required for utilizing Givens rotations to prepare multiconfigurational states on a quantum circuit. We compare this approach to an alternative technique that exploits the sparsity of the chemical state vector and find that the latter can outperform the method of externally controlled Givens rotations; highly reduced circuits can be obtained by taking advantage of the sparse nature (where the number of basis states is significantly less than 2^n_q for n_q qubits) of chemical wavefunctions. We demonstrate the benefits of these techniques in a range of applications, including the ground states of a strongly correlated molecule, matrix elements of the Q-SCEOM algorithm for excited states, as well as correlated initial states for a quantum subspace method based on quantum computed moments and quantum phase estimation.

---

1 Introduction

Quantum computation has the potential to make a large impact on quantum chemistry and condensed matter physics. This has motivated a recent rapid development of techniques for representing quantum chemical states on gate-based quantum computers.^1–12 From the point of view of general quantum state preparation, many research works have focused on maximizing the efficiency of loading classical data into a quantum processor and preparing a quantum state to represent those data.^13–18 It is interesting to consider applications of the latter approach to quantum chemistry, where the classical data represent previously obtained chemical information, and the quantum computer can be used to further evolve the quantum chemical state (for example, to increase the accuracy of the representation of the ground or excited states of a molecule).

The existence of gate sets that are universal for quantum chemistry (in the sense of spanning the space of all states that preserve fermionic symmetries) has recently been proved.³ These gates take the form of Givens rotations (GRs) to represent particle-number-conserving excitations that correlate the occupations of different fermionic modes. GRs have also been used to construct gate fabrics that can navigate fermionic Fock spaces.¹ Regarding the preparation of multiconfigurational states (where each configuration represents electrons distributed throughout a molecular orbital basis, in a second quantized framework), if the configurations are specified a priori then a sequence of GRs can be applied to the circuit to prepare the desired state vector.³ However, GRs will mix any basis states with qubit subspaces that match those involved in the rotation. To ensure that only the desired basis states are included in the state vector, GRs can be controlled by qubits outside the rotation space³ (which we refer to as “external” controls to distinguish these from internal CNOTs within the GR gate decomposition). However, externally controlling each GR on all qubits not directly involved in the rotation leads to very large circuits when decomposed into typical hardware gate sets. The question then remains as to how to apply external controls to guarantee the desired state vector is produced, while avoiding, when possible, externally controlling all GRs on all qubits outside the excitation. A specific example of this was presented in³ for a 6-qubit state consisting of 4 occupation number (ON) configurations. In this paper, we provide a general algorithm for finding external GR controls automatically for any particle-number-preserving fermionic state.

We also consider the previously proposed method of Gleinig and Hoefler to prepare sparse quantum states.¹⁹ Initially designed to overcome the problem of exponentially scaling resources for arbitrary quantum state preparation,^14,20,21 the sparse state preparation (SSP) method takes advantage of the fact that many interesting quantum states span a relatively small section of the full Hilbert space, where the number of non-zero basis coefficients is far less than O(2ⁿ). This is the case for typical states of interest in quantum chemistry; the full configuration interaction (FCI)²² wavefunction obeys fundamental fermionic symmetries (e.g. spin and particle number conservation) which rule out many basis states of the entire Hilbert space. Additionally, the success of methods such as semistochastic heat-bath configuration interaction (SHCI)²³ and perturbatively selected configuration interaction (CIPSI)²⁴ exemplify the fact that a relatively small number of configurations can often provide a good approximation to the full wavefunction.

The importance of the sparsity of chemical wavefunctions is emphasized when considering approximations to the ground state, particularly in situations requiring “warm start” states, i.e. those states with large overlap with the true ground state (relative to single configuration Hartree–Fock (HF) states) yet remain sparse enough to be efficiently prepared. Considering the example of quantum phase estimation (QPE) in which the success probability is proportional to the overlap between the initial state and the true ground state, and the fact that in many interesting cases (e.g. strongly correlated systems) the HF state has significantly low overlap, the ability to efficiently and conveniently prepare multiconfigurational states can be highly beneficial. A recently published work²⁵ noted the method of Gleinig and Hoefler as a non-variational approach to quantum chemical state preparation, and provided a comparison to other sparse state preparation approaches in terms of asymptotic scaling. Here, an implementation of Gleinig and Hoefler's approach suitable for quantum chemistry, along with its integration with various ground and excited state methods, is reported and applied to a strongly correlated chemical system.

In this paper we compare the techniques of externally controlled GRs³ and Gleinig and Hoefler's scheme¹⁹ for preparing quantum circuits corresponding to selected linear combinations of computational basis states, where the latter represent occupation number configurations (Slater determinants). These techniques are implemented in the InQuanto^26–28 software package, and we demonstrate their utility for a wide range of quantum algorithms useful for quantum chemistry.

2 Preparation of multiconfigurational states in quantum circuits

2.1 Givens rotations with external controls

Consider an arbitrary single electronic configuration specified by an ON configuration²² in a finite basis of spin orbitals |…, f_i, …, f_j, …〉 with fermionic occupation numbers f_i, f_j ∈ {0, 1}, where f_i ≠ f_j, and i, j label orbitals involved in the excitation. An excitation can be viewed as a mixing between ON configurations with different occupations. Adopting the Jordan–Wigner (JW) transformation between fermionic and qubit algebras^29,30 |…, f_i, …, f_j, …〉↦|…, q_i, …, q_j, …〉 where q_i, q_j ∈ {0, 1} are individual qubit states, a fermionic excitation can be encoded as a unitary operation applied to a register of qubits representing the ON configuration, rotating it to a convex combinationHere, is a fermionic excitation operator, and is a unitary operator (encoded as quantum gates) applied to the qubit register. Note that should preserve the total number of electrons (the qubit register Hamming weight).

Eqn (1) corresponds to a 1-body excitation. Using a 2-qubit state as a minimal example, the excitation can be expressed in the computational basis as a GR,^1–3 whose action on the 2^n_q=2 dimensional Hilbert space is to rotate in a 2 dimensional subspace of the state vector

The coefficients of the orthonormal basis are omitted for the sake of clarity (the state vector is normalized to unity). The real scalars are restricted by |c₁²| + |c₂²| = |c₃²| + |c₄²| = 1 with c₂ and c₃ having opposite signs, to ensure unitarity. This GR can be written aswhere θ parameterizes the quantum gates implementing the Givens rotation.^1,3 Eqn (2) and (3) show, in a small example, that individual GRs preserve the particle number within the subspace to which they are applied. Generalizations up to n_elec-body excitations (with n_elec the number of electrons) for n_q ≥ 2n_elec qubits are obtained by extending the mixing elements to other entries of the 2^n_q × 2^n_q matrix (explicit matrices for 2-body rotations on 4 qubits are given in ref. 1).

The Hamming distance between two binary strings x and y, defined as

plays a central role in the use of GRs to prepare multiconfigurational states. is related to the excitation order for λ-body excitations between two ON configurations. The qubits on which the excitation unitary acts can be specified using and its decomposition into an array of qubit-wise XOR values (x(0) ⊕ y(0), x(1) ⊕ y(1), …, x(n_q) ⊕ y(n_q)). The indexes of non-zero elements of this array correspond to the GR qubit indexes, while the value of specifies the excitation order.

Selected linear combinations of ON configurations with fixed coefficients can be prepared using sequences of GRs applied to the circuit. Each GR mixes a pair of basis states (i.e. excitation) and its rotation angle θ can be obtained from state vector coefficients via recursive normalization (see Arrazola et al.³). Alternatively, optimized coefficients can be found by varying the gate angles within each rotation in a variational algorithm. As mentioned in the introduction (and elaborated in previous work³), the gates that encode each GR may need to be externally controlled to ensure that the sequence of GRs does not yield unwanted basis states. Thus, all the steps necessary to prepare a linear combination of ON configurations on a quantum circuit are at hand, apart from a general method to specify the external controls.

In the following subsection, we outline a procedure to automatically accomplish this task given a selected set of ON configurations, which specifies the excitation unitaries and the external controls needed to accomplish each GR. This can be used to build the circuit representing the multiconfigurational state.

2.1.1 Algorithm to specify givens rotations and external controls. In Algorithm 1, a procedure is presented to find a sequence of GRs and their external controls required to prepare a selected linear combination of ON configurations. The input is an ordered set of D bit strings representing the ON configurations. Note that the reference for excitations is the first member of the input set, and the desired state vector is built up as sequential excitations relative to this reference. The ordered D – 1 Hamming distances and their qubit components are then stored. Gate angles of each GR, dependent on the ON coefficients, can be found by using recursive normalization of ON coefficients (a specific example is also discussed in Section 3.1.1).

Let represent a rotation that linearly mixes the reference with a basis state indexed by e = 2, …, D. The excitation indexes are obtained from a decomposition of the Hamming distance (non-zero contributions to the sum in eqn (4)), which specifies each GR unitary in the (D – 1)-length sequence, i.e. , where is a tuple of indexes of length 2λ corresponding to nonzero contributions to . Note that for the basis state index e = 2, external controls are not needed by definition.

For a given ordered set of input configurations, the search for external controls required for depends on whether a “previous” (p < e) basis state is “rotatable” by , and if so controls are applied to ensure that will not act on any basis state with index p = 2, …, e – 1.

For basis states x_e separated from the reference x₁ by or , Givens rotations or , respectively, can be used to mix x_e with x₁, and the external controls found by Algorithm 1 ensure that the circuit maintains the desired state vector. For and , we use the gate decompositions presented in previous works for 1-body and 2-body GRs^1,3 (compilations of these circuits in the H-series gate set³¹ are given in the Appendix Fig. 10). However, when , a unitary is required that encodes the (λ > 2)-body excitation. To this end, we adapt a gadget described by Arrazola et al.³ consisting of (a) a sequence of multi-controlled SWAPs, (b) a central , and (c) the reverse of (a). In the Appendix section A.1, Algorithm 2 describes how the indexes of these unitaries and their external controls are found. Compared to Arrazola et al.,³ our scheme achieves a lower overhead of external controls for SWAPs and the central to combine ON configurations separated by λ > 2. An example of this for a 3-body excitation of the |111000〉 configuration is shown in Fig. 12.

Once the for all x_e are found, along with their external controls (if required), the desired state |ψ〉 can be prepared as

2.2 Sparse state preparation for quantum chemistry

In this section, we briefly describe Gleinig and Hoefler's sparse state preparation (SSP) method. For more details, we refer the reader to the original article.¹⁹ Here, we provide a brief description for completeness and highlight those aspects most relevant to multiconfigurational state preparation for quantum chemistry. This method essentially solves the problem of finding some n-qubit circuit such that , where |ψ〉 is the desired n-qubit state to be prepared. Knowing , one can reverse this to obtain the state preparation unitary such that the desired state is . The input to the SSP method is an unordered set of bit strings representing ON configurations, with normalized coefficients c_d such that the state to be prepared is and .

In the ith iteration of the algorithm,¹⁹ an arbitrary angle (possibly controlled) 1-qubit rotation gate along with an arrangement of NOT and CNOT gates (collectively labeled here as ) are found, which if applied to the state at the ith iteration results in a new state hosting < D_i−1 non-zero basis coefficients (where ). This is accomplished by “merging” two bit strings x_1,i, x_2,i contained in , where merging amounts to applying such that (up to global phase). The arrangement of gates in each are determined by the distribution of binary values throughout the D bit strings.¹⁹ We note that for real-valued c_d the R_y(θ) gate parameterized by angle θ is sufficient to represent the 1-qubit rotation in each .

The gate angles of the SSP method, which determine the values of c_d, are obtained on-the-fly and enter in the 1-qubit rotations, each of which mixes 2 basis states x_1,i, x_2,i at a given iteration i of the SSP algorithm; the ith rotation angle is obtained by normalizing the x_2,i coefficient c_2,i to and taking the arcsin. For the next i + 1 iteration, the coefficient c_2,i is replaced by the previous normalizer , and the angle of the (i + 1)^th single qubit rotation is obtained for x_1,i+1, x_2,i+1 as in the previous iteration. For an example of this, see Gleinig and Hoefler's paper.¹⁹

An important consequence of this algorithm is that the sequence of applied to the initial circuit |0^⊗n〉 is independent of the order of bit strings provided by the user at input. This is a major difference with the GR method described in Section 2.1, and highlights the lack of a “reference” ON configuration for excitations in the SSP method (whereas in the GR method the excitation reference is fixed as the first member of the ordered set of input bit strings, and all GRs subsequently applied can be considered as rotations relative to the reference). While this does not affect the behavior of SSP-prepared circuits for fixed (non-variational) c_d, variational applications of this method (in which gate angles are variationally optimized, hence leading to variationally optimized c_d) can be affected, as discussed in Section 3.1.1.

2.3 Methods which benefit from multiconfigurational state preparation

In the following subsections we mention a non-exhaustive list of methods which we use to demonstrate the utility of circuit preparation to represent selected linear combinations of ON configurations. We also mention previous applications, namely the preparation of N ± 1-particle states for Green's functions,³² and quantum circuit representations of CASSCF (complete active space self-consistent field) wavefunctions to study orbital entanglement.³³

2.3.1 Variational quantum eigensolver. The variational quantum eigensolver (VQE)⁷ is a hybrid quantum-classical algorithm. As initially proposed, VQE utilizes a quantum computer to perform measurements of the expectation value of a Hamiltonian taken with respect to a parameterized ansatz for the wavefunction, which is optimized by a classical processor based on a minimization of the energy according to the variational principle.

On present-day noisy intermediate scale quantum (NISQ) hardware, the originally proposed form of VQE is significantly limited by various sources of noise (e.g. qubit decoherence, gate error rates, state preparation and measurement errors). Hence, in this work, VQE is utilized as a classical optimizer of gate angles in multiconfigurational state circuits, which directly translates into optimization of basis state (ON configuration) coefficients of the circuit state vector. This demonstrates how the methods for preparing multiconfigurational state circuits (Section 2.1 and Section 2.2) can be used to prepare variational ansatzes for VQE applied to quantum chemistry. In this context, VQE is run in an ideal setting in which quantum measurements are classically simulated so that consists of ideally optimized basis coefficients. can then be used as the initial (warm start) state in subsequent methods (such as quantum subspace methods, or quantum phase estimation). One can also measure expectation values of the Hamiltonian in a noisy setting using offline optimized to quantum compute estimates of the ground state energy corresponding to ideal VQE parameters (see Section 3.1.1).

2.3.2 Quantum computed moments. A wide range of quantum subspace methods currently exist,³⁴ which generally attempt to solve for accurate approximations to low-lying eigenvalues of a given Hamiltonian Ĥ by projecting onto a finite-dimensional subspace spanned by a set of basis states. Krylov basis states correspond to a subspace spanned by applying powers of the Hamiltonian to a trial state (reference), and various forms of Krylov subspace methods have recently been investigated (see a recent review³⁴ and references therein). The quantum computed moments (QCM) approach, recently introduced by Vallury et al.³⁵ and subsequently applied in various contexts,^32,36 also utilizes powers (or moments) of the Hamiltonian Ĥ. While essentially constituting a subspace method, the QCM approach does not explicitly construct the Krylov basis states on a quantum circuit. Rather, QCMs provide a correction to the energy that encodes electronic correlation beyond that included in 〈Ĥ〉 (where 〈〉 denotes the expectation value with respect to an input state that approximates the exact ground state). Based on the cumulant expansion of the Lanczos diagonalization method, the elements of the tridiagonalized representation of Ĥ can be expressed in terms of polynomials of 〈Ĥ_m〉, where the maximum value of m is related to the dimension l of the subspace expansion as 2l – 1.

By deriving a truncation of the Lanczos expansion at m = 4,³⁷ a non-perturbative approximation to the ground state energy referred to as QCM4 can be expressed using cumulants (also referred to as connected moments)

In this work, we also utilize the connected moments expansion^38–41 in which the “t-expansion” of Horn and Weinstein⁴² can be expressed in terms of . This leads to alternative expressions of the energy to various orders of moments. Here we consider the CMX2 formula, in which “2” refers to the dimension of the associated subspace l

hence this can be considered as an energy correction up to m = 2l – 1 = 3. We compare CMX2 to QCM4 as a lower-cost (generally less measurements due to lower-order moments), yet more approximate form of the energy.

To obtain expectation values 〈Ĥ_m〉, the fermionic Hamiltonian is JW encoded as (where Pauli strings P̂_r,m are tensor products of Pauli gates spanning the qubit register and h_r,m are real coefficients), then Ĥ_m > 1 are obtained by recursive multiplication of Ĥ. 〈Ĥ_m〉 can then be computed from quantum hardware measurements in the computational basis via Pauli string averaging. For ideal simulations, expectation values are computed classically (considered as the limit of infinite sampling in the absence of device noise).

We note that the use of multiconfigurational (or multireference) initial states for connected moment expansions of the energy has a long history in classical computational chemistry.³⁸ In Section 3.1.2 we show a quantum computational version of this approach, where the circuit corresponding to the initial state is tailored to linear combinations of specific ON configurations (Slater determinants).

2.3.3 Quantum phase estimation and Hamiltonian simulation. Quantum Phase Estimation (QPE) extracts ground or excited state energies when the prepared initial state has a large overlap with the desired target eigenstate. A promising alternative to canonical QPE uses Hamiltonian simulation with a Hadamard test to measure the complex overlap between the initial state and the time-evolved states at several time steps. The resulting time-series is then post-processed, e.g. by quantum complex exponential least-squares (QCELS),⁴³ to extract the energies. Similarly to the canonical QPE, the performance of the Hamiltonian simulation with QCELS depends sensitively on the initial state, therefore it is of interest to prepare the initial state in the most resource efficient way. Formally, the QCELS method requires the measurement of

Z(t) = 〈ψ|e^−itĤ|ψ〉, (9)

where if |ψ〉 = |ψ₀〉 is the exact ground state, (or an eigenstate), then it coincides with e^{−it〈ψ₀|Ĥ|ψ₀〉}, from which the energy E₀ = 〈ψ₀|Ĥ|ψ₀〉 can be easily extracted. Typically, however |ψ〉 is not an eigenstate. Following the QCELS procedure we sample the time t_n = nτ with n = 0, …, N – 1 and time step τ and obtain an estimate of Z(t_n) via the Hadamard test. We preparewhere U(τ) = e^−iτĤ, and the real and imaginary parts are estimated via the expectation values

Re Z_n = 〈Ψ(t_n)|X⊗I|Ψ(t_n)〉, (11)

and

Im Z_n = 〈Ψ(t_n)|Y⊗I|Ψ(t_n)〉. (12)

The energy is then extracted from the time-series by maximizing the QCELS objectiveand the location of the maximum yields the energy estimation. The time step is chosen such that , where E_max and E_min are the largest and smallest eigenvalues of the simulated Hamiltonian. This choice ensures τE ∈ (–π, π) for every eigenvalue E, therefore avoiding phase aliasing and enabling unambiguous extraction via eqn (13). In practice, the spectrum can be centered by using , so that the ground state energy corresponds to the smallest eigenphase.

In Section 3 we demonstrate how multiconfigurational state preparation can increase the overlap with the ground state and therefore impacts the performance of the QCELS method.

2.3.4 Quantum self-consistent equation of motion. The calculation of excited states has attracted the interest of the quantum computing community as a first step towards the calculation of response properties. Quantum self-consistent equation of motion (Q-SCEOM) was suggested recently for the calculation of excitation energies⁴⁴ and response properties.⁴⁵ Q-SCEOM uses self-consistent excitation operators in order to satisfy the vacuum annihilation condition even when an approximate set of excitation operators is used.^44,46,47 As a result of this, Q-SCEOM yields accurate results in cases where plain qEOM⁴⁸ fails qualitatively (see recent works^44,45,49,50 for comparisons between Q-SCEOM, qEOM, and quantum subspace expansion^51,52). As shown by Asthana et al.,⁴⁴ the use of self-consistent excitation operators simplifies the final EOM equations (cf. qEOM). The excitation energies are obtained by diagonalization of the M matrix whose elements are given by (see Asthana et al.⁴⁴ for a detailed derivation):where Ĝ_I are excitation operators, Ĥ is the system's Hamiltonian, and |ψ_HF〉 is the Hartree–Fock (HF) state. is a unitary ansatz with angles optimized for the ground state. A simplified form of the off-diagonal elements can be obtained⁴⁴whereAs shown in eqn (16), the calculation of the off-diagonal elements of the M matrix requires the construction of a multiconfigurational state, |ψ_I+J〉 = |(Ĝ_i + Ĝ_j)|ψ_HF〉 which is the most demanding step for the quantum computation part.

Here we discuss the application of the multiconfigurational state preparation methods presented in Section 2, comparing both approaches, each being a building block of the Q-SCEOM implementation in InQuanto.^26,28

2.4 Computational details

To demonstrate the use of multiconfigurational state preparation as a useful approach to prepare initial states for strongly correlated systems, we study the ground and excited states of twisted C₂H₄ in the minimal basis set STO-3G. The ground state corresponds to a singlet at torsion angles close to 0° or 180°, while at 90° the ground state is quasi-degenerate with the triplet having a slightly lower energy. In terms of classical preprocessing, the PySCF code⁵³ is used to build the fermionic (second quantized) Hamiltonian and to run complete active space configuration interaction (CASCI) and unrestricted HF calculations. Active spaces are selected based on the ordering of mutual information⁵⁴ between orbital pairs (see Fig. 15 in the Appendix for the plots of the molecular orbitals (MOs) that were used in the active spaces). For all active spaces and all torsion angles, the natural orbitals are dominated by the carbon 2p shell, with small admixtures of H 1s that grow (but remain minority) as the torsion angle tends to 90°. For larger (8-qubit and 12-qubit) active spaces, higher orbitals that contribute to correlation in the ground state also contain mixtures of 2s and 2p. Fermionic Hamiltonians are JW transformed to obtain Ĥ as a sum of Pauli strings, for each active space. For VQE optimizations, the L-BFGS-B method as implemented in the SciPy package⁵⁵ was used. The multiconfigurational state preparation methods are implemented in the InQuanto software package,^26–28 and all quantum calculations are carried out using InQuanto. Quantum circuits are compiled using the architecture agnostic software compiler TKET.⁵⁶

3 Results

3.1 Ground state of C₂H₄

Here we demonstrate the use of multiconfigurational state circuits to calculate the ground state energy of twisted C₂H₄, using variational and non-variational techniques. A common bottleneck in these eigensolvers is the large overhead in measurements needed to accurately represent the required expectation values (typically by averaging over Pauli string measurements). This issue is exacerbated in strongly correlated systems in which a single reference (single ON configuration) state does not have a large overlap with the true ground state. A warm start initial state can be highly beneficial in this regard, particularly for quantum subspace methods, as the measurement overhead can be reduced, for example by requiring a smaller number of Hamiltonian moments or a smaller dimensionality of the subspace.

Before discussing results for the ground state energies of specific active spaces, we mention a popular quantum computational method to approximate the ground state which is Unitary Coupled Cluster Singles and Doubles (UCCSD),¹² with excitation parameters θ_k optimized by VQE. This ansatz can be expressed as a series of unitary operators applied to a reference state (e.g. to the Hartree–Fock state |ψ_HF〉)

where the fermionic excitation operators T̂_k consist of single, double, triple, …, excitations, which excite 1, 2, 3, …, fermions between occupied and virtual orbitals. In UCCSD, the excitations are truncated to singles and doubles, corresponding to each T̂_k operating over a maximum four spin oribtal indexes. (In practise, a Trotterized form of UCCSD is used as the variational ansatz, since individual excitations may not commute. Here we omit the Trotterization from eqn (17) for brevity).

If one applies the second quantized fermionic operators T_k directly to |ψ_HF〉, a series of ON configurations can be generated which correspond to all possible single and double excitations relative to the HF state. The latter in fact correspond to the basis states of Configuration Interaction Singles and Doubles (CISD),²² and the CISD wavefunction is obtained once the coefficients of those basis states are found (typically by diagonalization of the Hamiltonian in the CISD basis). While UCCSD can be more accurate than CISD (due to size extensivity, for example), CISD contains a significant portion of electronic correlation and often yields reasonable approximations to the ground state. Due to the framework we developed for preparing multiconfigurational states, a variational circuit corresponding to the CI expansion can easily be generated once the ON configurations of the expansion are given. In Fig. 1, we compare CISD circuit sizes prepared using the SSP and GR methods to UCCSD circuits (with all gate angles represented symbolically) for a range of active spaces. For each active space, the ON configurations are obtained by applying T̂_k to the closed shell singlet HF state. Therefore, the states prepared using these ON configurations retain the same spin symmetry as the ground states of C₂H₄ for torsion angles close to 0° or 180°.

Fig. 1 Circuit resources obtained for UCCSD and CISD, where the latter are prepared using the GR (Section 2.1) or the SSP¹⁹ (Section 2.2) methods. Number of qubits n_q ∈ {4, 6, 8, 10, 12, 14, 16}. For a given n_q (equal to number of spin orbitals), the number of electrons is if is even or if is odd, for which the HF reference is a closed shell singlet configuration. For a given n_q and number of electrons, the number of configurations corresponds to the number of single and double excitations plus the HF reference. Circuits are compiled to the standard gate set using the Qiskit⁵⁷ extension of TKET.⁵⁶

We observe that the SSP method produces circuits significantly smaller than the GR method over the range of active spaces chosen (4–16 qubits correspond to 2–8 active molecular orbitals). This is primarily due to the scaling in the number of external controls required for the GR method as the state becomes more complex with size, in addition to the decomposition of GR 1-body and 2-body rotations. The circuits produced by SSP are also consistently smaller than those produced by UCCSD for the same active space, implying the lower susceptibility to gate fidelity errors and qubit decoherence of the SSP CISD circuits, despite the generally lower accuracy of CISD compared to that of UCCSD.

3.1.1 4-qubit active space: VQE. As a first demonstration of multiconfigurational state preparation methods, we consider ansatzes for VQE applied to twisted C₂H₄ in an active of space n_MO = 2 molecular orbitals (MOs), corresponding to 4 qubits in the JW representation, populated by 2 electrons. Representing the ON configurations as |q₁,q₂,…,q_nq〉 where n_q = 2n_MO, the HF state can be written in a bit string representation as |1100〉. The ground state at equilibrium is a closed shell singlet, which transitions to a triplet state 〈S²〉 = 2 at 90° torsion angle. Hence, to capture the lowest energy manifold throughout the 180° rotation the following ON configurations are selected and included in the variational ansatz for VQE optimizationThe relation between variational parameters and coefficients c_d are described as follows.

For the GR method, as described in Section 2.1 each GR linearly combines the reference (x₁) with another ON configuration in the input set. 3 GRs are required to linearly combine the configuration pairs . In this case, externally controlling on the second qubit (q₂) is required to maintain the desired state vector (see Fig. 2): since x₂ = |1001〉 is “rotatable” by , then the external control is placed so that only operates on basis states in which q₂ = 1 (and hence has no effect on |1001〉). The corresponding pairs of coefficients (elements of the matrices related to rotation angle, see Eq. (3)) are obtained through recursive normalization of the c_d state vector coefficients,³ which here become where and with α₀ = 1†. For the SSP method, coefficients are related to gate angles as described in sec. 2.2 and exemplified in Gleinig and Hoefler's paper.¹⁹

Fig. 2 Circuit corresponding to state c₁|1100〉 + c₂|1001〉 + c₃|0110〉 + c₄|0011〉 prepared using the GR method (see Section 2.1). Algorithm 1 found that externally controlling on the second qubit is required. Substituting the gate parameters for 90° torsion as an example and compiling to the H-series gate set,³¹ this circuit can be represented using 61 PhasedX gates, 4 Rz gates, and 44 2-qubit ZZMax gates.

Fig. 2 and 3 show the VQE ansatz circuits of 4 ON configurations (Eq. (18)), built using the SSP and GR methods. Depending on the optimized VQE parameters , this ansatz spans both singlet and triplet eigenstates (with spin number restricted to s_z = 0). For example, at torsion angle = 0°, optimization of results in (up to global phase) a doubly excited closed shell singlet with negligible contributions from singly excited configurations (with |c₂|, |c₃| < 2 × 10⁻⁶). Whereas for torsion angle = 90°, with s_z = 0 the ground state is dominated by a superposition of open shell configurations, with small but non-negligible contributions from the closed shell configurations. (as all 4 ON configurations have non-negligible weight, this state vector is used as the example for circuit resources reported in the captions of Fig. 2 and 3). We note that following the optimization of all VQE parameters for torsion angles between 0° and 180°, ideal energies match the lowest eigenvalues obtained from exact diagonalization of Ĥ in the 2-particle sector (see Fig. 4) (and in the absence of device or measurement noise, ideal energies obtained from the GR and SSP methods are identical). Hence, the multiconfigurational state ansatz reproduces the ground state manifold for all torsion angles after VQE optimization.

Fig. 3 Circuit corresponding to state c₁|1100〉 + c₂|1001〉 + c₃|0110〉 + c₄|0011〉 prepared using the SSP method¹⁹ (see Section 2.2). Substituting the gate parameters for 90° torsion as an example and compiling to the H-series gate set,³¹ this circuit can be represented using 10 PhasedX gates, 4 Rz gates, and 5 2-qubit ZZMax gates.

Fig. 4 Energies obtained from VQE-optimized multiconfigurational states for the 4-qubit (2 electrons in 4 spin orbitals) active space of C₂H₄. Simulated measurements correspond to 10⁴ shots per circuit, and the Hamiltonian consists of 14 Pauli operators. Top graph: GR method (see Section 2.1). Bottom graph: SSP method (see Section 2.2). C₂H₄ structure at torsion angles 0°, 90°, and 180° shown above graphs. H11E corresponds to emulations of hardware experiments with a noise model calibrated to the H1 trapped ion device.³¹

Emulations of quantum hardware experiments use a noise model calibrated to the H1 trapped ion quantum computer,³¹ labeled as “H11E”. As shown in Fig. 4, quantum computations of the expectation value are subject to greater amounts of noise for the GR method due to the larger circuits (compare circuit resources reported in Fig. 2 and 3), related to the decompositions of particle-conserving GRs^1,3 in addition to the external control (for ) required for this case. Despite the negligible coefficients of singly excited configurations for some torsion angles ≠90°, we retain all 4 ON configurations in the circuit state vector for all torsion angles for ease of comparison. However, we note that shorter circuits (which represent only 2 ON configurations) for some torsion angles (at this active space) will likely lead to similar ideal energies.

Despite the smaller circuits produced by SSP for the same state, and hence its lower susceptibility to device errors, the GR method does exhibit a conceptual advantage over SSP, which can be particularly useful when interpreting a variational state in a chemical context: when all gate angles of a circuit produced by the GR method are 0, the state vector falls back to the first ON configuration of the input set (x₁ in (x₁, …, x_D)). If x₁ is chosen to be the HF state, then the GR unitaries added to this circuit can be readily interpreted as excitations on top of the HF reference. This is not necessarily the case for the SSP method, whose state is non-trivial to predict as it depends on the distribution of binary values throughout the set of input ON configuration bit strings.¹⁹ Hence, the reference state for chemical excitations is not necessarily accessible as the state of SSP method. This issue can be exacerbated in variational optimizations of the gate angles, which can terminate in local minima at small gate angles: for the SSP method this local minimum state can be unpredictable, whereas for the GR method it is likely close or equal to the HF state. For the small VQE applications presented here, this issue can be bypassed by random initialization of SSP parameters (avoiding initialization). However, difficulties may arise in larger systems with more complicated parameter landscapes.

3.1.2 8-qubit active space: Hamiltonian moments. For the 8-qubit active spaces, ON configurations are selected from the largest weight determinants observed in CASCI calculations. For torsional angles of 0° to 60° (and 120° to 180°), the ground state is dominated by the HF determinant (|11110000〉), consisting of mostly 2p orbitals in all MOs, with small mixtures of H 1s in the first (q₁, q₂) and forth (q₇, q₈) MOs. A paired double excited ON configuration (|11001100〉) is the next largest, followed by spin-paired open shell configurations (|10011001〉, |01100110〉). For torsion angle = 80°, the first two largest determinants are similar, while the open shell configurations instead have one MO (the third) fully occupied (|10011100〉, |01101100〉). At torsion angle = 90°, the ground state is a triplet, which is represented in its s_z = 0 component using 4 ON configurations, each hosting open shell MOs (|11100100〉, |11011000〉, |10110100〉, |01111000〉). 8-qubit variational ansatzes of 4 ON configurations are then prepared for all torsion angles, which are optimized using ideal VQE. The resulting ON configurations and their coefficients are shown in Table 1.

Table 1 State vectors and corresponding circuit resources for all torsion angles studied in the 8-qubit active space of C₂H₄. Note the symmetry around the 90° torsion angle (hence the circuits for 0° and 180° are equivalent, etc.). Circuit resources represent the GR method (SSP method in parenthesis). Circuits are compiled to the H-series gate set.³¹

Torsion angle (°)	c1|x1〉	c2|x2〉	c3|x3〉	c4|x4〉	PhasedX	Rz	ZZMax
0/180	0.9690|11110000〉	−0.2345|11001100〉	0.0546|10011001〉	0.0547|01100110〉	174 (22)	7 (6)	128 (17)
20/160	0.9683|11110000〉	−0.2380|11001100〉	0.0533|10011001〉	0.0534|01100110〉	174 (22)	7 (6)	128 (17)
40/140	0.9617|11110000〉	−0.2648|11001100〉	0.0503|10011001〉	0.0503|01100110〉	174 (22)	7 (6)	128 (17)
60/120	0.9354|11110000〉	−0.3481|11001100〉	0.0441|10011001〉	0.0441|01100110〉	174 (22)	7 (6)	128 (17)
80/100	0.8281|11110000〉	−0.5522|11001100〉	−0.0681|10011100〉	0.0681|01101100〉	66 (18)	4 (3)	40 (13)
90	0.7044|11100100〉	0.7044|11011000〉	0.0615|10110100〉	0.0615|01111000〉	52 (16)	4 (3)	32 (11)

---

The multiconfigurational states from Table 1 are then used as initial states for the QCM4 and CMX2 methods, with energies obtained from expectation values of Hamiltonian moments. Ideal results (noiseless and infinite shot limit) are shown in Fig. 5. We first note that for torsion angles of 0° to 60° and 120°–180° (representing relatively weak electronic correlation), the ideal energy values are reasonably accurate for the D = 1 HF input state, and highly accurate for the D = 4 multiconfigurational states.

Fig. 5 Ideal QCM4 and CMX2 energies calculated using a single configuration HF (top, D = 1) and VQE-optimized multiconfigurational input states (bottom, D = 4) for the 8-qubit (4 electrons in 8 spin orbitals) active space of C₂H₄. VQE energies correspond to expectation values of Ĥ taken with respect to the optimized multiconfigurational state. For ideal simulations, externally controlled GRs (Section 2.1) and the SSP method (Section 2.2) yield identical results.

At torsion angles near the transition point (80°, 90°, 100°, strongly correlated as the singlet and triplet eigenstates become quasi-degenerate) we find that the D = 1 closed shell HF input state results in numerically very small cumulants , leading to the term approaching 0 or negative and an ill-defined QCM4 formula when taking the square root (see eqn (7)). Noting that the original derivation assumed the condition,,³⁷ we omit the QCM4 energies from these angles for D = 1 (also considering that the closed shell HF state insufficiently represents the strongly correlated ground state, leading to a poor description of the correlation represented by the moment expectation values). However, for the D = 4 multiconfigurational states, both the QCM4 and CMX2 calculations recover the exact lowest energies for all torsion angles. This shows the benefit of multiconfigurational input states for methods involving connected moments, as not only does the QCM4 formula avoid the issue of negative (we assume due to a higher quality input state), but the lower order theory (CMX2) can obtain a similar accuracy, hence necessitating a lower order of moments and ultimately less Pauli strings to measure. The latter point is further emphasized by considering that the VQE energies at torsion angle near 90° are < 1mHa above the exact values; at these geometries the expectation of Ĥ with respect to the multiconfigurational state already captures most of the electronic correlation, obviating the need for higher order moments with a small overhead in circuit depth (see circuit resources reported in Table 1).

3.1.3 12-qubit active space: quantum phase estimation and Hamiltonian simulation. To quantitatively examine how the state preparation impacts the ground state estimation of QCELS, we performed exact state vector time evolutions for the 12-qubit (6e, 6o) active space of C₂H₄ at a torsion angle of 80°. At this geometry the Hartree-Fock configuration, which has the largest weight in the exact ground state |ψ₀〉, has fidelity only |〈ψ₀|HF〉|² = 0.674, which is below the 0.71 threshold recommended in ref. 43 for reliable QCELS phase extraction. Moreover, it is expected that the required run time, therefore the circuit depth of the time evolution operator, becomes shorter as the initial state approaches the exact ground state.⁴³ Therefore we performed simulations with initial states containing increasing number of ON configurations. All together four initial states were used in the simulations, denoted as |ϕ_ON-1〉 = |HF〉, |ϕ_ON-2〉, |ϕ_ON-4〉, and |ϕ_ON-8〉, where the last three are derived from CASCI and correspond to the 2, 4, and 8 largest weight ON configurations, respectively. We also note the fidelities of the multiconfigurational states: |〈ψ₀|ϕ_ON-2〉|² = 0.986, and |〈ψ₀|ϕ_ON-4〉|² = 0.994, |〈ψ₀|ϕ_ON-8〉|² = 0.999.

Fig. 6 summarizes the simulations. The top panel shows the complex time series Z(t_n) (real and imaginary parts) for each initial state. All curves except the HF case lie almost on top of each other. The middle panel plots the relative error with respect to the exact value, e^{−it_n〈ψ₀|H|ψ₀〉}. We chose τ = 0.1 Ha⁻¹ and N = 31, giving a total evolution time T_max = (N – 1)τ, long enough to capture at least one full period of the exact Z(t). The right panel displays the QCELS objective, whose maximum was located with the BFGS method. Finally, the bottom-left panel shows the energy errors obtained by repeating the simulations with smaller T_max while keeping N = 31, effectively choosing a smaller time step. The energies derived from the initial states |ϕ_ON-1〉 and |ϕ_ON-2〉 do not reach the typical target precision ε = 1mHa within the time window, even at the largest T_max, whereas those from |ϕ_ON-4〉 and especially |ϕ_ON-8〉 do, and importantly, the latter hits the target precision at approximately half the evolution time.

Fig. 6 Top left panel shows 〈ψ(0)|ψ(t)〉 real (solid) and imaginary parts (dotted) for different initial states, for Ĥ representing the 12 qubit (6 electrons in 12 spin orbitals) active space for torsion angle = 80°. Middle left panel shows |〈ψ(0)|ψ(t)〉 – Exact| for different initial states, exact overlap is calculated with ground state as initial state. Bottom left panel shows the QCELS energies for different initial states relative to the exact ground state energy. The x-axis is the time period for QCELS, demonstrating that the energy computed for that period requires T_max simulation time. As the period increases, the energy approaches target accuracy (dashed horizontal line) with different rates for different initial state. Legend shows the fidelity between the initial state and the exact ground state. Right panel shows QCELS phase function for the largest T_max period with the exact ground state energy marked by a dashed vertical line.

Because the depth of the circuit representation of e^−itH scales linearly with t, compact, high fidelity state representation can save significant amount of gate operations by reducing the evolution time. Preparing |ϕ_ON-4〉 via the SSP requires 13 two-qubit gates, while |ϕ_ON-8〉 needs 40. Although the state preparation thus costs 27 additional 2-qubit gates, the shorter evolution time is expected to result a net saving well beyond the 27 gates.

3.2 Excited states of C₂H₄

Consequently, we look at the excited states of C₂H₄ using the Q-SCEOM method. The UCCSD ansatz, optimized by VQE for the ground state, was used for in eqn (14). We also used particle and spin conserving singles and doubles excitation operators. We compare the efficiency of the GR and SSP methods for the construction of the off-diagonal elements of the M matrix (see Sec. 2.3.4 for details). In Fig. 7 we present the total number of gates (upper panel) and the total number of 2-qubit gates (lower panel) (the circuits were compiled with the H-series emulator) for each element of the M matrix for the C₂H₄ molecule using the 8-qubit active space with torsion angle of 90°. More specifically, the circuits correspond to for the diagonal elements and for the off-diagonal elements. There are 26 excitation particle and spin conserving operators. Therefore the M matrix consists of 676 elements. Only the diagonal and the upper-half part of the matrix were calculated since M is a symmetric matrix. The GR method was used for the plots in the left side while the SSP method was used for plots in the right side. It is evident that the SSP method reduces significantly the number of gates. This is in accordance to results presented in Section 3.1 for the calculation of ground state energies. We note here that in the case of the GR method, the structure of the matrices shown in Fig. 7 is directly analogous to the Hamming distance of the associated states to be combined. In particular, the improvement of the SSP method over the GR method becomes greater for values of the Hamming distance larger than 4, for which the GR method requires the construction exemplified in Fig. 12. (see Section A.2 in the Appendix for more details).

Fig. 7 (Top) total number of gates and (bottom) number of 2-qubit gates required for the construction states needed for the evaluation of the M matrix ( for the diagonal elements and for the off-diagonal elements) for C₂H₄ with 8-qubit active space and 90° angle. The (left) GR, and (right) SSP methods were used for the off-diagonal terms. The circuits were compiled with the H-series emulator.

Additionally, the gate decompositions of the multiconfiguration states for Hamming distances ≤4 also contribute to the reduced circuit sizes obtained from the SSP for the M matrix elements. This is exemplified in Fig. 8, where the GR and SSP circuits are shown for the construction of the |ψ_I+J〉 state for the M₂₀ element for the Q-SCEOM calculation of the 8-qubit C₂H₄ at 90° torsion. We observe that 3 2-qubit entangling gates are needed for the SSP method, whereas the GR method requires 14 (see also Fig. 10).

Fig. 8 Circuits to construct the multiconfigurational state |(Ĝ₂ + Ĝ₀)|ψ_HF〉 for the M₂₀ element. Green vertical bars represent the ZZMax 2-qubit entangling gate, while red horizontal bars represent PhasedX(α, β) = R_z(β)R_x(α)R_z(−β) 1-qubit rotations.³¹ The circuits represent (up to global phase) the state . All gate angles are in units of π. (a) GR method, for which the rotation decomposes into 14 ZZMax gates (see bottom panel of Fig. 10 for the representation of in the H-series gate set). (b) SSP method, containing only 3 ZZMax gates.

Next, we look at the ground state and excited energies of C₂H₄ at various torsion angles. We present in Fig. 9 the results for the 4-qubit and 8-qubit case in the upper and lower parts, respectively. In the left panel we plot in red the correlation energy of the ground state energy as the absolute difference of the HF energy with respect to the exact result (first eigenvalue obtained with diagonalization of the Hamiltonian). As expected, the degree of correlation increases as we approach the 90° torsion. We also plot in black the absolute difference of the HF energy with respect to the VQE result obtained through optimization of the UCCSD ansatz. It is evident that at 90° torsion angle, VQE yields a higher energy compared to the exact result. We attribute this to the fact that the VQE optimization results in an unstable singlet state whereas the stable solution is a triplet, as discussed in Section 3.1. Note here, that for all the VQE calculations reported in this section, prior to the calculation of the M matrix, we limited our search to open-shell singlet states. Consequently, we plot in the middle panel the VQE result for the ground state energy and the energy of the first three excited states obtained with Q-SCEOM. We employed the SSP method for the construction of the off-diagonal elements of the M matrix for C₂H₄ at various torsion angles. We compare the Q-SCEOM results to the eigenvalues obtained through diagonalization of the Hamiltonian. For both active spaces, Q-SCEOM reproduces qualitatively the exact result. Interestingly, at the 90° angle, Q-SCEOM recovers the correct total energy (which matches the exact result) even though VQE (which affects the Q-SCEOM result through the optimized ansatz that is used for the calculation of the M matrix) yields a higher energy.

Fig. 9 Noiseless calculations of ground and excited states for various torsion angles of 4-qubit and 8-qubit active spaces of C₂H₄ in the Upper and Lower panel respectively. The absolute difference of the HF ground state energy with respect to the VQE ground state energy () and with respect to the exact energy () are plotted in the Left part. The ground state energies (− ⋆ –) and the first (), second (), and third () excited state energies obtained with VQE + Q-SCEOM and exact diagonalization of the Hamiltonian are shown in the Middle and Right part respectively. The SSP method was used in the Q-SCEOM calculations.

4 Discussion

In this work, quantum circuit preparation of multiconfigurational states for quantum chemistry has been demonstrated and compared as the state preparation step for a range of applications, including ground state energy calculations using VQE, QCM, and QPE, as well as excited state energies calculated using Q-SCEOM. The use-case to demonstrate these techniques is twisted C₂H₄, whose energy surface passes through a strongly correlated point as a function of torsion angle.

Multiconfigurational state preparation allows initial states that are more accurate than the single configuration HF state, which for QCM can facilitate a reduction in measurement overhead at a relatively small cost of extra circuit depth for the initial state, whereas for QPE the probability of accurately measuring the phase can be increased by a multiconfigurational initial state (relative to HF) since this probability is proportional to the overlap between the initial state and true ground state.⁵⁸ For Q-SCEOM, the off-diagonal elements of the M matrix can be generated using multiconfigurational state preparation, hence the latter facilitates a framework for automatically building and running the Q-SCEOM algorithm for a given molecule.

Overall, we observe significantly more efficient circuits using the SSP method compared to the GR method. This is largely due to the utilization of the sparsity of chemical states that can be exploited by the SSP, in addition to the gate decomposition of particle-conserving GRs and their required external controls leading to larger gate and depth overheads for the GR method. A conceptual advantage of the GR method can be seen when considering the state preparation as transformations of a reference state; setting the rotation angle of all GRs to 0 is guaranteed to result in the state |x₁〉 where x₁ is the first ON configuration in the ordered input set. Since x₁ can be chosen from HF, each GR can be considered as an excitation of the HF reference, preserving familiar notions from classical computational chemistry. In the context of variational searches of the ground state energy, this property of the GR method prevents unpredictable local minima at small values of variational parameters (which the SSP method can be susceptible to), and may be beneficial for analyzing contributions of specific basis states to the total electronic correlation.

In terms of future applications, our implementation easily allows for the composition of different ansatz circuit structures. For example, appending generalized⁵⁹ UCC unitaries to a multiconfigurational state circuit prepared using the GR or SSP methods provides a framework for multireference methods⁶⁰ at the quantum circuit level, such as multireference UCC. We also note a recent work that extended the reference state error mitigation scheme⁶¹ to multireference states⁶² using GRs. An interesting future direction would be to apply the SSP method to prepare the multireference states used for error mitigation, with the potential of significant reduction in circuit resources.

In addition, since multiconfigurational state preparation as presented here can be seen as an approach of loading classical (chemical) data to a quantum processor and preparing a quantum state to represent this data, this could be compared to preparation of multi-determinant states using quantum read-only memory (QROM),⁶³ or to sparse state preparation using quantum random access memory.¹⁸ We note that while these techniques achieve a similar asymptotic scaling to the SSP method in the number of 2-qubit gates (O(Dn_q), our implementation of multiconfigurational state preparation does not require controlled operations between auxiliary qubits and state qubits. We also note a recent work⁶⁴ that utilizes QROM to reduce the Toffoli count for matrix product state preparation.

In conclusion, this work provides a framework for preparing quantum circuit representations of multiconfigurational states, as implemented in InQuanto.^26–28 Overall, the impact of this work can be summarized by the following points.

(i) We provide a novel implementation for the use of externally controlled GRs for multiconfigurational state preparation, which automatically finds the required external controls of GRs for a given ordered set of ON configurations.

(ii) While the SSP algorithm was published in a previous article,¹⁹ this work demonstrates novel applications of the SSP method for a range of algorithms currently proposed for studying the ground and excited states of molecules quantum computationally, i.e. VQE, QCM4, QPE, and Q-SCEOM.

(iii) Specific to the application of the SSP method to VQE, our implementation allows for the use of SSP to generate a variational ansatz, in which the lowest energy state vector can be obtained by optimizing the angles of the 1-qubit unitaries within the SSP circuit.

The utility of multiconfigurational initial states is shown for various quantum computational approaches to quantum chemistry. In particular, when the state to be prepared exhibits sparsity in the sense of D ≪ 2^n_q, very efficient circuit representations can be achieved, therefore boosting the accuracy of QPE, for example, or reducing the measurement overhead of quantum subspace methods, without a large overhead in circuit depth. In addition, multiconfigurational state preparation is useful for enabling certain methods (such as Q-SCEOM) which require circuit constructions of selected excitations of a reference state, thus enhancing the tool set for quantum approaches to ground and excited states of molecules.

Author contributions

All authors contributed to the conceptualization of the work. G. G. D, G. P., and D. Z. M. contributed to the implementation of methods and analysis of results. All authors contributed to the manuscript.

Conflicts of interest

The authors declare no conflict of interest.

Note added after first publication

This article replaces the version published on 12^th December 2025, which contained a typographical error in the set of basis states where should have been as well as inconsistent italicisation of eqn (10)–(12).

Data availability

The data and binary version of the software supporting the conclusions of this paper are available from the corresponding author upon request. Any provision of the software or data is at the discretion of the providing company. A trial version of InQuanto is available at https://www.quantinuum.com/products-solutions/inquanto-trial. InQuanto version 5.0 (https://inquanto.quantinuum.com/index.html) is used in this work. Pseudocode for the GR method is presented in this article (see Algorithms 1 and 2). See also pseudocode for the SSP method reported in the original paper by Gleinig and Hoefler (https://doi.org/10.1109/DAC18074.2021.9586240). A repository is available that contains representative starting datasets and notebooks to allow similar conclusions to be drawn to those of the manuscript based on the workflow in the notebooks and self-written replacements for calls to proprietary code (https://doi.org/10.5281/zenodo.17466867).

A Appendices

A.1 Circuit decompositions for GR method

In Fig. 10 we show the gate decompositions of the 1-body and 2-body rotations. Compared to the corresponding circuits shown in previous work,³ we observe a similar number of 2-qubit entangling gates (ZZMax³¹ here and CNOTs in Arrazola et al.³). The Hadamard gates of the decompositions shown in³ are absorbed into the PhasedX(α, β) = R_z(β)R_x(α)R_z(−β)³¹ 1-qubit rotations, and additional R_z and PhasedX rotations are required to achieve a 2-qubit operation equivalent to a CNOT, resulting in differences in the total number of 1-qubit gates. For completeness, the action of the 2-qubit ZZMax is depicted in Fig. 11 in terms of CNOTs and R_z.

Fig. 10 Circuits for 1-body (2-qubit) and 2-body (4-qubit) rotations, where state vectors are obtained up to global phase. Green vertical bars represent the ZZMax 2-qubit entangling gate, while red horizontal bars represent PhasedX(α, β) = R_z(β)R_x(α)R_z(−β) 1-qubit rotations.³¹ All gate angles are in units of π. Subscripts on and are omitted.

Fig. 11 The 2-qubit entangling gate ZZMax³¹ in terms of CNOT and R_z gates. The equality is up to global phase. Gate angles are in units of π.

Fig. 12 Circuit to perform a (λ = 3)-body excitation of the ON configuration |111000〉, yielding cos θ|111000〉 + sin θ|000111〉, where θ is in units of π. The subscript e in is omitted in the figure. See top panel of Fig. 10 for the decomposition of the central in the H-series gate set.³¹

This section also describes our scheme to build the excitation gadget that connects basis states separated by Hamming distance (λ > 2 – body excitation), when preparing the multiconfigurational state using the GR method. Unlike Fig. 5 of in Arrazola et al.,³ here each SWAP is not controlled on all other qubits in the circuit, but only on “minority” qubits, i.e. those qubits with minority binary values (see q_minor in Algorithm 2). This results in lower overhead of external controls that scale with the number of minority qubits n_minor rather than n_q as in Arrazola et al.³ However, this leads to the possibility that (i) the control states (conditions of the external control) of a SWAP are matched by the qubits of one of the previous basis states, or (ii) a swapped version of the first basis state becomes equal to one of the previous basis states (note that qubits of at a given iteration of the WHILE loop lines 5–13 of Algorithm 2 form the control states for the next iteration). If (i) or (ii) occurs, it is not necessarily a problem once the central is also controlled on the minority qubits of (“fully swapped”) . The latter implies that the external controls of the central also scale with n_minor, while a cheaper alternative would be to use the usual controls from Algorithm 1. However, in order for the external controls found by Algorithm 1 to be usable for the central of the λ > 2 excitation gadget, we find that both (i) and (ii) must not occur.

A.2 Circuit resources and hamming distances for the M matrix in Q-SCEOM

In this section, we plot the Hamming distances between the |ψ_I+J〉 and |ψ_HF〉 states, used to build the M matrix for Q-SCEOM. The Hamming distances have a significant impact on the circuit resources used for each element of the M matrix, particularly for the GR method. Fig. 13 corresponds to the 4-qubit active space, where the left and middle panels show the total number of gates for the GR and SSP methods, respectively, while the right panel shows the Hamming distances. Fig. 14 shows the Hamming distances for the 8-qubit active space.

Fig. 13 Total number of gates using the (left) GR method, and the (middle) SSP method for the construction of the ’ket’ states of the M matrix for the 4-qubit active space for C₂H₄ at 90 °torsion. (Right) the hamming distance between the I, J states that are combined using the previous methods to construct the |ψ_I + ψ_J〉 state. The Ĝ_I|ψ_HF〉 states for each element are shown as x tick labels.

Fig. 14 Hamming distance between the I, J states that are combined to construct the |ψ_I + ψ_J〉 states for C₂H₄ using the 8-qubit active space.

A.3 Plots of orbitals

In this section we plot the spatial orbitals used in the active spaces of C₂H₄, shown in Fig. 15.

Fig. 15 UHF spatial orbitals of C₂H₄ at 0° torsion angle that were used for the construction of the active spaces. The HOMO and LUMO shown in (a) were used for the 4-qubit case. The 8-qubit case was augmented with the orbitals shown in (b) and the 12-qubit case used all six orbitals shown in the figure. Occupied and virtual orbitals are shown in the upper and lower part respectively. The isosurface level was set to 0.02. PySCF⁵³ was used for the HF calculation.

Acknowledgements

We thank Alec Owens and Carlo Gaggioli for helpful feedback and comments. We also acknowledge Josh J. M. Kirsopp and Iakov Polyak for useful discussions.

References

1. G.-L. R. Anselmetti, D. Wierichs, C. Gogolin and R. M. Parrish, New J. Phys., 2021, 23, 113010.
2. C. H. Chee, D. Leykam, A. M. Mak and D. G. Angelakis, Shallow quantum circuits for efficient preparation of Slater determinants and correlated states on a quantum computer, 2023, arXiv:2301.07477,  DOI:10.48550/arXiv.2301.07477.
3. J. M. Arrazola, O. D. Matteo, N. Quesada, S. Jahangiri, A. Delgado and N. Killoran, Quantum, 2022, 6, 742.
4. S. Fomichev, K. Hejazi, M. S. Zini, M. Kiser, J. F. Morales, P. A. M. Casares, A. Delgado, J. Huh, A.-C. Voigt, J. E. Mueller and J. M. Arrazola, Initial state preparation for quantum chemistry on quantum computers, 2024, arXiv:2310.18410,  DOI:10.48550/arXiv.2310.18410.
5. Y. S. Yordanov, D. R. M. Arvidsson-Shukur and C. H. W. Barnes, Phys. Rev. A, 2020, 102, 062612.
6. Y. S. Yordanov, V. Armaos, C. H. W. Barnes and D. R. M. Arvidsson-Shukur, Commun. Phys., 2021, 4, 228.
7. J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth and J. Tennyson, Phys. Rep., 2022, 986, 1–128.
8. I. G. Ryabinkin, T.-C. Yen, S. N. Genin and A. F. Izmaylov, J. Chem. Theory Comput., 2018, 14, 6317.
9. F. A. Evangelista, G. K.-L. Chan and G. E. Scuseria, J. Chem. Phys., 2019, 151, 244112.
10. K. Sugisaki, S. Nakazawa, K. Toyota, K. Sato, D. Shiomi and T. Takui, ACS Cent. Sci., 2019, 5, 167.
11. H. R. Grimsley, D. Claudino, S. E. Economou, E. Barnes and N. J. Mayhall, J. Chem. Theory Comput., 2020, 16, 1–6.
12. A. Anand, P. Schleich, S. Alperin-Lea, P. W. K. Jensen, S. Sim, M. Díaz-Tinoco, J. S. Kottmann, M. Degroote, A. F. Izmaylov and A. Aspuru-Guzik, Chem. Soc. Rev., 2022, 51, 1659–1684.
13. L. Grover and T. Rudolph, Creating superpositions that correspond to efficiently integrable probability distributions, 2002, arXiv:0208112,  DOI:10.48550/arXiv.0208112.
14. M. Möttönen, J. J. Vartiainen, V. Begholm and M. M. Salomaa, Quantum Inf. Comput., 2005, 5, 467–473.
15. V. Shende, S. Bullock and I. Markov, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 2006, vol. 25, pp. 1000–1010.
16. I. F. Araujo, D. K. Park, F. Petruccione and A. J. da Silva, Sci. Rep., 2021, 11, 6329.
17. X.-M. Zhang, T. Li and X. Yuan, Phys. Rev. Lett., 2022, 129, 230504.
18. T. M. L. de Veras, L. D. da Silva and A. J. da Silva, Quantum Inf. Process., 2022, 21, 204.
19. N. Gleinig and T. Hoefler, 58th ACM/IEEE Design Automation Conference (DAC), 2021, pp. 433–438.
20. M. Plesch and i. c. v. Brukner, Phys. Rev. A, 2011, 83, 032302.
21. P. Niemann, R. Datta and R. Wille, IEEE 46th International Symposium on Multiple-Valued Logic (ISMVL), 2016, pp. 247–252.
22. T. Helgaker, P. Jørgensen and J. Olsen, in Second Quantization, John Wiley & Sons, Ltd, 2000, pp. 1–33.
23. J. Li, M. Otten, A. A. Holmes, S. Sharma and C. J. Umrigar, J. Chem. Phys., 2018, 149, 214110.
24. M. Dash, S. Moroni, A. Scemama and C. Filippi, J. Chem. Theory Comput., 2018, 14, 4176–4182.
25. C. Feniou, O. Adjoua, B. Claudon, J. Zylberman, E. Giner and J.-P. Piquemal, J. Phys. Chem. Lett., 2024, 15, 3197–3205.
26. A. Tranter, C. Di Paola, D. Muñoz Ramo, D. Z. Manrique, D. Gowland, E. Plekhanov, G. Greene-Diniz, G. Christopoulou, G. Prokopiou, H. D. J. Keen, I. Polyak, I. T. Khan, J. Pilipczuk, J. J. M. Kirsopp, K. Yamamoto, M. Tudorovskaya, M. Krompiec, M. Sze, N. Fitzpatrick, R. J. Anderson and V. Bhasker, InQuanto: Quantum Computational Chemistry, 2022, https://www.quantinuum.com/products-solutions/inquanto.
27. A. Tranter, C. Di Paola, D. Muñoz Ramo, D. Z. Manrique, D. Gowland, E. Plekhanov, G. Greene-Diniz, G. Christopoulou, G. Prokopiou, H. D. J. Keen, I. Polyak, I. T. Khan, J. Pilipczuk, J. J. M. Kirsopp, K. Yamamoto, K. Sorathia, M. Tudorovskaya, M. Krompiec, M. Sze, N. Fitzpatrick, O. Backhouse, R. J. Anderson and V. Bhasker, Introduction to the InQuanto Computational Chemistry Platform For Quantum Computers, 2022.
28. A. Tranter, C. Di Paola, D. Muñoz Ramo, D. Z. Manrique, D. Gowland, E. Plekhanov, G. Greene-Diniz, G. Christopoulou, G. Prokopiou, H. D. J. Keen, I. Polyak, I. T. Khan, J. Pilipczuk, J. J. M. Kirsopp, K. Yamamoto, K. Sorathia, M. Tudorovskaya, M. Krompiec, M. Sze, N. Fitzpatrick, O. Backhouse, R. J. Anderson and V. Bhasker, InQuanto User Guide, 2024, https://inquanto.quantinuum.com/index.html.
29. P. Jordan and E. Wigner, Z. Phys., 1928, 47, 631–651.
30. S. Szalay, Z. Zimborás, M. Máté, G. Barcza, C. Schilling and Ö. Legeza, J. Phys. A: Math. Theor., 2021, 54, 393001.
31. Quantinuum H1-1, 2025, https://docs.quantinuum.com/h-series/.
32. G. Greene-Diniz, D. Z. Manrique, K. Yamamoto, E. Plekhanov, N. Fitzpatrick, M. Krompiec, R. Sakuma and D. M. Ramo, Quantum, 2024, 8, 1383.
33. G. Greene-Diniz, C. N. Self, M. Krompiec, L. Coopmans, M. Benedetti, D. Muñoz Ramo and M. Rosenkranz, Sci. Rep., 2025, 15, 28409.
34. M. Motta, W. Kirby, I. Liepuoniute, K. J. Sung, J. Cohn, A. Mezzacapo, K. Klymko, N. Nguyen, N. Yoshioka and J. E. Rice, Electron. Struct., 2024, 6, 013001.
35. H. J. Vallury, M. A. Jones, C. D. Hill and L. C. L. Hollenberg, Quantum, 2020, 4, 373.
36. M. A. Jones, H. J. Vallury, C. D. Hill and L. C. L. Hollenberg, Sci. Rep., 2022, 12, 8985.
37. L. C. L. Hollenberg and N. S. Witte, Phys. Rev. D, 1994, 50, 3382–3386.
38. J. Cioslowski, Phys. Rev. Lett., 1987, 58, 83–85.
39. P. J. Knowles, Chem. Phys. Lett., 1987, 134, 512–518.
40. K. Kowalski and B. Peng, J. Chem. Phys., 2020, 153, 201102.
41. D. Claudino, B. Peng, N. P. Bauman, K. Kowalski and T. S. Humble, Quantum Sci. Technol., 2021, 6, 034012.
42. D. Horn and M. Weinstein, Phys. Rev. D, 1984, 30, 1256–1270.
43. Z. Ding and L. Lin, PRX Quantum, 2023, 4, 020331.
44. A. Asthana, A. Kumar, V. Abraham, H. Grimsley, Y. Zhang, L. Cincio, S. Tretiak, P. A. Dub, S. E. Economou, E. Barnes and N. J. Mayhall, Chem. Sci., 2023, 14, 2405–2418.
45. A. Kumar, A. Asthana, V. Abraham, T. D. Crawford, N. J. Mayhall, Y. Zhang, L. Cincio, S. Tretiak and P. A. Dub, J. Chem. Theory Comput., 2023, 19, 9136–9150.
46. M. D. Prasad, S. Pal and D. Mukherjee, Phys. Rev. A, 1985, 31, 1287–1298.
47. B. Datta, D. Mukhopadhyay and D. Mukherjee, Phys. Rev. A, 1993, 47, 3632–3648.
48. P. J. Ollitrault, A. Kandala, C.-F. Chen, P. K. Barkoutsos, A. Mezzacapo, M. Pistoia, S. Sheldon, S. Woerner, J. M. Gambetta and I. Tavernelli, Phys. Rev. Res., 2020, 2, 043140.
49. M. Motta, W. Kirby, I. Liepuoniute, K. J. Sung, J. Cohn, A. Mezzacapo, K. Klymko, N. Nguyen, N. Yoshioka and J. E. Rice, Electron. Struct., 2024, 6, 013001.
50. P. F. Kwao, S. P. Sundar, B. Gupt and A. Asthana, On the generalized eigenvalue problem in subspace-based excited state methods for quantum computers, 2025, arXiv:2503.09670,  DOI:10.48550/arXiv.2503.09670.
51. J. R. McClean, M. E. Kimchi-Schwartz, J. Carter and W. A. de Jong, Phys. Rev. A, 2017, 95, 042308.
52. J. R. McClean, Z. Jiang, N. C. Rubin, R. Babbush and H. Neven, Nat. Commun., 2020, 11, 636.
53. Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters and G. K.-L. Chan, WIREs Comput. Mol. Sci., 2018, 8, e1340.
54. K. Boguslawski and P. Tecmer, Int. J. Quantum Chem., 2015, 115, 1289–1295.
55. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa and P. van Mulbregt, SciPy 1.0 Contributors, Nat. Methods, 2020, 17, 261–272.
56. S. Sivarajah, S. Dilkes, A. Cowtan, W. Simmons, A. Edgington and R. Duncan, Quantum Sci. Technol., 2020, 6, 014003.
57. A. Javadi-Abhari, M. Treinish, K. Krsulich, C. J. Wood, J. Lishman, J. Gacon, S. Martiel, P. D. Nation, L. S. Bishop, A. W. Cross, B. R. Johnson and J. M. Gambetta, Quantum Computing with Qiskit, 2024, arXiv:2405.08810,  DOI:10.48550/arXiv.2405.08810.
58. M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information: 10th Anniversary Edition, Cambridge University Press, 2010.
59. J. Lee, W. J. Huggins, M. Head-Gordon and K. B. Whaley, J. Chem. Theory Comput., 2019, 15, 311–324.
60. D. I. Lyakh, M. Musiał, V. F. Lotrich and R. J. Bartlett, Chem. Rev., 2012, 112, 182–243.
61. P. Lolur, M. Skogh, W. Dobrautz, C. Warren, J. Biznárová, A. Osman, G. Tancredi, G. Wendin, J. Bylander and M. Rahm, J. Chem. Theory Comput., 2023, 19, 783–789.
62. H. Zou, E. Magnusson, H. Brunander, W. Dobrautz and M. Rahm, Multireference error mitigation for quantum computation of chemistry, 2025, arXiv:2505.08291,  DOI:10.48550/arXiv.2505.08291.
63. N. M. Tubman, C. Mejuto-Zaera, J. M. Epstein, D. Hait, D. S. Levine, W. Huggins, Z. Jiang, J. R. McClean, R. Babbush, M. Head-Gordon and K. B. Whaley, Postponing the orthogonality catastrophe: efficient state preparation for electronic structure simulations on quantum devices, 2018, arXiv:1809.05523,  DOI:10.48550/arXiv.1809.05523.
64. D. W. Berry, Y. Tong, T. Khattar, A. White, T. I. Kim, G. H. Low, S. Boixo, Z. Ding, L. Lin, S. Lee, G. K.-L. Chan, R. Babbush and N. C. Rubin, PRX Quantum, 2025, 6, 020327.

---

Footnote

† This is consistent with where . The sequence of GRs applied to the circuit returns the desired state vector through this recursive mapping of the normalized coefficients, as discussed in Arrazola et al.³

---