Many-body techniques based on the double unitary coupled cluster (DUCC) ansatz can be used to downfold electronic Hamiltonians into low-dimensional active spaces. It can be shown that the resulting dimensionality reduced Hamiltonians are amenable for quantum computing. Recent studies performed for several benchmark systems using phase estimation (PE) algorithms for quantum computers demonstrated that these formulations can recover a significant portion of ground-state dynamical correlation effects that stem from the electron excitations outside of the active space. These results have also been confirmed in studies of ground-state potential energy surfaces using quantum simulators. In this letter, we study the effectiveness of the DUCC formalism in describing excited states. We also emphasize the role of the PE formalism and its stochastic nature in discovering/identifying excited states or excited-state processes in situations when the knowledge about the true configurational structure of a sought after excited state is limited or postulated (due to the specific physics driving excited-state processes of interest). In this context, we can view PE algorithms as an engine for verifying various hypotheses for excited-state processes and providing statistically meaningful results that correspond to the electronic state(s) with the largest overlap with a postulated configurational structure. We illustrate these ideas on examples of strongly correlated molecular systems, characterized by small energy gaps and high density of quasidegenerate states around the Fermi level.

There is significant interest in applying quantum computing techniques to describe and simulate chemical systems and processes.1–17 This approach brings hope to address the exponential barriers limiting the applicability of exact diagonalization procedures [or full configuration interaction (CI) methods] and also to provide access to complicated multiconfigurational electronic states, which often cannot be identified or described by vast classes of approximate methods used in routine simulations. Of special interest is the application to strongly correlated molecular systems characterized by small energy gaps between occupied and unoccupied orbitals where multiple electronic states that are defined by complex collective phenomena exist, usually involving higher than single excitations in the corresponding wavefunction expansions.

Even though impressive progress has been achieved in the development of wavefunction-based excited-state approaches such as complete-active-space perturbation theory (CASPT),18,19 multireference NEVPT,20,21 configuration interaction,22 equation-of-motion coupled cluster (EOMCC),23–29 multireference coupled cluster (MRCC),30–47 and the density matrix renormalization group (DMRG),48–51 problems with the description of complicated states dominated by high-rank excitations still exist. For example, in order to capture these states with EOMCC formalisms, one needs to include higher-than-double excitations.52 For doubly excited states, the “minimum” level of theory to tackle these states is EOMCC with singles, doubles, and triples (EOMCCSDT),53,54 although in several cases, it may not provide a quantitative level of accuracy.55 These problems commonly occur even for small molecular systems, and it is reasonable to expect that they may intensify for larger systems and strongly correlated systems (transition metal oxides, metal clusters, and actinides) with a high density of states located around the Fermi level. Recently, significant progress in addressing these problems has been achieved by integrating stochastic configuration interaction (CI) Quantum Monte Carlo (QMC) framework56,57 with deterministic EOMCC formulations.58 Another interesting aspect of modeling these complex excited states is the attainability of these states in situations where their initial configurational structure cannot be easily obtained to commence convergent iterative procedures.

Progress in the development of quantum algorithms may provide alternative solutions to these problems. The recent application of Variational Quantum Eigensolvers (VQEs)5,8,11,12,59–63 and Phase Estimation (PE)6,64–71 to excited states61,72–75 demonstrates that excited states can be effectively simulated on the quantum computers while at the same time bypassing the problems of conventional computing and approximate formulations.

In this paper, we present an excited state extension of recently developed techniques for active-space downfolding of the electronic Hamiltonian based on the double unitary coupled cluster (DUCC) transformation. This is combined with PE simulations available in the Microsoft Quantum Development Kit (QDK)74,76 to illustrate the excited-state version of DUCC formalism on the examples of H2 at equilibrium and stretched bond lengths and two H4 models: (1) trapezoidal H4, which is a popular benchmark system for studying quasidegenerate states,77 and (2) a linear H4 molecule, which was intensively studied in the context of singlet fission processes.78 For each of these systems, we perform PE simulations to characterize the structure of excited states using DUCC effective Hamiltonians. Moreover, we carry out a number of PE simulations to investigate the role and effect of different initial guesses, which are based on limited knowledge about the excited state of interest.

In Ref. 79, we introduced the unitary extension of the subsystem embedding subalgebra CC approach (SES-CC),80 which utilizes the double unitary CC expansion,

Ψ=eσexteσintΦ.
(1)

The characteristics of the expansion (1) is similar to the expansion discussed in the single-reference formulation of the active-space coupled cluster formalism81,82 (see also Refs. 83 and 84), which also utilizes the decomposition of the cluster operator into internal and external parts.

In analogy to Ref. 79, σint and σext are the anti-Hermitian operators (σint=σint and σext=σext) defined by excitations/de-excitations within and outside of active space, respectively. To be more precise, the amplitudes defining the σext operator must carry at least one inactive spin-orbital index. Using (1) in Schrödinger’s equation, one obtains equations for cluster amplitudes and the corresponding energy,

QeσinteσextHeσexteσintΦ=0,
(2)
ΦeσinteσextHeσexteσintΦ=E,
(3)

where Q is a projection operator on the space spanned by determinants orthogonal to the reference function Φ. In Eqs. (1)–(17) of this section, we consider the case of the exact limit (σint and σext include all possible excitations). In Ref. 79, we showed that when σint contains all possible excitations/de-excitations within the complete active space, the energy of system (2) can be obtained by diagonalizing the DUCC effective Hamiltonian,

H¯exteff(DUCC)eσintΦ=EeσintΦ,
(4)

where

H¯exteff(DUCC)=(P+Qint)H¯extDUCC(P+Qint)
(5)

and

H¯extDUCC=eσextHeσext.
(6)

In the above eigenvalue problem, the eσintΦ expansion defines a corresponding eigenvector and P and Qint are projection operators onto the reference function |Φ⟩ and excited determinants in the active space orthogonal to Φ, respectively.

To show this property, it is sufficient to introduce the resolution of identity eσinteσint to the left of the H¯extDUCC operator in

(P+Qint)H¯extDUCCeσint|Φ=E(P+Qint)eσint|Φ,
(7)

where we explicitly used the property of the eσint|Φ expansion

(P+Qint)eσint|Φ=eσint|Φ,
(8)

and to note that eσintH¯extDUCCeσint=eσinteσextHeσexteσint. Next, using matrix representation of the σint operator in the CAS space denoted as σint, this equation can be rewritten as

[eσint][y]=0,
(9)

where the first component of the [y] vector is equivalent to Φ|eσinteσextHeσexteσint|ΦE, while the remaining components correspond to projections of eσinteσextHeσexteσint|Φ onto excited configurations belonging to Qint. The [eσint] matrix is also nonsingular, which is a consequence of the formula

det(eσint)=eTr(σint)=1
(10)

and the anti-Hermitian characteristics of the σint matrix, i.e., Tr(σint) = 0 (where the real characteristics of σint cluster amplitudes is assumed). Given the nonsingular characteristics of the [eσint] matrix (see also Ref. 79), this proves the equivalence of these two representations.

The proof of the above property is not limited to the ground state and can be extended to any electronic state described by (1)–(3). Assuming that this ansatz can describe excited states, the DUCC effective Hamiltonian formalism can be used in the context of excited-state simulations. We will denote the general DUCC solution corresponding to the Kth state as

Ψ(K)=eσext(K)eσint(K)Φ,
(11)

where the Kth state energy can be obtained from diagonalizing the state-specific effective Hamiltonian,

H¯exteff(DUCC)(K)eσint(K)Φ=EKeσint(K)Φ,
(12)

where

H¯exteff(DUCC)(K)=(P+Qint)H¯extDUCC(K)(P+Qint)
(13)

and

H¯extDUCC(K)=eσext(K)Heσext(K).
(14)

Similar to the ground-state effective/downfolded Hamiltonians, the operators H¯exteff(DUCC)(K) are Hermitian and therefore amenable to the real-time simulation on a quantum computer. In analogy to the ground-state representation of DUCC, we will assume that

σint(K)=σint(K)    σext(K)=σext(K)
(15)

and

σint(K)=Sint(K)Sint(K),
(16)
σext(K)=Sext(K)Sext(K),
(17)

where Sint(K) and Sext(K) are CC-type cluster operators producing excitations within and outside the active space, respectively, when acting on the reference function Φ.

If the exact form of the operator σext(K) [or Sext(K)] is known, the effective Hamiltonian H¯exteff(DUCC)(K) can be diagonalized to find corresponding exact energy EK. In practice, in likeness to the ground-state formulation, we obtain an approximate model σext(K) operator by way of calculations with excited-state CC models. Additionally, we previously employed in ground-state DUCC79 the many-body form of the effective Hamiltonian H¯exteff(DUCC) driven by the perturbative analysis of the ground-state energy expansion. However, the same arguments cannot be invoked in the context of the excited-state variant of H¯extDUCC(K). Instead, in the analysis of the excited states, we use a basic expansion for H¯extDUCC(K) given by an expression involving a single commutator,

H¯extDUCC(K)=eσext(K)Heσext(K)H+[H,σext(K)].
(18)

In the future studies, we will explore the impact of the higher-order commutators on the quality of DUCC energies.

In this paper, we explore a simple strategy based on the utilization of the excited-state wavefunction in the equation-of-motion CC parameterization ΨKEOMCC,

ΨKEOMCC=RKeTΦ,
(19)

as a reference to extract the relevant information about σext(K). In Eq. (19), the cluster operator T satisfies the CC equations and the excitation operator RK (corresponding to the Kth excited state) is obtained through diagonalization of the similarity transformed Hamiltonian H¯=eTHeT. Since (11) represents a normalized state, in order to compare with the corresponding EOMCC expansion in the exact wavefunction limit, one needs to use a normalized form Ψ̃KEOMCC of (19),

Ψ̃KEOMCC=NKRKeTΦ=eσext(K)eσint(K)Φ,
(20)

where

NK=1ΨKEOMCCΨKEOMCC.
(21)

In the simplest approximate variants, we use the EOMCCSD approximations (EOMCC with singles and doubles24) to extract singly excited [Sext, 1(K)] and doubly excited [Sext, 2(K)] components of the Sext(K) operator,

Sext,1(K)ΦQext,1Ψ̃KEOMCC,
(22)
Sext,2(K)ΦQext,2Ψ̃KEOMCC,
(23)

where Qext, 1 and Qext, 2 are projection operators on subspaces of singly and doubly excited external excitations, respectively. In (22) and (23), we approximate Ψ̃KEOMCCΨKEOMCCSD(A), which is defined in the following way:

ΨKEOMCCSD(A)=(P+Q1+Q2)(RK,0+RK,1+RK,2)eT1+T2Φ,
(24)

where Q1 and Q2 operators are projection operators on spaces of singly and doubly excited configurations, RK,i (i = 0, 1, 2) represent EOMCC excitation operators producing i-tuply excited configurations when acting onto reference functions, and T1 and T2 are the singly and doubly excited cluster operators, respectively. In our approximation, referred to as the DUCC-ex(K) approach, we use T1 and T2 from standard CCSD calculations and RK,i (i = 0, 1, 2) from the standard EOMCCSD diagonalization procedure. Consequently, Sext, 1(K) and Sext, 2(K) take the following form:

Sext,1(K)ΦNK(A)Qext,1(RK,0T1+RK,1)Φ,
(25)
Sext,2(K)ΦNK(A)Qext,2RK,0T2+12T12+RK,1T1+RK,2Φ,
(26)

where

NK(A)=1ΨKEOMCCSD(A)ΨKEOMCCSD(A).
(27)

As in the ground-state case, the DUCC-ex(K) approximation is defined by the length of the commutator expansion and the source of the σext(K) amplitudes. Since the EOMCCSD approximation is used to define σext(K), one should expect that this scheme works in cases where the EOMCCSD approach delivers reliable results and active-space used in the state-selective DUCC-ex(K) formalism capable of capturing main configurations corresponding to the Kth state.

Our simulation of the DUCC Hamiltonian consists of two broad components: (1) prepare a suitable classical representation of the DUCC Hamiltonian and the CCSD and EOMCCSD state to be simulated and (2) perform a quantum simulation of the system using quantum algorithms. The proposed workflow is analogous to the one introduced by Wang et al.,72 up to a difference in initial state(s), the form of the diagonalized Hamiltonian, and the numerical simulation of the quantum algorithm.

The CCSD and EOMCCSD calculations are performed for the ground and excited states using Tensor Contraction Engine (TCE)85 CCSD and EOMCCSD implementations available in NWChem.86 The quantum simulations are performed using the QDK-NWChem interface and targets a classical emulation of a quantum computer. Here, NWChem generates a yaml file containing the orbital approximation of the Hamiltonian and quantum state. Subsequently, eigenenergy estimates are obtained using QDK implementations of QPE and simulation based on either Trotterization66 or qubitization.76 

We now describe these components and estimate the quantum resources required. Note that the number of quantum gates executed is polynomial in the number of orbitals. This is efficient on a quantum computer but of course requires exponential time on a classical simulation.

The following steps 1–4 are performed using classical resources.

  1. Form|Ψ̃KEOMCC. This state approximates the targeted excited eigenstate, call it ΨK, and is obtained by normalizing the |ΨKEOMCCSD(A) vector and forming Sext,1(K) and Sext,2(K) amplitudes as defined by (25) and (26).

  2. Form theH¯extDUCC(K)Hamiltonian in the particle-hole representation. As in Ref. 79, we calculate the normal product form of H¯extDUCC(K)(18), which is approximated by the scalar, one-, and two-body interactions, i.e.,

H¯extDUCC(K)(Γ(K))scalar+P,Qγ(K)QPN[aPaQ]ph+14P,Q,R,Sγ(K)RSPQN[aPaQaSaR]ph,
(28)

where N[⋯]ph denotes the particle-hole (ph) normal product form of a given operator expression and P, Q, R, S, … denote general active spin-orbital indices. The algebraic form of the γ(K)QP and γ(K)RSPQ tensors for the situation when all occupied orbitals are deemed active is given in Tables I and II.

  1. TransformH¯extDUCC(K)from particle-hole representation(28)to Fermi vacuum representation.

H¯extDUCC(K)PQχ(K)QPaPaQ+14P,Q,R,Sχ(K)RSPQaPaQaSaR
(29)

using identities from Table IV of Ref. 79. In analogy to Ref. 79 as an input for the QDK, we are using a symmetric orbital approximation of the spin orbital χ tensor, where it is assumed that Mulliken orbital-type dressed integrals (PQ|RS) (as well as (PQ|SR), (QP|RS), …, (SR|QP)) are obtained from χRαSβPαQβ (where P, Q, R, and S are the orbitals corresponding to active spin orbitals P, Q, R, and S, respectively). In general, the number M of unique integrals (accounting for the symmetries (PQ|RS) = (PQ|SR), (QP|RS), …, (SR|QP)) is much fewer than the worst-case of M=O(N4), where N is the number of active-space spin-orbitals.

  1. Transform Fermi vacuum representation to Jordan-Wigner representation. The targeted quantum computer applies quantum gates on qubits rather than fermions. Thus, we map fermion operators to Pauli operators {XP, YP, ZP} acting on qubits P ∈ [0, N − 1]. It suffices for our purposes to use the Jordan-Wigner representation,

aPZ0ZP1(XpiYp)/2,aPZ0ZP1(Xp+iYp)/2.
(30)

One may verify that the transformed operators satisfy the expected fermion anticommutation relationships. Thus, the DUCC Hamiltonian is of the form

H¯extDUCC(K)j=0Θ(M)αjPj,
(31)

where Pj are Pauli operators with coefficients αj.

TABLE I.

The algebraic form of the γQP(K) and γ(K)RSPQ stemming from the normal product form of the bare Hamiltonian in the active space HN=PQfQPN[aPaQ]ph+14P,Q,R,SvRSPQ[aPaQaSaR]ph. For simplicity, γQP(K) and γ(K)RSPQ are denoted as γQP and γRSPQ.

AmplitudeExpressionAmplitudeExpression
HN term 
γAB=fAB γIJ=fIJ 
γAI=fAI γIA=fIA 
γIABC=vIABC γIJKA=vIJKA 
γABCI=vABCI γKAIJ=vKAIJ 
γIAJB=vIAJB γABCD=vABCD 
γIJKL=vIJKL γABIJ=vABIJ 
γIJAB=vIJAB   
AmplitudeExpressionAmplitudeExpression
HN term 
γAB=fAB γIJ=fIJ 
γAI=fAI γIA=fIA 
γIABC=vIABC γIJKA=vIJKA 
γABCI=vABCI γKAIJ=vKAIJ 
γIAJB=vIAJB γABCD=vABCD 
γIJKL=vIJKL γABIJ=vABIJ 
γIJAB=vIJAB   
TABLE II.

Algebraic form of γ(K)QP and γ(K)RSPQ amplitudes in Eq. (28) (continued). Amplitudes defining Sext,1(K) and Sext,2(K) operators are denoted by sia and sijab (i, j, … and a, b, … are generic occupied and virtual spin-orbitals, respectively), while amplitudes defining Sext,1(K) and Sext,2(K) are denoted as sai and sabij, respectively. By definition of the external parts of S(K) and S(K), all s-amplitudes that carry active spin-orbital indices only disappear. These terms pertain to active spaces that contain all correlated occupied orbitals and a subset of the virtual ones. For simplicity we assume a restricted/unrestricted Hartree-Fock (RHF/UHF) reference |Φ⟩, where all nondiagonal Fock matrix elements disappear. The Einstein summation convention is invoked. For simplicity, γQP(K) and γ(K)RSPQ are denoted as γQP and γRSPQ.

AmplitudeExpression
(HNText)C,open term 
γAB+=veAMBsMe12veAMNsMNeB 
γIJ+=veIMJsMe+12vefMJsMIef 
γAI+=veAMIsMe 
γIA+=feAsIe+veIMAsMe12veIMNsMNeA+12vefMAsMIef 
γIABC+=veABCsIeveAMBsMIeC+veAMCsMIeB 
γIJKA+=veJKAsIeveIKAsJe+12vefKAsIJefveJMKsMIeA+veIMKsMJeA 
γKAIJ+=veAIJsKe 
γIAJB+=veAJBsIeveAMJsMIeB 
γIJKL+=veJKLsIeveIKLsJe+12vefKLsIJef 
γIJAB+=veJABsIeveIABsJe+feAsIJeBfeBsIJeA+12vefABsIJef 
veJMAsMIeB+veIMAsMJeB+veJMBsMIeAveIMBsMJeA 
(TextHN)C,open term 
γBA+=vMBeAseM12vMNeAseBMN 
γJI+=vMJeIsMe+12vMJefsefMI 
γIA+=vMIeAsMe 
γAI+=fAeseI+vMAeIseM12vMNeIseAMN+12vMAefsefMI 
γBCIA+=vBCeAseIvMBeAseCMI+vMCeAseBMI 
γKAIJ+=vKAeJseIvKAeIseJ+12vKAefsefIJvMKeJseAMI+vMKeIseAMJ 
γIJKA+=vIJeAseK 
γJBIA+=vJBeAseIvMJeAseBMI 
γKLIJ+=vKLeJseIvKLeIseJ+12vKLefsefIJ 
γABIJ+=vABeJseIvABeIseJ+fAeseBIJfBeseAIJ+12vABefsefIJ 
vMAeJseBMI+vMAeIseBMJ+vMBeJseAMIvMBeIseAMJ 
AmplitudeExpression
(HNText)C,open term 
γAB+=veAMBsMe12veAMNsMNeB 
γIJ+=veIMJsMe+12vefMJsMIef 
γAI+=veAMIsMe 
γIA+=feAsIe+veIMAsMe12veIMNsMNeA+12vefMAsMIef 
γIABC+=veABCsIeveAMBsMIeC+veAMCsMIeB 
γIJKA+=veJKAsIeveIKAsJe+12vefKAsIJefveJMKsMIeA+veIMKsMJeA 
γKAIJ+=veAIJsKe 
γIAJB+=veAJBsIeveAMJsMIeB 
γIJKL+=veJKLsIeveIKLsJe+12vefKLsIJef 
γIJAB+=veJABsIeveIABsJe+feAsIJeBfeBsIJeA+12vefABsIJef 
veJMAsMIeB+veIMAsMJeB+veJMBsMIeAveIMBsMJeA 
(TextHN)C,open term 
γBA+=vMBeAseM12vMNeAseBMN 
γJI+=vMJeIsMe+12vMJefsefMI 
γIA+=vMIeAsMe 
γAI+=fAeseI+vMAeIseM12vMNeIseAMN+12vMAefsefMI 
γBCIA+=vBCeAseIvMBeAseCMI+vMCeAseBMI 
γKAIJ+=vKAeJseIvKAeIseJ+12vKAefsefIJvMKeJseAMI+vMKeIseAMJ 
γIJKA+=vIJeAseK 
γJBIA+=vJBeAseIvMJeAseBMI 
γKLIJ+=vKLeJseIvKLeIseJ+12vKLefsefIJ 
γABIJ+=vABeJseIvABeIseJ+fAeseBIJfBeseAIJ+12vABefsefIJ 
vMAeJseBMI+vMAeIseBMJ+vMBeJseAMIvMBeIseAMJ 

We obtain estimates of excited state energies using a phase estimation algorithm on a quantum circuit approximating the unitary time-evolution operator eiHt, acting on a trial wavefunction ΨTrial. As abstractly illustrated in Fig. 1, the trial wavefunction may have overlap with several states (Ψ1 and Ψ2) in Fig. 1 and be orthogonal to others (Ψ3) in Fig. 1. Each run of PE returns an approximate energy estimate for a random eigenstate with probability equal to the amount of overlap |⟨ΨjTrial⟩|2 of the trial wavefunction with that state’s corresponding wavefunction. As an example, Fig. 2 illustrates just such a distribution of energies for a strongly correlated system. Importantly, one will always obtain an energy estimate even if this overlap is small. This is sufficient for the purpose of characterizing the excited state spectrum.

FIG. 1.

An abstract representation of how a trial wavefunction may have overlap Pj = |⟨ΨjTrial⟩|2 with several states (Ψ1 and Ψ2) and be orthogonal to others (Ψ3). In this figure, Pi represents the probability of sampling the ith state in a simulation with the PE algorithm.

FIG. 1.

An abstract representation of how a trial wavefunction may have overlap Pj = |⟨ΨjTrial⟩|2 with several states (Ψ1 and Ψ2) and be orthogonal to others (Ψ3). In this figure, Pi represents the probability of sampling the ith state in a simulation with the PE algorithm.

Close modal
FIG. 2.

A typical distribution of energies for a strongly correlated system obtained from several simulations with the PE algorithm. This particular spread of energies corresponds to the DUCC-ex(21Ag) results for H2 at a stretched geometry (RH–H = 10 a.u.) (see Table IV).

FIG. 2.

A typical distribution of energies for a strongly correlated system obtained from several simulations with the PE algorithm. This particular spread of energies corresponds to the DUCC-ex(21Ag) results for H2 at a stretched geometry (RH–H = 10 a.u.) (see Table IV).

Close modal

We now outline the quantum gate cost of these algorithms.

1. Trial wavefunction preparation

As described in Ref. 74, the distribution of energies for the ground and excited states from PE is determined by the Hamiltonian and a trial wavefunction composed of a superposition of m Slater determinants. In the QDK, this state can be prepared nondeterministically with success probability O(1/N2) using O(mN) quantum gates. With quantum amplitude amplification,87 the success probability can be boosted close to unity using a total of O(mN2) quantum gates. Note that alternate techniques may efficiently synthesize on a quantum computer trial wavefunction with exponentially many components, such as using the UCC formalism.79 However, this is unnecessary for PE methods when the exact wavefunction is dominated by only a few Slater determinants. Unlike VQE techniques, where obtaining energy estimates that are correct to chemical accuracy requires the overlap to also be unity up to at least chemical accuracy, PE methods only need the overlap to be polynomially large, say, |⟨ΨjTrial⟩|2 = Ω(1/poly(N)).

For this purpose, a wide array of excited-state EOMCC formulations implemented in NWChem can be used to provide necessary information about the configurational structure of excited states in molecular systems and therefore complement a strategy based on the utilization multiconfigurational self-consistent field (MCSCF) wave functions (see Ref. 72). The following part of the discussion in this section will however deal with the hypothetical situation when our knowledge of the initial state, for whatever reason, is not perfect and we cannot precisely determine its structure, instead, we can only hypothesize about the structure of the sought state. This is a situation which may occur in studies of strongly correlated systems. Through the repeated application of PE where a carefully chosen trial wavefunction is reinitialized each time, one can obtain a distribution of energies for several states of interest.

The stochastic nature of PE algorithms in providing information on a distribution of targeted states is unlike VQE algorithms that only provide energy estimates for a single targeted state. As which single state to target may not be obvious a priori, PE methods enable the discovery and investigation of exotic and novel states that are unobtainable with conventional computing and current approximate formulations.

2. Phase estimation

Phase estimation algorithms, such as quantum phase estimation,87 return a random eigenphase θj of an arbitrary unitary U=jeiθjΦjΦj with eigenvectors Φj. The dominant cost of this algorithm is in applying controlled-U some number T times. The dominant error in the PE energy estimate arises from an analog of the time-energy uncertainty relation. Each application of phase estimation estimates θj with a standard deviation of σθ = Θ(T−1), and the robust PE algorithm68 we use in the QDK empirically achieves σθ ≈ 2/T and uses only one ancillary qubit.

Implementing U exactly can be difficult as it must be expressed in terms of primitive quantum gates that may be implemented on a universal quantum computer. Thus, we permit a maximum simulation error ϵθ and instead output an operation U′ such that

UUϵθ.
(32)

This implies that the eigenvalues and eigenphases θj of U′ differ from U by at most ϵθ. The overall error in estimating θj is then a random error of σθ and a systematic shift of ϵ.

3. Unitary time-evolution

In order to estimate energy with a standard deviation σE, PE algorithms query the unitary time-evolution operator eiHΔ=jeiEjΔψjψj by the DUCC Hamiltonian. In order to avoid ambiguity in the principle value of the phase, Δ should be chosen so that for the eigenstate of energy, the range of energies |EmaxEmin|Δ ≤ 2π. Thus, obtaining an energy estimate with standard deviation σE requires σEΔ ≈ 2/T. In many cases, the step-size Δ determines the error of approximation. Suppose we approximate eiHΔ by a unitary operator with error,

eiHΔUϵθ.
(33)

Then, this contributes a systematic error in energy of ϵE = ϵθ/Δ. At present, the two most popular methods for constructing the operator U are Trotter-Suzuki methods66 and qubitization,76 both of which are available in the QDK.

In the first-order Trotter-Suzuki method, real-time evolution by the DUCC Hamiltonian is approximated by

U1(Δ)=jΘ(M)eiαjΔPj,eiHΔU=O|Δ|2α12,
(34)

where α1=j|αj|. The second-order method with error O|Δ|3α13 applies a symmetric sequence U2(Δ) that is U1(Δ/2) followed by its component exponentials in reverse order. In order for ϵE to be commensurate with σE, this implies that

Δ=OσE/α13,T=O1ΔσE.
(35)

Higher-order j methods with error O5j/21|Δ|j+1α1j+1 are known but can be impractical due to the exponential prefactor. It is also worth noting that these error bounds represent worst-case behavior. In empirical studies,13 the step-size can be much larger than suggested by Eq. (35).

The overall gate cost of obtaining a single energy estimate is then T times the gate cost of U2(Δ). The overall space requirement is N qubits for representing a system with N spin-orbitals in the Jordan-Wigner representation and one more for phase estimation. Within fault-tolerant models of quantum computation, the Clifford gates that conjugate to Pj to a single-qubit Pauli Z are considered cheap and dominant cost of U1(Δ) is in applying Θ(M) arbitrary single-qubit Z-rotations. Alternate methods based on quantum signal processing88 or qubitization76,89,90 achieve even smaller error. We do not consider them further as they require O(logN) additional ancillary qubits, which is not conducive to a numerical simulation.

To summarize, this is the scheme of phase estimation combined with unitary time-evolution and is efficient as the number of quantum gates is polynomial in the number of spin-orbitals. Furthermore, speedups in gate depth and count could be obtained by using more compressed Hamiltonians. However, progress in this field is contingent on the progress achieved in understanding the many-body structure of downfolded Hamiltonians and quantifying its errors, such as in the representation of the σ(K)ext operator.

We performed a series of numerical tests with the QDK simulator74 for three systems characterized by strong ground-state correlation effects and low-lying singlet excited states defined by complicated multireference configurational structures: (1) H2 systems for H–H separation corresponding to the equilibrium (RH–H = 1.4008 a.u.) and a stretched geometry (RH–H = 10 a.u.), (2) H4 system in the trapezoidal configuration corresponding to geometrical parameter α equal to 0.001 [see Fig. 3(a)]—for this geometry, one can observe a strong quasidegeneracy of low-lying electronic states, and (3) linear form of H4 used in studies of singlet fission processes [see Fig. 3(b)]. Special attention is paid to singlet doubly excited states that pose a significant challenge to existing many-body methodologies. For systems considered here we employ the cc-pVTZ basis set91 and we correlate all electrons in all calculations. Since the QDK cannot exploit spatial symmetry, all PE simulations were performed using C1 symmetry. For simplicity, we will refer to the D2h, C2v, and D2h irreducible representations of the H2, H4 trapezoidal, and H4 linear systems, respectively, when reporting full configuration interaction (FCI) numbers obtained with NWChem86 (or equivalently, EOMCCSD or EOMCCSDTQ results for the H2 and H4 models, respectively).

FIG. 3.

The trapezoidal (a) and linear (b) H4 models employed in excited-state simulations (see text for details).

FIG. 3.

The trapezoidal (a) and linear (b) H4 models employed in excited-state simulations (see text for details).

Close modal

For H2 models, we examine the performance of the bare and DUCC-ex transformed Hamiltonians for excited states, which have a significant (A) or partial (PA) amount of leading characteristic excitations in the active space, which contains one occupied and three lowest-lying RHF orbitals. Energies of these states for RH–H = 1.4008 a.u. and RH–H = 10.0 a.u. are shown in Tables III and IV, respectively. Since the lowest-lying doubly excited states are of Ag symmetry, in the forthcoming analysis, we will entirely focus on the states corresponding to this symmetry. In Table III, we consider three states 21Ag, 31Ag, and 51Ag, which are either comprehensively described within the active space (21Ag) or have a significant overlap with it (31Ag and 51Ag states). The 21Ag is dominated by Φ13 and Φ1¯3¯ single excitations, while the two other states have non-negligible out-of-active-space components. For example, the 31Ag state is dominated by Φ11¯22¯ excitation but contains important component originating in the Φ17 and Φ1¯7¯ excitations. In a similar way, the 51Ag state is dominated by Φ11¯33¯ excitations, but it also contains important contributions from Φ11¯37¯ and Φ11¯73¯ configurations, which do not belong to the active space.

TABLE III.

Total energies of low-lying singlet excited states of the H2 system (RH–H = 1.4008 a.u.) of the Ag symmetry (all classical calculations for H2 have been performed using the D2h symmetry group). The values in parentheses are errors relative to FCI. In all calculations, the restricted Hartree-Fock orbitals were used.

StateDUCC-ex(21Ag)DUCC-ex(31Ag)DUCC-ex(51Ag)Bare Hamiltonian
(characteristics)aFCIb SD: 12(Φ13+Φ1¯3¯)SD: Φ11¯22¯SD: Φ11¯33¯ SD: Φ11¯22¯
21Ag −0.5487 −0.5478 ± 0.0066 −0.5197 ± 0.0019 −0.5034 ± 0.0030 −0.5306 ± 0.0017 
(A)  (0.0009) (0.0290) (0.0453) (0.0181) 
31Ag −0.1210 … 0.0644 ± 0.0017 −0.0696 ± 0.0016 −0.0622 ± 0.0018 
(PA)   (0.0566) (0.0514) (0.0588) 
51Ag 0.2310 0.2914 ± 0.0021 0.2819 ± 0.0065 0.2671 ± 0.0031 0.2842 ± 0.0014 
(PA)  (0.0604) (0.0508) (0.0361) (0.0532) 
StateDUCC-ex(21Ag)DUCC-ex(31Ag)DUCC-ex(51Ag)Bare Hamiltonian
(characteristics)aFCIb SD: 12(Φ13+Φ1¯3¯)SD: Φ11¯22¯SD: Φ11¯33¯ SD: Φ11¯22¯
21Ag −0.5487 −0.5478 ± 0.0066 −0.5197 ± 0.0019 −0.5034 ± 0.0030 −0.5306 ± 0.0017 
(A)  (0.0009) (0.0290) (0.0453) (0.0181) 
31Ag −0.1210 … 0.0644 ± 0.0017 −0.0696 ± 0.0016 −0.0622 ± 0.0018 
(PA)   (0.0566) (0.0514) (0.0588) 
51Ag 0.2310 0.2914 ± 0.0021 0.2819 ± 0.0065 0.2671 ± 0.0031 0.2842 ± 0.0014 
(PA)  (0.0604) (0.0508) (0.0361) (0.0532) 
a

Characteristics of the electronic states are determined by the active-space contribution: A, dominated by active-space configurations; PA, dominated by configurations not belonging to active space.

b

Full configuration interaction calculations were performed using all 30 molecular orbitals.

TABLE IV.

Total energies of low-lying singlet excited states of the H2 system (RH–H = 10 a.u.) of the Ag symmetry (all classical calculations for H2 have been performed using the D2h symmetry group). The values in parentheses are errors relative to FCI. In all calculations, restricted Hartree-Fock orbitals were used.

StateDUCC-ex(21Ag) DUCC-ex(31Ag) DUCC-ex(51Ag) DUCC-ex(71Ag)Bare Hamiltonian
(characteristics)aFCIb SD: Φ11¯22¯ SD: 12(Φ13+|Φ1¯3¯) SD: Φ11¯33¯ SD: Φ11¯44¯ SD: Φ11¯22¯
21Ag −0.5981 −0.6047 ± 0.0018 −0.5847 ± 0.0018 −0.5517 ± 0.0018 −0.5610 ± 0.0013 −0.5803 ± 0.0019 
(A)  (−0.0065) (0.0134) (0.0464) (0.0371) (0.0178) 
31Ag −0.4873 −0.4825 ± 0.0045 −0.4848 ± 0.0050 −0.4747 ± 0.0043 −0.4785 ± 0.0035 −0.4795 ± 0.0043 
(A)  (0.0048) (0.0025) (0.0126) (0.0088) (0.0078) 
51Ag −0.1040 −0.0812 ± 0.0018 −0.0779 ± 0.0017 −0.0625 ± 0.0018 −0.0592 ± 0.0020 −0.0763 ± 0.0017 
(A)  (0.0228) (0.0261) (0.0415) (0.0448) (0.0277) 
61Ag 0.0238 … … 0.0405 ± 0.0016 0.0365 ± 0.0017 0.0368 ± 0.0024 
(A)    (0.0167) (0.0127) (0.0130) 
71Ag 0.2326 0.2525 ± 0.0023 0.2545 ± 0.0019 0.2623 ± 0.0019 0.2584 ± 0.0017 0.2556 ± 0.0020 
(A)  (0.0199) (0.0219) (0.0297) (0.0258) (0.0230) 
StateDUCC-ex(21Ag) DUCC-ex(31Ag) DUCC-ex(51Ag) DUCC-ex(71Ag)Bare Hamiltonian
(characteristics)aFCIb SD: Φ11¯22¯ SD: 12(Φ13+|Φ1¯3¯) SD: Φ11¯33¯ SD: Φ11¯44¯ SD: Φ11¯22¯
21Ag −0.5981 −0.6047 ± 0.0018 −0.5847 ± 0.0018 −0.5517 ± 0.0018 −0.5610 ± 0.0013 −0.5803 ± 0.0019 
(A)  (−0.0065) (0.0134) (0.0464) (0.0371) (0.0178) 
31Ag −0.4873 −0.4825 ± 0.0045 −0.4848 ± 0.0050 −0.4747 ± 0.0043 −0.4785 ± 0.0035 −0.4795 ± 0.0043 
(A)  (0.0048) (0.0025) (0.0126) (0.0088) (0.0078) 
51Ag −0.1040 −0.0812 ± 0.0018 −0.0779 ± 0.0017 −0.0625 ± 0.0018 −0.0592 ± 0.0020 −0.0763 ± 0.0017 
(A)  (0.0228) (0.0261) (0.0415) (0.0448) (0.0277) 
61Ag 0.0238 … … 0.0405 ± 0.0016 0.0365 ± 0.0017 0.0368 ± 0.0024 
(A)    (0.0167) (0.0127) (0.0130) 
71Ag 0.2326 0.2525 ± 0.0023 0.2545 ± 0.0019 0.2623 ± 0.0019 0.2584 ± 0.0017 0.2556 ± 0.0020 
(A)  (0.0199) (0.0219) (0.0297) (0.0258) (0.0230) 
a

Characteristics of the electronic states are determined by the active-space contribution: A, dominated by active-space configurations and PA, dominated by configurations not belonging to active space.

b

Full configuration interaction calculations were performed using all 30 molecular orbitals.

When the bare Hamiltonian for H2 at equilibrium is diagonalized in the active space, all three excited states can be observed. To improve the distribution and get a better statistical sampling for the excited states (particularly 31Ag), the initial wavefunction was chosen to be the Φ11¯22¯ configuration, although it is not necessary to track these excited states. Unfortunately, there is a significant difference between the FCI energies and energies obtained by diagonalizing the bare Hamiltonian with errors of over 18 mhartree for 21Ag and 50 mhartree for 31Ag and 51Ag.

The DUCC-ex formalism for the 21Ag state (DUCC-ex(21Ag)), where the similarity transformation is driven by σext(K) corresponding to the 21Ag state [see (23) and (24)], provides an excellent agreement with the FCI energies with an error less than 1 mhartree, a substantial improvement over the energy obtained with the bare Hamiltonian. From Table III, one can note that by using the trial wavefunction 12(Φ13+Φ1¯3¯), there is a chance to describe not only target 21Ag (which is the physical solution given the state-specific nature of the DUCC-ex approach) but also other states, which have nonzero overlap with the trial wavefunction (in this case 51Ag). For the 31Ag and 51Ag states, the effect of important out-of-active-space Slater determinants needs to be determined perturbatively and is only approximately captured by the corresponding σext(K) operators. Consequently, the DUCC-ex(31Ag) and DUCC-ex(51Ag) results are less accurate compared to the DUCC-ex(21Ag) case. However, the results obtained with the respective DUCC-ex(31Ag) and DUCC-ex(51Ag) Hamiltonians are still better compared to the diagonalization of the bare-Hamiltonian. In addition, in the simulations for 31Ag and 51Ag, all three states are observed when the initial wavefunctions are taken to be the leading excitation for the corresponding state.

Compared to the equilibrium geometry, the H2 system corresponding to the RH–H = 10.0 a.u. is characterized by a larger number of excited states with a significant overlap with active space. For RH–H = 10.0 a.u., one can identify five excited states of Ag symmetry falling into this category. As shown in Table IV, the low-lying excited states can be efficiently described by the DUCC-ex formalism. For example, the DUCC-ex(21Ag) and DUCC-ex(31Ag) Hamiltonians provide very accurate estimates for the 21Ag and 31Ag states, which is quite remarkable given the size of the active space. In both cases, one can observe a significant improvement of the results obtained by the diagonalization of the bare Hamiltonian in the same active space, which provides a good illustration of the efficiency of the DUCC-ex formalism, even in a simple case corresponding to single commutator expansion of the DUCC-ex downfolded Hamiltonian. For higher excited states (51Ag, 61Ag, and 71Ag), the DUCC-ex results are less accurate, but still comparable to the results obtained through the diagonalization of the bare Hamiltonian in the active space. This behavior is justified given the complex excitation manifold describing the higher-lying states and ought to be resolved with improved approximations for σext(K).

Table IV also provides an excellent illustration of the fact that the trial wavefunction can be used to probe various electronic states in quantum simulations (see also Fig. 2). For example, using the trial state Φ11¯22¯ for PE simulations of the DUCC-ex(21Ag) Hamiltonian, one can obtain a statistically meaningful presence of states different from a challenging doubly excited 21Ag state, i.e., 11Ag (not shown in Table IV), 31Ag, 51Ag, and 71Ag. This fact is associated with the strong quasidegeneracy of the ground-state and low-lying excited states with contribution from the Φ11¯22¯ Slater determinant. It is also worth mentioning that all states considered here are either purely doubly excited states or states of mixed single- and doubly excited characteristics, which usually pose a significant challenge for approximate EOMCC formulations. Analogous situations are naturally occurring in strongly correlated systems, and the fact that PE can deal with these challenging problems provides yet another argument in favor of developing quantum algorithms for excited states. Similar behavior can also be seen in the case of the active-space representation of the bare Hamiltonian. Although the DUCC-ex approach is state-specific in its nature, the appearance of other physically interpretable excited states in the spectra of DUCC-ex Hamiltonians should be explored further.

We exploited the state-specific nature of the DUCC-ex approach in an attempt to describe the lowest-lying fully symmetric 21A1 and 21Ag states (vide infra) of the H4 system in trapezoidal (α = 0.001) [H4(a) system] and linear configuration [H4(b) system] shown in Fig. 3. For both systems, FCI results were obtained by running EOMCCSDTQ calculations using C2v and D2h symmetries, respectively. The active space used to construct DUCC effective Hamiltonians in both configurations includes the two occupied orbitals and the five lowest-lying virtual orbitals. While the ground state for H4(a) discloses strong quasidegeneracy effects between Φ and Φ22¯33¯ Slater determinants, for the H4(b) system, these effects slightly weaker. In the description of the 21A1 state of the H4(a) system, the dominant role is played by the Φ22¯33¯ determinant, which almost entirely dictates the corresponding wavefunction expansion. The 21Ag state of H4(b) reveals a more multiconfigurational structure where determinants Φ13, Φ1¯3¯, Φ24, Φ2¯4¯, Φ1234, Φ1¯2¯3¯4¯, Φ11¯33¯, Φ22¯33¯, and Φ22¯44¯ all have a non-negligible contribution to the wavefunction. As shown in Table V, we are able to inspect the lowest-lying doubly excited state for both systems with the bare Hamiltonian. It is important to note that unlike the H2 systems, one must use an initial guess that contains the Φ22¯33¯ determinant with some significant weight. Otherwise, only the ground state will be observed. Even though these states can be tracked with the bare Hamiltonians, they come with errors of 73 mhartree for the H4(a) system and 38 mhartree for the H4(b) system. As shown in Table V, the DUCC-ex(21A1) approach can very efficiently reproduce the FCI 21A1 energy with relative error less than 10 mhartree. As in the case of H2 system using Φ22¯33¯ as an initial state for quantum simulations, one can observe a statistical presence of states which have a non-negligible overlap with this determinant. It is interesting to note that for the H4(b) system, the DUCC-ex results for the 21Ag are less accurate and they provide comparable results to the bare Hamiltonian energy. In our opinion, this is associated with the strong multiconfigurational characteristics of the corresponding wavefunction, where simple approximations (25) and (26) may not provide the desired level of accuracy. This statement can be supported by the analysis of lowest-lying orbitals (shapes of the active orbitals used in our calculations are shown in Fig. 4). As shown in Fig. 5, orbital energies of orbitals not included in the active space (orbitals 8,9, …) for H4(a) system grow faster compared to the H4(b) case, where energies of nonactive orbitals reach a plateau involving orbitals 8–13. This behavior suggests that much larger active spaces are needed for modeling H4(b) with the DUCC-ex formalism.

TABLE V.

Total energies of low-lying singlet states of the H4 system of the A1 [H4(a)] or Ag [H4(b)] symmetry. The classical calculations for H4 have been performed using either the C2v [H4(a)] or D2h [H4(b)] symmetry group. The values in parentheses are errors relative to FCI. In all calculations, the restricted Hartree-Fock orbitals were used.

StateBare HamiltonianDUCC-ex(21A1)
(characteristics)aFCIbSD: Φ22¯33¯SD: Φ22¯33¯
α = 0.001 [H4(a) system] 
21A1 −2.0280 −1.9546 ± 0.0037 −2.0184 ± 0.0033 
(A) (0.0734) (0.0096) 
Linear [H4(b) system] 
21Ag −1.9901 −1.9513 ± 0.0052 −1.9458 ± 0.002894 
(A) (0.0388) (0.0443) 
StateBare HamiltonianDUCC-ex(21A1)
(characteristics)aFCIbSD: Φ22¯33¯SD: Φ22¯33¯
α = 0.001 [H4(a) system] 
21A1 −2.0280 −1.9546 ± 0.0037 −2.0184 ± 0.0033 
(A) (0.0734) (0.0096) 
Linear [H4(b) system] 
21Ag −1.9901 −1.9513 ± 0.0052 −1.9458 ± 0.002894 
(A) (0.0388) (0.0443) 
a

Characteristics of the electronic states are determined by the active-space contribution: A, dominated by active-space configurations and PA, dominated by configurations not belonging to active space.

b

Full configuration interaction calculations were performed using all 60 molecular orbitals.

FIG. 4.

The two occupied orbitals and five lowest-lying virtual orbitals for the H4 system in the trapezoidal (α = 0.001) [H4(a) system] and linear configuration [H4(b) system] shown in Fig. 3.

FIG. 4.

The two occupied orbitals and five lowest-lying virtual orbitals for the H4 system in the trapezoidal (α = 0.001) [H4(a) system] and linear configuration [H4(b) system] shown in Fig. 3.

Close modal
FIG. 5.

Orbital energies of lowest-lying molecular orbitals of H4(a) and H4(b) systems described by the cc-pVTZ basis set.

FIG. 5.

Orbital energies of lowest-lying molecular orbitals of H4(a) and H4(b) systems described by the cc-pVTZ basis set.

Close modal

In this paper, we discussed the excited-state extension of the DUCC formalism and its amenability for quantum computing. Using simple approximation schemes based on the utilization of the EOMCCSD wavefunction representation and lowest-rank contribution stemming from a single commutator electronic Hamiltonian and external σext(K) operator, we were able to demonstrate that the active-space representation of downfolded Hamiltonians can be used to reproduce a large portion of excited-state correlation effects. For H2 and H4(a) models, we observed significant improvements in excited-state energies compared to the diagonalization of the bare Hamiltonian in the active space, which was especially true for low-lying states dominated by active-space contributions (usually attributed to singly excited states and low-lying doubly excited states). As expected, for active-space dominated states characterized by higher excitation energies and more complicated configurational structure, the efficiency of simple approximation schemes discussed here deteriorates, which is indicative of the need for the inclusion of higher-order commutators and more efficient estimates of the σext(K) operators. This issue will be explored in future research by coupling DUCC-ex formalism with higher-rank EOMCC formulations. We have also demonstrated that quantum phase estimation can provide an efficient tool for testing various “excited-state” hypotheses for strongly correlated systems, which usually pose a significant challenge for existing many-body formalism due to the need of inclusion of higher rank excitation effects and high density of states in a narrow energy gap, which renders the numerical identification of state of interest numerically unfeasible. Instead, using PE techniques, one can define hypothesis or trial state and obtain in the course of calculation statistical footprint of all states that have non-negligible overlap with the hypothesis state. Although this fact is a well-known foundation of quantum computing, it deserves broader exposure. The stochastic characteristics of PE may be very useful in studies of excited-state processes, especially in the strongly correlated or metallic-like system.

This work was supported by the “Embedding Quantum Computing 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. A portion of this research was funded by the Quantum Algorithms, Software, and Architectures (QUASAR) Initiative at the Pacific Northwest National Laboratory (PNNL). It was conducted under the Laboratory Directed Research and Development Program at PNNL. All calculations have been performed using the Molecular Science Computing Facility (MSCF) in the Environmental Molecular Sciences Laboratory (EMSL) at the Pacific Northwest National Laboratory (PNNL). EMSL is funded by the Office of Biological and Environmental Research in the U.S. Department of Energy. PNNL is operated for the U.S. Department of Energy by the Battelle Memorial Institute under Contract No. DE-AC06-76RLO-1830.

1.
G.
Ortiz
,
J. E.
Gubernatis
,
E.
Knill
, and
R.
Laflamme
,
Phys. Rev. A
64
,
022319
(
2001
).
2.
A.
Aspuru-Guzik
,
A. D.
Dutoi
,
P. J.
Love
, and
M.
Head-Gordon
,
Science
309
,
1704
(
2005
).
3.
B. P.
Lanyon
,
J. D.
Whitfield
,
G. G.
Gillett
,
M. E.
Goggin
,
M. P.
Almeida
,
I.
Kassal
,
J. D.
Biamonte
,
M.
Mohseni
,
B. J.
Powell
,
M.
Barbieri
 et al,
Nat. Chem.
2
,
106
(
2010
).
4.
J. T.
Seeley
,
M. J.
Richard
, and
P. J.
Love
,
J. Chem. Phys.
137
,
224109
(
2012
).
5.
A.
Peruzzo
,
J.
McClean
,
P.
Shadbolt
,
M.-H.
Yung
,
X.-Q.
Zhou
,
P. J.
Love
,
A.
Aspuru-Guzik
, and
J. L.
O’brien
,
Nat. Commun.
5
,
4213
(
2014
).
6.
D.
Wecker
,
M. B.
Hastings
, and
M.
Troyer
,
Phys. Rev. A
92
,
042303
(
2015
).
7.
R.
Babbush
,
D. W.
Berry
,
I. D.
Kivlichan
,
A. Y.
Wei
,
P. J.
Love
, and
A.
Aspuru-Guzik
,
New J. Phys.
18
,
033032
(
2016
).
8.
J. R.
McClean
,
J.
Romero
,
R.
Babbush
, and
A.
Aspuru-Guzik
,
New J. Phys.
18
,
023023
(
2016
).
9.
V.
Havlíček
,
M.
Troyer
, and
J. D.
Whitfield
,
Phys. Rev. A
95
,
032332
(
2017
).
10.
J. R.
Fontalvo
,
R.
Babbush
,
J.
McClean
,
C.
Hempel
,
P. J.
Love
, and
A.
Aspuru-Guzik
, preprint arXiv:1701.02691 (
2017
).
11.
J. R.
McClean
,
I. D.
Kivlichan
,
D. S.
Steiger
,
Y.
Cao
,
E. S.
Fried
,
C.
Gidney
,
T.
Häner
,
V.
Havlíček
,
Z.
Jiang
,
M.
Neeley
 et al, preprint arXiv:1710.07629 (
2017
).
12.
Y.
Shen
,
X.
Zhang
,
S.
Zhang
,
J.-N.
Zhang
,
M.-H.
Yung
, and
K.
Kim
,
Phys. Rev. A
95
,
020501
(
2017
).
13.
M.
Reiher
,
N.
Wiebe
,
K. M.
Svore
,
D.
Wecker
, and
M.
Troyer
,
Proc. Natl. Acad. Sci. U. S. A.
114
,
7555
(
2017
).
14.
K.
Setia
and
J. D.
Whitfield
,
J. Chem. Phys.
148
,
164104
(
2018
).
15.
R.
Babbush
,
N.
Wiebe
,
J.
McClean
,
J.
McClain
,
H.
Neven
, and
G. K.-L.
Chan
,
Phys. Rev. X
8
,
011044
(
2018
).
16.
M.
Motta
,
E.
Ye
,
J. R.
McClean
,
Z.
Li
,
A. J.
Minnich
,
R.
Babbush
, and
G. K.
Chan
, preprint arXiv:1808.02625 (
2018
).
17.
H. R.
Grimsley
,
S. E.
Economou
,
E.
Barnes
, and
N. J.
Mayhall
,
Nat. Commun.
10
,
3007
(
2019
).
18.
K.
Andersson
,
P.-Å.
Malmqvist
,
B. O.
Roos
,
A. J.
Sadlej
, and
K.
Wolinski
,
J. Phys. Chem.
94
,
5483
(
1990
).
19.
K.
Andersson
,
P.-Å.
Malmqvist
, and
B. O.
Roos
,
J. Chem. Phys.
96
,
1218
(
1992
).
20.
C.
Angeli
,
R.
Cimiraglia
, and
J.-P.
Malrieu
,
Chem. Phys. Lett.
350
,
297
(
2001
).
21.
C.
Angeli
,
R.
Cimiraglia
,
S.
Evangelisti
,
T.
Leininger
, and
J.-P.
Malrieu
,
J. Chem. Phys.
114
,
10252
(
2001
).
22.
H.-J.
Werner
and
P. J.
Knowles
,
J. Chem. Phys.
89
,
5803
(
1988
).
23.
R. J.
Bartlett
and
M.
Musiał
,
Rev. Mod. Phys.
79
,
291
(
2007
).
24.
J.
Geertsen
,
M.
Rittby
, and
R. J.
Bartlett
,
Chem. Phys. Lett.
164
,
57
(
1989
).
25.
D. C.
Comeau
and
R. J.
Bartlett
,
Chem. Phys. Lett.
207
,
414
(
1993
).
26.
J. F.
Stanton
and
R. J.
Bartlett
,
J. Chem. Phys.
98
,
7029
(
1993
).
27.
P.
Piecuch
and
R. J.
Bartlett
,
Advances in Quantum Chemistry
(
Elsevier
,
1999
), Vol. 34, pp.
295
380
.
28.
A. I.
Krylov
,
Chem. Phys. Lett.
338
,
375
(
2001
).
29.
S.
Hirata
,
J. Chem. Phys.
121
,
51
(
2004
).
30.
D.
Mukherjee
,
R. K.
Moitra
, and
A.
Mukhopadhyay
,
Mol. Phys.
30
,
1861
(
1975
).
31.
U.
Kaldor
,
Theor. Chim. Acta
80
,
427
(
1991
).
32.
C.
Rittby
and
R.
Bartlett
,
Theor. Chim. Acta
80
,
469
(
1991
).
33.
L.
Meissner
,
J. Chem. Phys.
108
,
9227
(
1998
).
34.
M.
Musiał
and
R. J.
Bartlett
,
J. Chem. Phys.
129
,
134105
(
2008
).
35.
B.
Jeziorski
and
H. J.
Monkhorst
,
Phys. Rev. A
24
,
1668
(
1981
).
36.
L.
Meissner
,
K.
Jankowski
, and
J.
Wasilewski
,
Int. J. Quantum Chem.
34
,
535
(
1988
).
37.
J.
Paldus
,
P.
Piecuch
,
L.
Pylypow
, and
B.
Jeziorski
,
Phys. Rev. A
47
,
2738
(
1993
).
38.
P.
Piecuch
and
J.
Paldus
,
Phys. Rev. A
49
,
3479
(
1994
).
39.
P.
Piecuch
and
J.
Paldus
,
Theor. Chem. Acc.
83
,
69
(
1992
).
40.
U. S.
Mahapatra
,
B.
Datta
, and
D.
Mukherjee
,
Mol. Phys.
94
,
157
(
1998
).
41.
U. S.
Mahapatra
,
B.
Datta
,
B.
Bandyopadhyay
, and
D.
Mukherjee
,
Adv. Quantum Chem.
30
,
163
(
1998
).
42.
F. A.
Evangelista
,
W. D.
Allen
, and
H. F.
Schaefer
 III
,
J. Chem. Phys.
127
,
024102
(
2007
).
43.
J.
Pittner
,
J. Chem. Phys.
118
,
10876
(
2003
).
44.
M.
Hanauer
and
A.
Köhn
,
J. Chem. Phys.
134
,
204111
(
2011
).
45.
Y. A.
Aoto
and
A.
Köhn
,
J. Chem. Phys.
144
,
074103
(
2016
).
46.
A.
Köhn
,
M.
Hanauer
,
L. A.
Mück
,
T.-C.
Jagau
, and
J.
Gauss
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
3
,
176
(
2013
).
47.
D. I.
Lyakh
,
M.
Musiał
,
V. F.
Lotrich
, and
R. J.
Bartlett
,
Chem. Rev.
112
,
182
(
2012
).
48.
49.
U.
Schollwöck
,
Rev. Mod. Phys.
77
,
259
(
2005
).
50.
O.
Legeza
and
J.
Sólyom
,
Phys. Rev. B
68
,
195116
(
2003
).
51.
G. K.-L.
Chan
and
S.
Sharma
,
Annu. Rev. Phys. Chem.
62
,
465
(
2011
).
52.
J. D.
Watts
and
R. J.
Bartlett
,
Spectrochim. Acta, Part A
55
,
495
(
1999
).
53.
K.
Kowalski
and
P.
Piecuch
,
J. Chem. Phys.
115
,
643
(
2001
).
54.
S. A.
Kucharski
,
M.
Włoch
,
M.
Musiał
, and
R. J.
Bartlett
,
J. Chem. Phys.
115
,
8263
(
2001
).
55.
K.
Bhaskaran-Nair
and
K.
Kowalski
,
J. Chem. Phys.
137
,
216101
(
2012
).
56.
G. H.
Booth
,
A. J.
Thom
, and
A.
Alavi
,
J. Chem. Phys.
131
,
054106
(
2009
).
57.
D.
Cleland
,
G. H.
Booth
, and
A.
Alavi
,
J. Chem. Phys.
132
,
041103
(
2010
).
58.
J. E.
Deustua
,
S. H.
Yuwono
,
J.
Shen
, and
P.
Piecuch
,
J. Chem. Phys.
150
,
111101
(
2019
).
59.
J.
Romero
,
R.
Babbush
,
J. R.
McClean
,
C.
Hempel
,
P. J.
Love
, and
A.
Aspuru-Guzik
,
Quantum Sci. Technol.
4
,
014008
(
2018
).
60.
A.
Kandala
,
K.
Temme
,
A. D.
Corcoles
,
A.
Mezzacapo
,
J. M.
Chow
, and
J. M.
Gambetta
, preprint arXiv:1805.04492 (
2018
).
61.
J. I.
Colless
,
V. V.
Ramasesh
,
D.
Dahlen
,
M. S.
Blok
,
M. E.
Kimchi-Schwartz
,
J. R.
McClean
,
J.
Carter
,
W. A.
de Jong
, and
I.
Siddiqi
,
Phys. Rev. X
8
,
011021
(
2018
).
62.
Y.
Cao
,
J.
Romero
,
J. P.
Olson
,
M.
Degroote
,
P. D.
Johnson
,
M.
Kieferová
,
I. D.
Kivlichan
,
T.
Menke
,
B.
Peropadre
,
N. P.
Sawaya
 et al,
Chem. Rev.
119
,
10856
(
2019
).
63.
P. J.
Ollitrault
,
A.
Kandala
,
C.-F.
Chen
,
P. K.
Barkoutsos
,
A.
Mezzacapo
,
M.
Pistoia
,
S.
Sheldon
,
S.
Woerner
,
J.
Gambetta
, and
I.
Tavernelli
, preprint arXiv:1910.12890 (
2019
).
64.
A.
Luis
and
J.
Peřina
,
Phys. Rev. A
54
,
4564
(
1996
).
65.
R.
Cleve
,
A.
Ekert
,
C.
Macchiavello
, and
M.
Mosca
,
Proc. R. Soc. London, Ser. A
454
,
339
(
1998
).
66.
D. W.
Berry
,
G.
Ahokas
,
R.
Cleve
, and
B. C.
Sanders
,
Commun. Math. Phys.
270
,
359
(
2007
).
67.
A. M.
Childs
,
Commun. Math. Phys.
294
,
581
(
2010
).
68.
S.
Kimmel
,
G. H.
Low
, and
T. J.
Yoder
,
Phys. Rev. A
92
,
062315
(
2015
).
69.
N.
Wiebe
and
C.
Granade
,
Phys. Rev. Lett.
117
,
010503
(
2016
).
70.
T.
Häner
,
D. S.
Steiger
,
M.
Smelyanskiy
, and
M.
Troyer
, in
SC’16: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis
(
IEEE
,
2016
), pp.
866
874
.
71.
D.
Poulin
,
A.
Kitaev
,
D. S.
Steiger
,
M. B.
Hastings
, and
M.
Troyer
,
Phys. Rev. Lett.
121
,
010501
(
2018
).
72.
H.
Wang
,
S.
Kais
,
A.
Aspuru-Guzik
, and
M. R.
Hoffmann
,
Phys. Chem. Chem. Phys.
10
,
5388
(
2008
).
73.
J. R.
McClean
,
M. E.
Kimchi-Schwartz
,
J.
Carter
, and
W. A.
de Jong
,
Phys. Rev. A
95
,
042308
(
2017
).
74.
G. H.
Low
,
N. P.
Bauman
,
C. E.
Granade
,
B.
Peng
,
N.
Wiebe
,
E. J.
Bylaska
,
D.
Wecker
,
S.
Krishnamoorthy
,
M.
Roetteler
,
K.
Kowalski
 et al, preprint arXiv:1904.01131 (
2019
).
75.
O.
Higgott
,
D.
Wang
, and
S.
Brierley
,
Quantum
3
,
156
(
2019
).
76.
G. H.
Low
and
I. L.
Chuang
,
Quantum
3
,
163
(
2019
).
77.
K.
Jankowski
and
J.
Paldus
,
Int. J. Quantum Chem.
18
,
1243
(
1980
).
78.
T.
Minami
and
M.
Nakano
,
J. Phys. Chem. Lett.
3
,
145
(
2011
).
79.
N. P.
Bauman
,
E. J.
Bylaska
,
S.
Krishnamoorthy
,
G. H.
Low
,
N.
Wiebe
,
C. E.
Granade
,
M.
Roetteler
,
M.
Troyer
, and
K.
Kowalski
,
J. Chem. Phys.
151
,
014107
(
2019
).
80.
K.
Kowalski
,
J. Chem. Phys.
148
,
094104
(
2018
).
81.
P.
Piecuch
,
N.
Oliphant
, and
L.
Adamowicz
,
J. Chem. Phys.
99
,
1875
(
1993
).
83.
N.
Oliphant
and
L.
Adamowicz
,
J. Chem. Phys.
94
,
1229
(
1991
).
84.
N.
Oliphant
and
L.
Adamowicz
,
J. Chem. Phys.
96
,
3739
(
1992
).
85.
S.
Hirata
,
J. Phys. Chem. A
107
,
9887
(
2003
).
86.
M.
Valiev
,
E.
Bylaska
,
N.
Govind
,
K.
Kowalski
,
T.
Straatsma
,
H. V.
Dam
,
D.
Wang
,
J.
Nieplocha
,
E.
Apra
,
T.
Windus
, and
W.
de Jong
,
Comput. Phys. Commun.
181
,
1477
(
2010
).
87.
M. A.
Nielsen
and
I. L.
Chuang
,
Quantum Computation and Quantum Information: 10th Anniversary Edition
, 10th ed. (
Cambridge University Press
,
New York, NY, USA
,
2011
).
88.
G. H.
Low
and
I. L.
Chuang
,
Phys. Rev. Lett.
118
,
010501
(
2017
).
89.
D. W.
Berry
,
M.
Kieferová
,
A.
Scherer
,
Y. R.
Sanders
,
G. H.
Low
,
N.
Wiebe
,
C.
Gidney
, and
R.
Babbush
,
npj Quantum Inf.
4
,
22
(
2018
).
90.
D.
Poulin
,
A.
Kitaev
,
D. S.
Steiger
,
M. B.
Hastings
, and
M.
Troyer
, preprint arXiv:1711.11025 (
2017
).
91.
T. H.
Dunning
, Jr.
,
J. Chem. Phys.
90
,
1007
(
1989
).