Spin-projection for quantum computation: A low-depth approach to strong correlation

Although spin is a core property in fermionic systems, its symmetry can be easily violated in a variational simulation, especially when strong correlation plays a vital role therein. In this study, we will demonstrate that the broken spin-symmetry can be restored exactly in a quantum computer, with little overhead in circuits, while delivering additional strong correlation energy with the desired spin quantum number. The proposed scheme permits drastic reduction of a potentially large number of measurements required to ensure spin-symmetry by employing a superposition of only a few rotated quantum states. Our implementation is universal, simple, and, most importantly, straightforwardly applicable to any ansatz proposed to date.


I. INTRODUCTION
Recent advancements on quantum devices have created widespread interest in the development of efficient quantum algorithms.The Variational Quantum Eigensolver (VQE) is one of the most pursued approaches in the NISQ (noisy intermediate-scale quantum) era, where a trial state |ψ(θ) parameterized by θ is variationally optimized to minimize its Hamiltonian expectation value.[1,2] Among the several candidates for |ψ , the unitary coupled cluster with singles and doubles (UCCSD) [3][4][5] ansatz has been extensively used as an entangler in the preparation of a trial state from Hartree-Fock (HF) |Φ , as shown by the following equation: Here, . are real parameters and anti-hermitian pairs of the kth excitation and de-excitation operators in a spin-orbital basis.In this work, we use convention for the orbital indices: i, j for the occupied orbitals of |Φ HF , a, b for the virtual orbitals, and p, q, r, s for the general orbitals.To make the unitary exponential operator programmable on a quantum device, Trotterization is required in practice for non-commutative exponents, i.e., e T1+ T2 ≈ ai e t a i τ a i /µ abij e t ab ij τ ab ij /µ µ .[6][7][8] However, applying UCCSD to strongly correlated systems often triggers large t-amplitudes to account for higher excitations, which, in turn, necessitates a large Trotter number µ, thereby requiring a high-depth quantum circuit.Since Trotterization is nothing but an artifact needed to respect the ansatz Eq. ( 1), many authors have studied the efficacy of a single Trotter step (µ = 1).[6,9,10] More flexible ansätze have been also proposed, where, * tsuchimochi@gmail.com in the product of exponential operators, the same antihermitian excitations may be repeated albeit with different amplitudes.[11,12] To our knowledge, these recent studies on UCC and its variants have vastly neglected the spin properties of the obtained solutions.Thus, their discussions remained somewhat ambiguous.We argue that it is extremely important to monitor Ŝ2 because, even with a spin-restricted HF (RHF) reference, UCC amplitudes can spontaneously violate spin-symmetry, thereby variationally lowering its energy as opposed to traditional CC.This often happens especially in strongly correlated systems, such as bond dissociations, as will be demonstrated below.Such broken-symmetry (BS) solutions are not physical for non-relativistic Hamiltonians.Moreover, the problems caused by spin-contamination are well known and documented.[13][14][15] This so-called symmetry dilemma poses a challenge in quantum computation.As will be shown, except for singles, using spin-free generators cannot fix this problem because each term is not necessarily commutative.Thus, such "spin-adapted" (SA) methods [16,17] still incur Trotterization for large t ab ij .Some of the previous studies on the matter have suggested the use of a constrained approach, [14,[18][19][20] wherein the Hamiltonian is augmented with a penalty term λ Ŝ2 − s(s + 1) 2 , where λ → ∞ allows the elimination of spin-contamination.However, whereas Ŝ2 can be evaluated with the same effort as the energy by measuring O(n 4 ) terms in the operator, the measurement of the expectation value of Ŝ4 has a steep O(n 8 ) scaling, where n is the number of qubits (i.e., spin orbitals), compared to the O(n 4 ) scaling of the bare Hamiltonian.Ryabinkin et al. have attempted to avoid this prohibitive scaling by constraining the average spin Ŝ2 .[20] Nevertheless, the constrained approach does not seem particularly suited to spin as the Trotterized UCCSD ansatz is not flexible enough, provoking the symmetry dilemma.The resulting energy can become unreliable depending on the choice of λ, if the ansatz is not flexible enough to deal with strong correlation.
In this work, we propose an alternative way of preserving the spin of an arbitrary quantum state in a numerically exact, but much less costly, manner, while simultaneously treating strong correlation.To this end, we will let a quantum state break its spin symmetry, which is subsequently recovered by means of symmetry projection.By writing a spin-adapted state as a superposition of rotated broken-symmetry states, the quantum circuit we develop here fully appreciates the advantage of VQE by replacing a potentially long quantum circuit to describe entanglements with many measurements with short circuits.
We organize the present work as follows.In Section II A, we introduce spin-adapted UCCSD for quantum computing, and in Section II B, we demonstrate the practical difficulty in the spin-constrained approach to handle the Lagrange multiplier λ.Section II C then introduces the spin-projection technique for VQE as well as the required quantum circuit.Section III provides some illustrative calculations that are challenging for standard UCC.Finally, we draw our conclusions in Section IV.

II. THEORY A. Spin-adapted UCCSD
To circumvent the variational collapse via spinsymmetry breaking in UCCSD, one can introduce spinadapted (SA)-UCCSD, where the T1 and T2 operators are constructed with unitary group generators Êa i = a † a a i + a † āaī, with a bar indicating the β spin (α without a bar).That is, we use the same amplitude for spincomplement operators, or, equivalently, introduce the following restrictions,[17] We note that SA-UCCSD uses exactly the same circuit as BS-UCCSD but with less number of amplitudes.This inexactness also holds for generalized doubles, which uses general orbitals.Hence, we assume that, for example, the adaptive approach of Grimsley et al. [12] can suffer from spin-contamination (incidentally, we should also point out that the operator set employed in Ref. [12] does not constitute spin-free operators anyway, treating {τ pq rs , τ p q rs } and {τ pq rs , τ pq rs } differently as opposed to the above equation).The dilemma is that the extent of spin-contamination varies by size of amplitudes; so, for strongly correlated systems where t amplitudes can be significantly large to describe higher excitation effects, one should pay attention to Ŝ2 and Trotterization may be eventually required to correctly obtain the desired spin state.

B. Spin-constrained UCCSD
A constraint on the spin s can be applied to (Trotterized) UCCSD by augmenting the energy functional with a penalty term, becomes less spincontaminated.However, when |ψ UCCSD is Trotterized, this becomes a very challenging task.Figure 1 shows the total energy of Trotterized UCCSD (µ = 1) and its spin expectation value Ŝ2 of the singlet nitrogen molecule at a bond length of 2.8 Å by varying λ.We have used the STO-6G basis set with 1s and 2s orbitals frozen.While λ is small, the energy is almost unchanged but so is Ŝ2 , suffering from spin-contamination.However, requiring the latter to be sufficiently precise with a large λ introduces an enormous increase in the evaluated energy, which can be on the order of tens of mHartree.With such a significant energy change, it is extremely difficult to determine the appropriate λ.This test case clearly demonstrates the Trotterized UCCSD ansatz is not flexible enough to capture both the correct spin and strong correlation.To eliminate the strong dependency of energy on λ, one is required to employ large µ, with which we presume Eq. 4 would at least give the result of a similar quality to SA-UCCSD.

VQE algorithm
The spin-projection operator P(s) has proven itself useful not only in removing spin contamination but also in capturing strong correlations with its inherent high excitation effects when used with BS ansätze.[21][22][23] This is manifested from the well-known form due to Löwdin, which is written as a product of Ŝ2 operators.[24] From this perspective, one can expect that the role of exponential doubles e t ab ij τ ab ij in the spin-projected UCCSD, | ψUCCSD = P|ψ UCCSD , should be to mainly describe weak dynamical correlation, with moderately small amplitudes t ab ij .Hence, the introduction of P can lead to tremendous error reduction in the Trotter approximation in UCCSD.However, despite its formal simplicity and appealing properties, the many-body nature of Löwdin's P makes it intractable, even with a quantum computer, due to the need of evaluating Ĥ Ŝ2 , Ĥ Ŝ4 , • • • , which results in the exponential increase in the number of Pauli operators, [19] hindering its practical usage.
Hence, we seek the possibility of the different representation of P, which is particularly suitable for quantum computation.It is known that a projector onto the Hilbert space of spin s and spin magnetization m is written in the following integral form, [25][26][27] where Ω = (α, β, γ) are the Euler angles, D s mm (Ω) = s; m| R(Ω)|s; m is the Wigner D-matrix, and Û (Ω) is known as the spin-rotation operators, whose resemblance to the Pauli-rotation operations is striking.Since each of the exponents of Û(Ω) is an anti-hermitian one-body operator, the many-body effect of the Löwdin operator has been translated to a (infinite) linear combination of orbital rotations.Stated differently, a spin-projected state | ψ(θ) = P|ψ(θ) with |ψ(θ) being a broken-symmetry state, is understood as a superposition of broken-symmetry nonorthogonal states Û(Ω)|ψ(θ) .For a spin-free observable Ô, its expectation value is simplified as P † Ô P ≡ Ô P .Because of the non-unitary operation of P, |ψ does not preserve its norm upon projection.Therefore, the energy expectation value for projected VQE becomes It is noteworthy that |ψ can be any quantum state.If |ψ is an eigenfunction of Ŝz with an eigenvalue of m such as UCC, one may a priori integrate out α and γ and use the operator P instead of P, where we have introduced the Gauss-Legendre quadrature using the N g sampling points {0 ≤ β g ≤ π : g = 1, • • • , N g } for rotations Ûg = e −iβg Ŝy along with accordingly defined weights w g .With this form, one typically requires only a few grid points to achieve Ŝ2 = s(s + 1) that is accurate to several decimal points, and the convergence is exponentially fast.[26,28] To be precise, in this work, we propose to evaluate Ô Ûg θ for UCC with a quantum computer and then use a classical adder to evaluate E, instead of constructing a symmetry-projected ansatz.
With a classical computer, the cost of evaluating Ô Ûg θ scales exponentially if |ψ(θ) is an exponential ansatz, such as coupled-cluster, [29][30][31][32] since the orbital rotations Ûg are still many-body operators that include de-excitations.On the other hand, with a quantum computer, it is reduced to a polynomial cost by the Hadamard test with a Ûg gate controlled by an ancilla qubit |q anc (see Appendix A). [33][34][35] 2. Quantum Circuit Now, we consider an efficient circuit that is optimal for the controlled-U g gate.We start by writing Ŝy in the second quantization as a linear combination of imaginary spin-flip excitations, where S pq = φ p |φ q is the overlap between the spatial orbitals of α and β spins.It is important to note that although almost all spin-projection methods are based on a broken-symmetry basis of spin-unrestricted HF (UHF), their orbital set is spatially nonorthogonal, i.e. −1 ≤ S pq ≤ 1. [15,21] While the Trotterization can be circumvented for singles, [38] the implementation with UHF can add a non-negligible amount of CNOT gates, which result in an overhead to a quantum circuit.We wish to minimize it.In contrast, a quantum circuit for the controlled-U g gate can be made substantially simpler in an RHF orbital basis, since S pq ≡ δ pq .This gives rise to no Trotter error since each component is commutative with one another.Therefore, It can be shown that not only is the number of gates reduced, but also that the controlled-U g gate with Eq. ( 10) can be efficiently implemented with the Jordan-Wigner transformation.Suppose spin-orbitals are mapped onto qubits with spins alternating as Then, we have, in an RHF orbital basis, which does not entail a sequence of nearest-neighbor CNOT gates that would be needed for a product of σ z between p and q to ensure the anti-commutation relation of Fermion operators.[6,8] Thus, the controlled-U g gate has a nice block structure with local (controlled) spinflips e βg 2 τ p p (Figure 2).A quantum circuit that performs a local spin-flip is given in Figure 3, along with a controlled-R z gate decomposed to two R z and CNOT gates to facilitate actual implementations.With this circuit, the overall overhead introduced by the controlled-U g gate is negligible, with the number of additional CNOT gates being only 4n.Therefore, a large number of measurements required for the purposes of preserving spin-symmetry in both the constrained approach and Löwdin's spinprojection scheme can be completely avoided with our scheme, which needs only O(N g n 4 ) measurements.
Implementing e −iα Ŝz is straightforward independent of the basis, thanks to the diagonal nature of Ŝz .A full quantum circuit for Û(Ω) is also available in Appendix B.

Broken-Symmetry Ansatz
While the foregoing discussion holds true with regard to the use of an RHF orbital basis, spin-projection can be even more beneficial if combined with physicallymotivated broken-symmetry ansätze.[39,40] One could use BS-UCCSD for this purpose, but here we describe a general procedure to deliberately break spin-symmetry by applying spin-dependent orbital rotations to a quantum state |ψ in a way similar to Ûg .To achieve this task, we conveniently write the orbital rotation operator K defined as a product of spin-dependent Givens rotations {e κ p q τ p q , e κ p q τ p q }.This procedure has been used to generate arbitrary Slater determinants.[38,41] Applying K to RHF generates any UHF state represented by the RHF orbitals.This means that the projected HF (PHF) can be simply given by It is worth noting that the orbital basis is still that of RHF, whereas UHF is expressed by a linear combination of excited determinants from |Φ RHF -the Thouless theorem.[42] We can extend this scheme to any other VQE ansatz including UCC without loss of generality.[43][44][45][46] For UCCSD, however, T1 and K play similar roles and are considered largely redundant.Therefore, in our study we omit T1 and consider UCCD combined with K and spin-projection: where both the κ and t amplitudes are fully optimized.We call this scheme projected UCCD (PUCCD).Note that the orbital-optimization effect is encoded in K.This way, the amplitudes in T2 are expected to be small, and higher excitation effects needed for strong correlation are now shifted to the orbital rotation followed by spinprojection.Accordingly, in many cases, e T2 in Eq. ( 13) is well approximated by µ = 1, which is referred to as the disentangled PUCCD (dPUCCD), following the work of Evangelista and co-workers.[9] HF and UCC are both invariant with respect to orbital rotations within occupied and virtual spaces.[47] Applying spin-projection does not change this property.[31] Therefore, it is enough to take into account only the occupied-virtual orbital rotations with {κ a i , κ ā ī }.However, introducing the Trotter approximation in e T2 of Eq. ( 13) no longer guarantees the said invariance.Nevertheless, our experiences indicate that occupied-occupied and virtual-virtual rotations are very much redundant even for dPUCCD, considering that their energy derivatives were found to be numerically zero.Thus, they will not be taken into consideration below.

III. ILLUSTRATIVE CALCULATIONS
We have used PySCF [48] to generate molecular integrals, OpenFermion [49] to perform Jordan-Wigner transformation of Ĥ and Ŝ2 , and Qulacs [50] to construct quantum circuits for simulation on a classical computer.For optimization, the L-BFGS method was used.[51]N g required for the numerically exact projection varies depending on the degree of spin-contamination; we have used N g = 2 for N 2 and H 2 O and N g = 3 for O, which are found to be sufficient for the Ŝ2 values to be accurate to at least 10 −10 .The exponential operators in Trotterization are applied in the following order: abij → āb īj → ābīj → ai → āī , where, in each spinblock, the leftmost virtual label is the outermost loop while the rightmost occupied label is the innermost loop.
The Trotter approximation is considered as having been converged with µ = 10.
We first assess the improvement that spin-projection has to offer by computing the potential energy curves of the triple bond dissociation of the singlet N 2 (m = 0).We have used the STO-6G basis set and frozen 1s and 2s core electrons of atoms, resulting in an active space of (6e, 6o).Plotted in Figure 4 are the energy errors from FCI (kcal/mol) for N 2 in a logarithmic scale.The shaded area corresponds to the "chemical accuracy," i.e., errors less than 1 kcal/mol.In the lower panel, Ŝ2 of UHF and BS-UCCSD are likewise depicted.All other methods numerically preserve the correct singlet value Ŝ2 = 0 and are, therefore, not shown.The energy of BS-UCCSD starts deviating from that of SA-UCCSD by breaking the symmetry restriction in amplitudes.It is found that this so-called Coulson-Fischer point happens at a much longer distance for UCCSD (2.0 Å) compared to that of HF (1.2 Å).This is because strong correlation can be partially captured via exponential excitations in UCCSD, thereby mitigating spin contamination.However, the ansatz only accounts for disconnected higher excitations, e.g., T 2 2 , resulting in larger errors.On the other hand, it is worth noting that dPUCCD delivers essentially exact results for these systems while guaranteeing the singlet spin; the largest deviation is 0.007 kcal/mol at 1.5 Å.The error is at least one order of magnitude smaller than those of SA-UCCSD and BS-UCCSD.
Both in PHF and PUCCD, higher excitations are treated via spin-projection and orbital-rotations K,  which do not need Trotterization.Hence, it is expected that the doubles amplitudes of PUCCD are small, thereby allowing for a small Trotter error in dPUCCD compared to those of the disentangled UCCSD (dUCCSD).This is demonstrated in Figure 5, where we have plotted the Trotter error ∆ of each disentangled scheme in the dissociation of N 2 .In Table I, we tabulated the norms of singles (either |t| or |κ|) and doubles amplitudes obtained at R N−N = 2.2 Å with µ = 10 as well as the distances from those of µ = 1, as another indicator of the Trotter error.As can be clearly seen, PUCCD gives the minimal T 2 amplitudes, thereby allowing for an efficient treatment of strong correlation.In contrast, the converged T 1 -amplitudes of BS-UCCSD are considerable and require a re-optimization of amplitudes for different µ.This indicates that BS-dUCCSD is a whole different ansatz from BS-UCCSD, even for a large bond length, R N−N = 2.8 Å, where both yield almost identical energies, as is evident from the totally different Ŝ2 , i.e., 1.37 for µ = 1 and 2.55 for µ = 10.The above test system employs the minimal basis and admittedly neglects most of the dynamical correlation effect.To show that the accuracy of PUCCD is not a fortuitous artifact resulting from the chosen basis, we have also carried out the calculations for the symmetric bond dissociation of H 2 O with a fixed angle of 104.5 • using 6-31G and STO-6G for oxygen and hydrogen, respectively.Table II summarizes the energy errors of each disentangled method from FCI, where the numbers in parentheses imply those obtained with the STO-6G basis for all the atoms for comparison.The presence of non-negligible dy-namical correlation leads to substantial errors if strong correlation is not treated appropriately.This is the case for SA-and BS-UCCSD, while PUCCD still remains to achieve the chemical accuracy.
Lastly, we computed the singlet-triplet energy gap of the oxygen atom to show that neglecting spin-symmetry can often result in a catastrophic error.The basis set used is 6-31G with 1s frozen core, constituting another example with a balanced mixture of dynamical and strong correlation effects.Table III lists the singlet and triplet energies in Hartree, computed with m = 0 and 2, respectively, along with their difference ∆ ST in kcal/mol.It is readily clear that both SA-and BS-UCCSD schemes failed to describe strong correlations of the open-shell singlet oxygen.BS-UCCSD with m = 0 spontaneously breaks spin and ends up in mainly the low-spin triplet state ( Ŝ2 = 1.96).Spin-projection in dPUCCD enables to capture strong correlation in the singlet oxygen.As a result, its ∆ ST is found to be in excellent agreement with FCI.

IV. CONCLUSIONS
In this Letter, we have shown that numerically exact spin-projection can be made feasible within the framework of VQE.It should be emphasized that spinprojection attains remarkable accuracy in exchange for one ancilla qubit, albeit with little overhead in circuit depth (additional 4n CNOT gates).We note that one disadvantage of our scheme is the increase in the number of measurements, but it grows only linearly with N g as opposed to standard nonorthogonal methods that show a quadratic scaling.[37] Also, the scheme can potentially suffer from errors when a projected state is not suitable for dealing with strong correlation in the system [52,53].Nonetheless, we should once again stress that the spinprojection operator lends itself straightforwardly to any ansatz in order to guarantee the desired spin-symmetry, and therefore is a promising, versatile tool.
For the full-spin projection for a spin-generalized quantum state (both Ŝ2 and Ŝz symmetries are broken), one has to implement Û (Ω) (Eq.( 6)) instead of the simplified rotator Ûg .This necessitates a quantum circuit that performs e −iθ Ŝz for θ = α, γ, where This means that e −iθ Ŝz only requires single qubit operations, where all qubits that represent up-spin and downspin orbitals are rotated by R z − θ 2 and R z θ 2 , respectively.Therefore, a quantum circuit that perform the full spin-projection with Û (α, β, γ) simply becomes as follows:

FIG. 1 .
FIG. 1. Spin-constrained UCCSD energy and its spin expectation value Ŝ2 of the singlet nitrogen molecule at a fixed bond length of 2.8 Å for different λ.The dashed lines indicate the values computed with λ = 0.

FIG. 2 .
FIG.2.The block structure of the controlled-Ug circuit in spin-restricted orbital basis.

FIG. 3 .
FIG. 3. Quantum circuits for (a) the spin-flip operation by angle βg and (b) controlled-Rz.H and Y = Rx(−π/2) respectively transform a qubit to the σ x and σ y basis.

2 pp ap e i θ 2 a † p a p = p e −i θ 4 (Ip−σ z p ) e i θ 4 (
of Ŝz makes it particularly easy to implement its exponential form on a quantum computer.Namely,e −iθ Ŝz = e −i θ Ip−σ z p )

TABLE I .
Norms of singles and doubles amplitudes for N2 at 2.