Quantum simulations employing connected moments expansions

Quantum computing (QC) techniques attract much attention in many areas of mathematics, physics, and chemistry by providing a means to address insurmountable computational barriers for simulating quantum systems on classical computers. One of the best illustrations of this fact is associated with solving Schrödinger equations for many-electron systems in quantum chemistry, where the attempt of including all configurations spanning Hilbert space (or at least configurations from the subspace relevant to a problem of interest) quickly evolves into numerical problems characterized by exponentially growing numerical cost. Unfortunately, for a large class of the problems commonly referred to as the strongly correlated systems, the inclusion of all configurations is necessary to obtain a desired level of accuracy. To bypass these problems, several classes of many-body methods targeting the so-called full configuration interaction (FCI) accuracy have been developed including high-rank coupled-cluster (CC),^1–3 selected CI,^4,5 density matrix renormalization group (DMRG),^6–9 stochastic and semi-stochastic CI,^10–12 hybrid stochastic CI and CC,^13–15 and virtual orbitals many-body expansions¹⁶ methodologies. Recently, accuracies of these formalisms have been evaluated using a benzene system,¹⁷ where considerable differences in the calculated energies with high-accuracy formalisms were reported.

Although QC has not reached its maturity yet, intensive effort has been directed toward developing algorithms for quantum simulations that can take advantage of early noisy quantum registers. Among several algorithms that are being developed for this purpose, one should mention the hybrid quantum/classical Variational Quantum Eigensolver (VQE) approach,^18–25 which has been intensively tested from the point of view of its utilization on noisy intermediate-scale quantum devices (NISQ). Other quantum algorithms such as quantum phase estimation (QPE),^26–33 and more recently quantum algorithms for imaginary time evolution (ITE),^34,35 quantum filter diagonalization,³⁶ quantum inverse iteration algorithms,³⁷ and quantum power methods³⁸ seem to be more hardware demanding than the VQE algorithm. Some of the above mentioned algorithms may draw heavily on the utilization of Trotter formulas for evaluating the exponential forms of the operators corresponding to unitary coupled cluster Ansatz, unitary time evolution of the system, and the solution of the imaginary-time Schrödinger equation, respectively. For example, for the unitary time evolution defined by qubit-encoded Hamiltonian H (various types of encoding algorithms for many-body Hamiltonians—including Jordan–Wigner³⁹ and Bravyi–Kitaev⁴⁰ methods—are discussed in Refs. 30 and 41), where H is the sum of multiple (generally non-commuting) terms H_j, $H=∑j=1MHj$ (here, H_j = h_jP_j with h_j being complex numbers and P_j representing tensor products of Pauli matrices and/or identity matrices), it is hard to encode e^−iHt directly into quantum gates. Instead, in quantum simulations, one resorts to the Trotter formula^42,43 with K Trotter steps,

where Δ_t = t/K, and for one Trotter step, U(Δ_t) is further approximated as $∏j=1Me−ihjPjΔt$⁠. Note that even in its simplest case (K = 1), Eq. (1) results in arbitrary connectivity between qubits. With the number of terms in Hamiltonian being proportional to $O(N4)$ (where N stands for the number of one-particle basis functions employed), the optimization of the corresponding circuits plays a crucial role in the quantum simulations (especially when employing the techniques related to the fermionic swap network⁴⁴). Nevertheless, even for formulations using a reduced number of terms [e.g., proportional to $O(N2)$] in many-body Hamiltonian, the resulting gate depth still depends polynomially on N and the total number of CNOT gates is still proportional to $O(N)$ (letting alone the multiple CNOT gates for multi-qubit operators needed, for example, in QPE).

In this Communication, we propose an alternative approach for quantum simulations of many-body systems, where the number of instances for using the Trotter formula is significantly reduced. For this purpose, we utilize connected moments expansion introduced by Horn and Weinstein⁴⁵ and further advanced by Cioslowski⁴⁶ and others,^47–54 where the energy of a quantum system can be calculated using trial wave function |Φ⟩ (non-orthogonal to the exact ground state |Ψ⟩) and the expectation values of the Hamiltonian powers ⟨Φ|Hⁿ|Φ⟩. Here, ⟨Φ|Hⁿ|Φ⟩’s are calculated using simple variants of the Hadamard test, where CNOT gates are used to control the products of the unitaries P_j. Additional CNOT gates are optionally required for strongly correlated systems where correlation effects need to be captured in the process of trial state preparation, which in a realistic setting requires only the inclusion of a relatively small number of parameters in the VQE formulations. We also evaluate the effect of noise on the accuracy of quantum connected moment expansion (CMX) simulations.

The connected moments expansion is derived from the Horn–Weinstein (HW) theorem⁴⁵ that relates the function E(τ) to a series expansion in τ with coefficients being represented by connected moments I_k,

Here, the connected moments I_k are defined through a recursive formula (see Ref. 45 for details),

It can be seen that E(τ) is monotonically decreasing and lim_τ→∞E(τ) corresponds to the exact ground-state energy E₀, which is a consequence of the fact that e^−τH contracts any trial wave function |Φ⟩ toward the true ground state |Ψ⟩ (assuming ⟨Φ|Ψ⟩ ≠ 0). In practical applications based on the truncated moments expansion, Padé approximants need to be used to reproduce the proper behavior in the τ = ∞ limit. The behavior of this expansion and some techniques for maintaining size-extensivity are discussed in Ref. 48 (see also Ref. 55).

By applying re-summation techniques to Eq. (2), Cioslowski⁴⁶ derived an analytical form for exact energy in the τ = ∞ limit,

where S_k,1 = I_k (k = 2, 3, …) and $Sk,i+1=Sk,1Sk+2,i−Sk+1,i2$⁠. The truncation of the CMX series (4) after the first K terms leads to the CMX(K) approximations. For example, CMX(2) and CMX(3) energies are defined as follows:

The CMX formalism provides a trade-off between the rank of the connected moments included in the approximation and the quality of the trial wave function. Specifically, it allows one to use the multi-configurational⁴⁶ or even correlated representations (in the form of truncated CI or CC wave functions⁵⁶) of the trail wave functions for quasi-degenerate sought-for-states.

The analytical properties of CMX expansions and their relations to Lanczos methods have been discussed in the literature.^{47,51,52,57–60} Several techniques have been introduced to counteract possible problems associated with singular behavior of the CMX expansions (as pointed out by Mancini, there exists an infinite number of valid CMX expansions⁵²) including Knowles’s generalized Padé approximation,⁴⁷ alternate moments expansion (AMX),⁵² and generalized moments expansion.⁵⁹ Additionally, as suggested by Noga et al. in Ref. 56, a proper definition of the trial wave function may significantly alleviate possible issues associated with singular behavior. An interesting yet less known CMX formulation has been proposed by Peeters and Devreese⁶¹ (further analyzed by Soldatov⁶²) where the upper bound of the ground-state energy in n-th order approximation is calculated as a root of the polynomial P_n(x),

where a₀ = 1 and coefficients a_i’s for 1 ≤ i ≤ n (forming vector a) are obtained by solving linear equations

with matrix elements M_ij = ⟨Φ|H^2n−(i+j)|Φ⟩ and vector component b_i = ⟨Φ|H²ⁿ⁻ⁱ|Φ⟩. Remaining roots of the polynomial correspond to the upper bounds of the excited-state energies. While both CMX and PDS (PDS refers to the formulation of Peeters, Devreese, and Soldatov) methods are invariant under orbital rotations, the PDS approach, in contrast to the CMX expansions, does not furnish size-extensive energies.

Recently, QC algorithms for ITE and quantum Lanczos algorithm^34,35,63 have attracted a lot of attention. One of the most appealing features of these methods is the avoidance of a large number of ancillae qubits and complex circuits. In this context, we believe that the CMX framework provides another way for the utilization of near-term quantum architectures, where one can almost entirely eliminate Trotter steps and significantly reduce the number of CNOT gates. Below, we will describe a simple quantum algorithm for calculating moments K_n = ⟨Φ|Hⁿ|Φ⟩ for an arbitrary trial state.

In this section, we describe a general structure of the quantum algorithm for calculating moments K_n. Further simplifications and optimization of the circuit can be achieved by adding more ancillae qubits and/or merging with other quantum algorithms.

The main difference between our algorithm and ITE algorithms (for example, the QITE approach of Ref. 35) is the fact that, in the latter approach, the infinitesimal ITE generated by unitary H_j, i.e., $e−ΔτHj$⁠, is mirrored by a unitary evolution generated by the $e−iΔτAj$ operator acting onto properly normalized states, while in the former method, the standard moments K_n = ⟨Φ|Hⁿ|Φ⟩ are directly calculated. Several quantum algorithms have been proposed for computing K_n. For example, since $Hn=in∂n∂tneiHt|t=0$⁠, Seki and Yunoki³⁸ showed that K_n can be generally approximated by a linear combination of the time-evolution operators at n + 1 different time variables. Moreover, to estimate ground state properties of a programmable quantum device, Kyriienko^37,64 proposed a quantum inversion iteration algorithm where the inverse of Hamiltonian powers is decomposed as a sum of unitary operators. Here, in this work, different from the above-mentioned approaches, we directly employ the Hadamard test to evaluate the contributions to K_n for constructing connected moments and computing CMX energies.

The Quantum CMX (QCMX) algorithm works as follows: (1) Preparing an initial state |Φ⟩. In our numerical tests, we used a simple single Slater determinant for |Φ⟩. Alternatively, one can envision the use of a “non-aggressive” variant of VQE where a small number of amplitudes that define basic static correlation effects in the sought-for-wave-function are optimized. (2) Performing the Hadamard test to evaluate p_J = ⟨Φ|U_J|Φ⟩, where $UJ=∏k=1lPjk$ [for CMX(K) formulation, l < 2K − 1] (see Fig. 1) and cumulative index J designates l-tuple {j₁, …, j_l}. The contribution from p_J to K_l is equal to $pJ×∏k=1lhjk$⁠. (3) Compute connected moments I’s from moments K’s. (4) Once all I_l (l = 1, …, 2K − 1) are known, we choose the optimal form of the CMX expansion. For example, if I₃ is close to 0, then we choose a CMX form that is not using I₃ in the denominator (this will provide an optimal utilization of the information obtained from QC). Additionally, using identity σ_iσ_j = δ_ijI + iϵ_ijkσ_k (i, j, k = x, y, z), one can reduce entire U_J to an effective unitary $P¯J$ corresponding to a tensor product of Pauli matrices and/or identity matrices with an appropriate product of phase factors (Fig. 2). It means that the gate depth is exactly the same irrespective of the rank l of the calculated moment K_l (of course, the number of terms to be evaluated in this way increases with the increase in rank). Summarizing, the gate depth of the quantum CMX algorithm is mainly determined by the wave function preparation step and the practical realization of multi-qubit CNOT in the Hadamard test. In the case where the state preparation corresponds to flipping states from |0⟩ to |1⟩ using X_i operators (for quantum registers corresponding to occupied spin-orbitals), the resulting gate depth is extremely shallow. The utilization of VQE for |Φ⟩ state preparation increases the gate depth. Although if it includes only a small number of amplitudes, the resulting gate depth should not be excessive. Moreover, employing a trial wave function encapsulating basic correlation effects may have positive effects on the accuracy of finite-order CMX formulations. When an entangled trial wave function is used, one can also see advantages of using quantum computers—complexity and numerical overhead of the expressions defining expectation values of Hⁿ operators for correlated trial wave function(s) grow rapidly making classical calculation infeasible.

A large number of measurements is a common bottleneck of quantum algorithms, including VQE, ITE, and QCMX formalisms. Due to its simplicity, for the QCMX formalism, a simple reduced-measurement variant can be proposed. For example, one can select only the leading terms in the Hamiltonians’ second-quantized form and perform the fermion-to-qubit mapping, which results in the Hamiltonian H′ characterized by a much smaller number of terms, i.e.,

where M_R < M. One should also notice that a single measurement of specific matrix element

(where for generality, we assume that the trial wave function can be optimized with respect to some variational parameters θ) can be used to determine classes of contributions to moments of arbitrary rank. For example, for “diagonal” products, one has

where we used the fact that ⟨Φ(θ)|P(R)_jP(R)_j|Φ(θ)⟩ = 1 for normalized |Φ(θ)⟩. In a similar way, one can calculate contributions from ⟨Φ(θ)|P(R)_iP(R)_j|Φ(θ)⟩, i.e., ⟨Φ(θ)|P(R)_iP(R)_jP(R)_iP(R)_j|Φ(θ)⟩, ⟨Φ(θ)|P(R)_iP(R)_jP(R)_iP(R)_jP(R)_iP(R)_j|Φ(θ)⟩, etc. When combined, these two techniques can significantly reduce the number of measurements and enable a simple evaluation of higher-order moments (see the supplementary material for more discussion). In the quantum CMX extensions, in order to reduce the effort associated with calculating ⟨Φ|Hⁿ|Φ⟩, we are also planning to utilize quantum power methods introduced in Ref. 38.

For building quantum circuits to calculate the moments, K_l’s, we used Qiskit software.⁶⁵ As benchmark systems, we chose H₂ molecules in minimum basis (for various geometries corresponding to situations characterized by weak and strong correlation effects)⁶⁶ and a two-site single-impurity Anderson model (SIAM) (for a broad range of hybridization strength V) described by model Hamiltonians. In our studies, we used various orders of Cioslowski CMX(K), Knowles CMX(K),⁴⁷ and K-th order Peeters, Devreese, Soldatov expansion, PDS(K). In our quantum simulations, individual spin orbitals are directly assigned to qubits. For each qubit, |1⟩ and |0⟩ correspond to occupied and unoccupied states of the spin orbital. We applied simple X gates on the vacuum state to generate our trial states. Fidelities of trial states are given in the supplementary material. The results of quantum simulations are shown in Figs. 3 and 4.

In Fig. 3, all three CMX variants with up to fourth order expansions are able to reproduce the FCI energy in at least the mHartree level for R_H–H being up to 1.50 Å. The only exception is PDS(1), which corresponds to I₁ = ⟨Φ|H|Φ⟩. The higher-order PDS expansions perform exceptionally well, and the deviations between the PDS(K) (K = 2–4) and FCI energy are below $10−14$ a.u. The original Cioslowski CMX results start to diverge and show singularity on the potential energy surface when R_H–H approaches the dissociation limit (⁠$RH−H>1.50$ Å).

The singularity problem originates from the near-zero connected moments defining the denominators in the CMX energy expansion [see Eqs. (5) and (6)]. Mancini et al.⁵² have shown how to mitigate these singularities by re-summing the CMX expansion (i.e., the AMX approach) and introducing a new class of denominators corresponding to non-zero higher order moments. By a similar substitution derived from a Padé approximant, Knowles⁴⁷ showed that the singularity could be largely eliminated. This can be observed from H₂ results as shown in Figs. 3(a) and 3(b), where at the dissociation limit, the Knowles’s approach [Fig. 3(b)] is able to provide a smooth curve and move toward the FCI limit in an oscillatory manner, while divergence appears in the original Cioslowski’s CMX formalism after the second order, and additionally the singularity emerges in the original CMX(4) for R_H–H being ∼1.75 Å and ∼2.50 Å. Applications of QCMX to the H4 system are included in the supplementary material.

We further examine the performance of the CMX variants in terms of computing the ground state energy ranging from weak to strong correlation in the context of a single-impurity Anderson model (SIAM). As can be seen from Fig. 4, for fixed Hubbard repulsion U (U = 8), the exact correlation energy curve monotonically declines as the hybridization V becomes large. Among the low order CMX(K) (K = 1, 2) results, all the curves behave very similarly showing the same trend as the FCI curve [except the PDS(1) curve where the correlation contribution from the single-particle V terms operating on the present trial wave function is zero]. Among high order CMX(K) (K > 2) results, the original CMX diverges as the V becomes large, while the Knowles’s approach improves the original CMX results by showing slow convergence as the expansion order increases. On the other hand, the PDS approach shows fast convergence and basically reproduces the FCI results after the second order.

Note that in all the CMX calculations shown in Figs. 3 and 4, a simple form of the trial wave function (i.e., single determinant wave function) is used in either weak or strong correlation scenarios, and the commonly used unitary coupled-cluster Ansätze or, in general, the unitary operations for the state preparation as discussed in previous practices (see, e.g., Refs. 22 and 66-69), were not invoked. In other words, a crude trial wave function may be still useful by employing the proposed quantum algorithm to obtain highly accurate energies. This does not defy the importance of the state preparation at the dissociation or strong correlation limits but rather provides an alternatively great simplification for consideration when dealing with some quantum applications. Remarkably, different from the other two CMX variants, the PDS is, in principle, able to provide upper bounds for all the energy levels,⁶² and the upper bound for the ground state energy from the PDS approach is relatively tight even at the dissociation and strong correlation limits. Therefore, the PDS can further work with some unitary parameterization to improve the accuracy of the ground state energy computed at the low expansion order. As shown in Fig. 5, a unitary parameterization suggested in Ref. 68, U(θ) = exp(iθY₀X₁X₂X₃), can improve the fidelity of the trial wave function (see the supplementary material), thus helping the PDS(2) to well reproduce the FCI curve. As can be seen in Fig. 5(d), the energy deviation can be reduced on average by about three orders of magnitude when variational PDS(2) energies are minimized through rotating the trial wave function using the selected unitary operation U(θ) over the studied range of V values. Note here that employing the selected unitary parametrization alone cannot produce FCI energy for the four-qubit, 2-electron problem [here $EU(θ)=I1=Φ|U(θ)†HSIAMJWU(θ)|Φ=−4$ for |Φ⟩ = |0110⟩ and arbitrary θ, see Figs. 5(b) and 5(c)]. Similar results can be obtained with alternative parametrizations of U(θ) given by exp(iθX₀Y₁X₂X₃) or exp(iθY₀X₁). We have performed the CMX/PDS simulation for larger systems. For example, in the supplementary material, we have shown that the CMX/PDS quantum simulations using up to 12 qubits are able to target the full CI energies of the singlet and triplet states of the H₄ system described by different geometries and active spaces of various sizes, especially for the strongly correlated planar H₄ system.

We have also preliminarily studied the noise impact on the ideal CMX results employing the noise models that adopt reported values of Rigetti Aspen-1 QPUs and typical parameters for benchmarking error mitigation algorithms.^35,70 The results are shown in Figs. 6 and 7. As can be seen, the noise models slightly shifted the CMX(2)/PDS(2) results with respect to the idealized ones. Remarkably, for H₂ at larger R_H–H, the noise can “regularize” the idealized CMX(2) curve when the latter starts to disclose singularity due to the close-to-zero I₃ moments. For the Anderson model, even though the effect of noise is slightly bigger in the simulation with the correlated trial wave function, the quality of noisy PDS(2) results with the correlated reference is still better than the noisy PDS(2) results employing the |0110⟩ reference. Note that we do not include all sources of error in our emulation, and more comprehensive testing will be pursued in our future work.

The striking advantage offered by quantum CMX algorithms is their “insensitivity” with respect to the choice of the initial state. As long as the initial state is non-orthogonal to the exact state, the CMX expansions (especially the PDS approach) are capable of reproducing exact or near-to-exact energies. This feature could largely eliminate the challenging state-preparation process in quantum algorithms. As can be seen from the performance of the CMX variants in Figs. 3 and 4, instead of constructing a better trial wave function, the moment matrix M can be constructed to avoid singularity problems for quasi-degenerate situations. Finally, it is worth mentioning that, based on the same Krylov subspace, we have also performed the Lanczos approach (for up to fourth order) for the present studied systems, and did not find any observable difference between PDS(K) and Lanczos(K) results. However, noticing the numerical performance differences between different eigensolvers (or linear solvers) that are based on a same Krylov subspace in the classical computing, as well as previous classical comparisons between CMX and Lanczos results,^{47,51,52,57–60} we suppose performance differences between PDS and Lanczos methods would manifest for certain systems at certain scenarios. We believe that PDS and Lanczos methods might be complementary to each other for challenging many-body QC studies. Extensive testing for other models and molecular systems is currently underway.

We demonstrated the feasibility of a new QC algorithm based on calculating various types of CMX expansions. The discussed algorithm is robust in the sense of possible gate depth reduction and could be highly scalable with the increasing number of available qubits. Since the CMX algorithm offers a trade-off between the quality of the trial wave function and the rank of the moments required to achieve a high level of accuracy, we believe that its combination with the VQE approach, that is, utilizing the unitary coupled-cluster representation of the trail wave function, may provide a much needed formulation for improving the quality of VQE energies. Additionally, the flexibility of CMX re-summation techniques allows one to define the customized (or optimal) form of the expansion that avoids possible singularities associated with the near-zero values of the calculated connected moments. In our studies, we demonstrated that the PDS expansions provide superior results compared to other CMX expansions for the systems considered here, especially for the strongly correlated situations (characterized by a large R_H–H distance or large hybridization strength) studied with a relatively poor choice of trial wave function corresponding to a single Slater determinant. Our studies also indicate that the flexibility in choosing the reference function may be a regulatory factor to establish a compromise between accuracy and the effect of noise. In our future studies, we will explore the applicability of the PDS expansions for describing multi-configuration states and excited states. Of special interest will be integrating quantum CMX algorithms with downfolded Hamiltonians^71,72 and tensor decomposition techniques for defining effective interactions.^73–76 We will also explore the effect of noise in CMX/PDS simulations for larger systems.

See the supplementary material for the second quantized form of the H₂ and SIAM Hamiltonians, results of noisy CMX executions for H₂ and SIAM models, fidelity analysis of the trial state, and results for H₄ systems obtained in simulations with and without the noise.

This work was supported by the “Embedding QC into Many-body Frameworks for Strongly Correlated Molecular and Materials Systems” project, which is funded by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, the Division of Chemical Sciences, Geosciences, and Biosciences. All calculations have been performed at the Pacific Northwest National Laboratory (PNNL). PNNL is operated for the U.S. Department of Energy by the Battelle Memorial Institute under Contract No. DE-AC06-76RLO-1830.

The data that support the findings of this study are available from the corresponding author upon reasonable request.