Bridging physical intuition and hardware efficiency for correlated electronic states: the local unitary cluster Jastrow ansatz for electronic structure

A prominent goal in quantum chemistry is to solve the molecular electronic structure problem for ground state energy with high accuracy. While classical quantum chemistry is a relatively mature field, the accurate and scalable prediction of strongly correlated states found, e.g., in bond breaking and polynuclear transition metal compounds remains an open problem. Within the context of a variational quantum eigensolver, we propose a new family of ansatzes which provides a more physically appropriate description of strongly correlated electrons than a unitary coupled cluster with single and double excitations (qUCCSD), with vastly reduced quantum resource requirements. Specifically, we present a set of local approximations to the unitary cluster Jastrow wavefunction motivated by Hubbard physics. As in the case of qUCCSD, exactly computing the energy scales factorially with system size on classical computers but polynomially on quantum devices. The local unitary cluster Jastrow ansatz removes the need for SWAP gates, can be tailored to arbitrary qubit topologies (e.g., square, hex, and heavy-hex), and is well-suited to take advantage of continuous sets of quantum gates recently realized on superconducting devices with tunable couplers. The proposed family of ansatzes demonstrates that hardware efficiency and physical transparency are not mutually exclusive; indeed, chemical and physical intuition regarding electron correlation can illuminate a useful path towards hardware-friendly quantum circuits.

In Table S1, we report additional calculations with the unitary coupled cluster ansatz (qUCCSD).In the case of ethene, unrestriction is necessary to approach the FCI energy E = −77.629025Hartree, although it produces a spincontaminated state with ⟨S 2 ⟩ ≃ 0.98.In the cases of cyclobutadiene and benzene, we observe no spin contamination, and good agreement between all the implementations of qUCCSD and with the FCI energies of E = −153.339314and E = −230.238284148Hartree, respectively.S1. qUCCSD calculations carried out in this study."flavor" refers to a restricted closed-shell (R) or unrestricted (U) implementation."formula" refers to whether the exponential of the cluster operator is approximated with a Trotter or Suzuki decomposition."ordering" refers to the ordering of operators (exponentials of doubles followed by singles, d+s, or viceversa, s+d) in the product formula approximation.
In Figure S1 we study the dissociation of ethene using the LUCJ ansatze with L = 1, 2, 3, 4. For L = 1, the accuracy of LUCJ is very limited, with all-to-all, square, and hex ansatzes overestimating the dissociation energy by 50 milliHartree, and heavy-hex by more than 100 milliHartree.The accuracy of LUCJ improves with increasing L. With the exception of heavy-hex, all ansatzes with L = 2 differ from FCI by less than 1.6 milliHartree.

S3
S3. ADDITIONAL LUCJ CALCULATIONS In Figures S2 and S3 we report additional LUCJ calculations for the H 2 and LiH molecules, in a basis of two Lowdin orbitals and three natural orbitals in the A 1 irrep of the C ∞v group, respectively.
-1.20 -1.10 -1.00 -0.90 -0.80 -0.70  FIG.S3.Potential energy curve of LiH in a (2e,2o) active space of natural orbitals, using the LUCJ ansatz with all-to-all, square-lattice, hex-lattice connectivity (top, middle, bottom), with or without J αα and X terms.All curves except those without J αα and X terms and hex-lattice connectivity (green lines, bottom panels) agree with FCI within 10 −6 Hartree.Because these unitaries conserve particle number and the Z-component of spin, we can restrict the simulation to the subspace spanned by basis vectors with the same particle number and spin as the initial state.The dimension of this subspace is D = N Nα N N β where N is the number of orbitals, N σ is the number of particles of spin σ, and n k is the binomial coefficient.The basis vectors are indexed by pairs of bitstrings (I α , I β ), where I σ is a bitstring of length N whose j-th bit is 1 if and only if orbital j of spin σ is occupied, and we use the notation j ∈ I σ to indicate this is the case.The coefficient corresponding to the basis vector (I α , I β ) is denoted γ Iα,I β , or simply γ I if we wish to ignore the spin.This wavefunction representation is the same as that used in classical full configuration-interaction (FCI) calculations, so existing FCI code can be repurposed.To this end, we used the open-source software package PySCF in our simulations [1].To simulate the LUCJ ansatz, we considered two methods.

A. Taylor series method
The exponential action of the operators (S1) can be approximated using a Taylor series expansion.This method requires only the ability to apply the operators (S1) to a wavefunction.Because similar operations are also performed in classical FCI calculations, it was relatively straightforward to implement these operations using PySCF.The application of a single term in a sum incurs a computational cost scaling as O(DN 2 ).To achieved a fixed precision, a constant number of terms in the Taylor series needs to be included.In practice, we found that fewer than 100 terms were needed to achieve a convergence threshold of 10 −12 as measured by the Euclidean norm of the residual vector.

B. Exact method via unitary decompositions
We also considered simulation methods inspired by quantum circuit implementations of the unitaries.These methods utilize a decomposition of the unitary into a sequence of basic operations and directly apply these operations by updating the coefficients stored in the state vector.While these methods have the same asymptotic complexity as the Taylor series method, the constant factors are much smaller, and we found that in practice they were on the order of a hundred times faster.In the following sections, we describe how to implement each of the two kinds of unitaries.

Orbital rotation
The unitary e Kµ represents a rotation of the orbital basis.It can be decomposed as where each P j is a phase operator of the form and each G j is a Givens rotation operator acting on neighboring fermionic modes [2,3].A Givens rotation operator G acting on modes p and q has the action where G is a 2 × 2 unitary matrix of the form The number of Givens rotations is N G = N (N − 1)/2.The Givens rotation matrices and phase angles can be calculated from the matrix e K µ using the procedure described in Ref. [3].Once these quantities are calculated, the orbital rotation can be achieved by applying the sequence of operations given by Eq. (S2) directly on the state vector.
In the following, we will describe how to apply the operations considering a single spin.The full operation is achieved by applying the single-spin operation for both spins in sequence.Let γ I denote a coefficient of the original state vector and γ ′ I a coefficient of the transformed state vector.Applying a phase operator P j amounts to multiplying by a phase the coefficients corresponding to basis states in which the corresponding orbital is occupied: Applying a Givens rotation G acting on neighboring modes p and q is achieved by rotating the subspaces spanned by basis vectors in which exactly one of these modes is occupied: cγ I + sγ Ip →q if p ∈ I and q / ∈ I −s * γ Iq →p + cγ I if q ∈ I and p / ∈ I γ I otherwise. (S7) Here, the notation I p →q denotes the bitstring obtained from I by setting bit p to 0 and bit q to 1. Since there are O(N 2 ) total operations to apply, the compute time for performing a full orbital rotation is O(DN 2 ).

Diagonal Coulomb operator
The unitary e i Ĵµ is equal to a product of commuting terms where each term has the form exp(iφn pσ nqτ ) for some angle φ.This term is a diagonal operator with the following action: