Construction of Green's functions on a quantum computer: applications to molecular systems

TL;DR

This paper presents a quantum computing scheme for constructing one-particle Green's functions of molecular systems, enabling spectral analysis and energy calculations with probabilistic state preparation and quantum circuits.

Contribution

It introduces a novel quantum algorithm for Green's function construction using minimal ancillary qubits and demonstrates its effectiveness through molecular simulations.

Findings

Successfully simulated GFs for LiH and H2O molecules.

Achieved accurate spectral results compared to exact UCC and FCI calculations.

Validated sampling method via Galitskii--Migdal energy estimates.

Abstract

We propose a scheme for the construction of one-particle Green's function (GF) of an interacting electronic system via statistical sampling on a quantum computer. Although the non-unitarity of creation and annihilation operators for the electronic spin orbitals prevents us from preparing specific states selectively, probabilistic state preparation is demonstrated to be possible for the qubits. We provide quantum circuits equipped with at most two ancillary qubits for obtaining all the components of GF. We perform simulations of such construction of GFs for LiH and H$_2$O molecules based on the unitary coupled-cluster (UCC) method to demonstrate the validity of our scheme by comparing the spectra exact within UCC and those from full configuration interaction calculations. We also examine the accuracy of sampling method by exploiting the Galitskii--Migdal formula, which gives the total energy only from the GF.

Equations

$$U(\mathbf{\theta })=\exp [T(\mathbf{\theta })-T(\mathbf{\theta }{)}^{†}],$$

$$G_{m{m}^{\prime }}(z)=G_{m{m}^{\prime }}^{(\mathrm{e})}(z)+G_{m{m}^{\prime }}^{(\mathrm{h})}(z)$$

$$G_{m{m}^{\prime }}^{(\mathrm{e})}(z)=⟨{\Psi }_{\mathrm{gs}}^N|a_m\frac{1}{z+E_{\mathrm{gs}}^N-\mathcal{H}}{a}_{m^{\prime }}^{†}|{\Psi }_{\mathrm{gs}}^N⟩$$

$$=\sum _{\lambda \in N+1}\frac{B_{\lambda m{m}^{\prime }}^{(\mathrm{e})}}{z+E_{\mathrm{gs}}^N-E_{\lambda }^{N+1}}$$

$$G_{m{m}^{\prime }}^{(\mathrm{h})}(z)=⟨{\Psi }_{\mathrm{gs}}^N|a_{m^{\prime }}^{†}\frac{1}{z-E_{\mathrm{gs}}^N+\mathcal{H}}{a}_m|{\Psi }_{\mathrm{gs}}^N⟩$$

$$=\sum _{\lambda \in N-1}\frac{B_{\lambda m{m}^{\prime }}^{(\mathrm{h})}}{z+E_{\lambda }^{N-1}-E_{\mathrm{gs}}^N}$$

$$B_{\lambda m{m}^{\prime }}^{(\mathrm{e})}\equiv ⟨{\Psi }_{\mathrm{gs}}^N|a_m|{\Psi }_{\lambda }^{N+1}⟩⟨{\Psi }_{\lambda }^{N+1}|a_{m^{\prime }}^{†}|{\Psi }_{\mathrm{gs}}^N⟩$$

$$B_{\lambda m{m}^{\prime }}^{(\mathrm{h})}\equiv ⟨{\Psi }_{\mathrm{gs}}^N|a_{m^{\prime }}^{†}|{\Psi }_{\lambda }^{N-1}⟩⟨{\Psi }_{\lambda }^{N-1}|a_m|{\Psi }_{\mathrm{gs}}^N⟩$$

$$A(\omega )=-\frac{1}{\pi }\mathrm{ImTr}G(\omega +i\delta )$$

$$a_m^{†}=\frac{U_{0m}-U_{1m}}{2}$$

$$a_m=\frac{U_{0m}+U_{1m}}{2}$$

$$|0⟩\otimes |\psi ⟩⟼|0⟩\otimes \frac{U_{0m}+U_{1m}}{2}|\psi ⟩+|1⟩\otimes \frac{U_{0m}-U_{1m}}{2}|\psi ⟩$$

$$=|0⟩\otimes a_m|\psi ⟩+|1⟩\otimes a_m^{†}|\psi ⟩\equiv |{\Phi }_m⟩.$$

$$|{\Phi }_m⟩\stackrel{|0⟩\mathrm{observed}}{⟼}|0⟩\otimes \frac{a_m}{\sqrt{p_m(\mathrm{h})}}|\psi ⟩$$

$$prob.⟨\psi |a_m^{†}{a}_m|\psi ⟩\equiv p_m(\mathrm{h})$$

$$|{\Phi }_m⟩\stackrel{|1⟩\mathrm{observed}}{⟼}|1⟩\otimes \frac{a_m^{†}}{\sqrt{p_m(\mathrm{e})}}|\psi ⟩$$

$$prob.⟨\psi |a_m{a}_m^{†}|\psi ⟩\equiv p_m(\mathrm{e})$$

$$a_{m{m}^{\prime }}^{\pm }\equiv \frac{a_m\pm e^{-i\pi /4}{a}_{m^{\prime }}}{2}$$

$$a_{m{m}^{\prime }}^{\pm †}\equiv \frac{a_m^{†}\pm e^{i\pi /4}{a}_{m^{\prime }}^{†}}{2},$$

$$|{\Psi }_{m{m}^{\prime }}^{(\mathrm{e})\pm }⟩\equiv a_{m{m}^{\prime }}^{\pm †}|{\Psi }_{\mathrm{gs}}^N⟩$$

$$D_{\lambda m{m}^{\prime }}^{(\mathrm{e})\pm }\equiv |⟨{\Psi }_{\lambda }^{N+1}|{\Psi }_{m{m}^{\prime }}^{(\mathrm{e})\pm }⟩{|}^2$$

$$=\frac{B_{\lambda mm}^{(\mathrm{e})}+B_{\lambda m^{\prime }{m}^{\prime }}^{(\mathrm{e})}}{4}\pm \frac{e^{i\pi /4}{B}_{\lambda m{m}^{\prime }}^{(\mathrm{e})}+e^{-i\pi /4}{B}_{\lambda m^{\prime }m}^{(\mathrm{e})}}{4}.$$

$$\displaystyle B_{\lambda mm^{\prime}}^{\mathrm{(e)}}=\pm\Bigg{(}2e^{-i\pi/4}D_{\lambda mm^{\prime}}^{\mathrm{(e)}\pm}+2e^{i\pi/4}D_{\lambda m^{\prime}m}^{\mathrm{(e)}\pm}$$

$$\displaystyle-\frac{B_{\lambda mm}^{\mathrm{(e)}}+B_{\lambda m^{\prime}m^{\prime}}^{\mathrm{(e)}}}{\sqrt{2}}\Bigg{)}.$$

$$B_{\lambda m{m}^{\prime }}^{(\mathrm{e})}=e^{-i\pi /4}(D_{\lambda m{m}^{\prime }}^{(\mathrm{e})+}-D_{\lambda m{m}^{\prime }}^{(\mathrm{e})-})$$

$$+e^{i\pi /4}(D_{\lambda m^{\prime }m}^{(\mathrm{e})+}-D_{\lambda m^{\prime }m}^{(\mathrm{e})-}).$$

$$|{\Psi }_{m{m}^{\prime }}^{(\mathrm{h})\pm }⟩\equiv a_{m{m}^{\prime }}^{\pm }|{\Psi }_{\mathrm{gs}}^N⟩,$$

$$|q_1^{\mathrm{A}}=0⟩\otimes |q_0^{\mathrm{A}}=0⟩\otimes |\psi ⟩$$

$$⟼|0⟩\otimes |0⟩\otimes \frac{U_{0m}+U_{1m}+e^{i\pi /4}(U_{0m^{\prime }}+U_{1m^{\prime }})}{4}|\psi ⟩$$

$$+|0⟩\otimes |1⟩\otimes \frac{U_{0m}-U_{1m}+e^{i\pi /4}(U_{0m^{\prime }}-U_{1m^{\prime }})}{4}|\psi ⟩$$

$$+|1⟩\otimes |0⟩\otimes \frac{U_{0m}+U_{1m}-e^{i\pi /4}(U_{0m^{\prime }}+U_{1m^{\prime }})}{4}|\psi ⟩$$

$$+|1⟩\otimes |1⟩\otimes \frac{U_{0m}-U_{1m}-e^{i\pi /4}(U_{0m^{\prime }}-U_{1m^{\prime }})}{4}|\psi ⟩$$

$$=|0⟩\otimes |0⟩\otimes e^{i\pi /4}{a}_{m^{\prime }m}^{+}|\psi ⟩+|0⟩\otimes |1⟩\otimes a_{m{m}^{\prime }}^{+†}|\psi ⟩$$

$$-|1⟩\otimes |0⟩\otimes e^{i\pi /4}{a}_{m^{\prime }m}^{-}|\psi ⟩+|1⟩\otimes |1⟩\otimes a_{m{m}^{\prime }}^{-†}|\psi ⟩$$

$$\equiv |{\Phi }_{m{m}^{\prime }}⟩.$$

$$|{\Phi }_{m{m}^{\prime }}⟩\stackrel{|0⟩\otimes |0⟩\mathrm{observed}}{⟼}|0⟩\otimes |0⟩\otimes \frac{a_{m^{\prime }m}^{+}}{\sqrt{p_{m{m}^{\prime }}(\mathrm{h},+)}}|\psi ⟩$$

$$prob.⟨\psi |a_{m^{\prime }m}^{+†}{a}_{m^{\prime }m}^{+}|\psi ⟩\equiv p_{m{m}^{\prime }}(\mathrm{h},+)$$

$$|{\Phi }_{m{m}^{\prime }}⟩\stackrel{|0⟩\otimes |1⟩\mathrm{observed}}{⟼}|0⟩\otimes |1⟩\otimes \frac{a_{m{m}^{\prime }}^{+†}}{\sqrt{p_{m{m}^{\prime }}(\mathrm{e},+)}}|\psi ⟩$$

$$prob.⟨\psi |a_{m{m}^{\prime }}^{+}{a}_{m{m}^{\prime }}^{+†}|\psi ⟩\equiv p_{m{m}^{\prime }}(\mathrm{e},+)$$

$$|{\Phi }_{m{m}^{\prime }}⟩\stackrel{|1⟩\otimes |0⟩\mathrm{observed}}{⟼}|1⟩\otimes |0⟩\otimes \frac{a_{m^{\prime }m}^{-}}{\sqrt{p_{m{m}^{\prime }}(\mathrm{h},-)}}|\psi ⟩$$

Full text

††preprint: APS/123-QED

I Introduction

Quantum chemistry calculationsMcArdle et al. (2018) as a kind of quantum simulationFeynman (1982) have been drawing attention increasingly since they serve as lucid proof-of-principle for quantum computationNielsen and Chuang (2011) and, at the same time, are directly related to the state-of-the-art quantum hardware. The electronic states and the operators acting on them for a given Hamiltonian have to be mapped to the qubits comprising a quantum computer and the operators for them by an appropriate transformation.Seeley et al. (2012) The Jordan–Wigner (JW)Jordan and Wigner (1928) and Bravyi–Kitaev (BK)Bravyi and Kitaev (2002) transformations are often used for quantum chemistry calculations. Various approaches for obtaining the energy spectra of a many-electron system have been proposed. The earliest oneAspuru-Guzik et al. (2005) employs the quantum phase estimation (QPE) algorithmAbrams and Lloyd (1997, 1999) and the Suzuki–Trotter decompositionSuzuki (1992) of the qubit Hamiltonian into a sequence of one- and two-qubit logic gatesWhitfield et al. (2011) for unitary operations. This approach was realizedO’Malley et al. (2016) by using superconducting qubits. Variational quantum eigensolver (VQE) is a newer approach, in which a trial many-electron state is prepared via a quantum circuit with parameters to be optimized aiming at the ground state. It uses a classical computer for updating the parameters based on the measurement results of the qubits, which is why it is also called a quantum-classical hybrid algorithm.McClean et al. (2016) This approach was first realizedPeruzzo et al. (2014) by using a quantum photonic device. It has also been realized by superconductingO’Malley et al. (2016); Kandala et al. (2017) and ion trapHempel et al. (2018) quantum computers. Another approach for obtaining the energy spectra is the imaginary-time evolution. It was recently proposedJones et al. (2019); McArdle et al. (2019); Yuan et al. (2019) as a quantum-classical hybrid algorithm based on the McLachlan’s variational principle.McLachlan (1964)

An experiment of photoelectron spectroscopy (PES) irradiates light to a sample and measures the energy of photoelectrons coming out of the sample. An experiment of inverse PES is for the reverse process of PES. Particularly for angle-resolved photoemission spectroscopy (ARPES), the measured spectra of an interacting electronic system are often explained via the one-particle Green’s function (GF).Damascelli (2004); Moser (2017); Kosugi et al. (2017) Since the GF contains rich information about the correlation effects in an electronic systemKosugi and Matsushita (2018), the GFs in the context of quantum chemistryNooijen and Snijders (1992, 1993, 1995); Kowalski et al. (2014); Bhaskaran-Nair et al. (2016) (on classical computers) have been intensively studied recently for isolatedKosugi et al. (2018); Nishi et al. (2018); Peng and Kowalski (2018a); Peng et al. (2019); Peng and Kowalski (2018b) and periodicFurukawa et al. (2018); Kosugi and Matsushita (2019) systems. The reliable calculation of GFs is hence as important as that of the ground-state energies for molecular and solid-state systems. It is, however, essentially expensive for classical computation since it demands large memory and, often simultaneously, large storage for description of an electronic state made up of lots of Slater determinants. A quantum computer allows for, on the other hand, representation of such an electronic state using the qubits thanks to the principle of superposition. It is thus worth developing tools for electronic-structure calculations on quantum computers, which are coming to practical usage.

In this study, we propose a scheme for the construction of one-particle GF of an interacting electronic system via statistical sampling on a quantum computer. We introduce quantum circuits for probabilistic state preparation which allow us to calculate the GF from the histogram obtained via measurements on the qubits. Our scheme exploits the probabilistic preparation of electron-added and -removed states, in contrast to the existing methods for GFs.Wecker et al. (2015); Bauer et al. (2016) For demonstrating the validity of our scheme, we perform simulations of such construction of molecular GFs based on the unitary coupled-cluster (UCC) methodRomero et al. (2018) by referring to the spectral functions exact within UCC. We also examine the accuracy of sampling method by calculating the correlation energies from the GFs.

This paper is organized as follows. In Section II, we explain the theoretical perspective of our scheme. In particular, we describe the quantum circuits in detail for obtaining GFs via statistical sampling. In Section III, we describe the computational details for our simulations on a classical computer. In Section IV, we show the simulation results of quantum computation for LiH and H2O molecules. In Section V, we provide the conclusions.

II Methods

We describe the scheme for constructing the GF using a quantum computer below in detail. Although we use the ground states obtained in UCC calculations in the present study, the scheme is applicable as long as the ground state can be prepared as the qubits.

II.1 Unitary coupled-cluster method

Let us consider an interacting $N$-electron system whose second-quantized Hamiltonian is $\mathcal{H}$. A VQE calculationPeruzzo et al. (2014) of quantum chemistry using the UCC methodRomero et al. (2018) starts from an ansatz of the form

[TABLE]

where $T(\bm{\theta})$ is an appropriately chosen cluster operator that depends on parameter(s) $\bm{\theta}$. The transformation $U(\bm{\theta})$, which is unitary by definition, is used to construct a trial ground state $|\Psi(\bm{\theta})\rangle\equiv U(\bm{\theta})|\Psi_{\mathrm{ref}}\rangle$ for given $\bm{\theta}$ from a reference state $|\Psi_{\mathrm{ref}}\rangle$. In a practical VQE process, the unitary transformation is implemented as parametrized operations on the qubits comprising a quantum computer. The expected total energy $E(\bm{\theta})=\langle\Psi(\bm{\theta})|\mathcal{H}|\Psi(\bm{\theta})\rangle$ is obtained via measurements (Hamiltonian averagingO’Malley et al. (2016); Hempel et al. (2018); Peruzzo et al. (2014)), which is then used to update $\bm{\theta}$ iteratively according to an optimization scheme on a classical computer so that the measured energy at the next iteration is lower. The optimized trial state will be used as the ground state $|\Psi_{\mathrm{gs}}^{N}\rangle$ for our scheme described below.

Although we introduce the electronic cluster operators for ansätze and rewrite them into the qubit representation in the present study, one can instead start directly from ansätze given as qubit operators. The qubit coupled-clusterRyabinkin et al. (2018) is an approach in this direction.

II.2 One-particle GFs

II.2.1 Definition

Although we assume the ground state $|\Psi_{\mathrm{gs}}^{N}\rangle$ to be non-degenerate and to be at zero temperature for simplicity, the expressions provided below are easily extended for systems having degenerate ground states at nonzero temperature. The one-particle GFFetter and Walecka (2003); Stefanucci and van Leeuwen (2013) of the system in frequency domain is given by

[TABLE]

for a complex frequency $z$, where

[TABLE]

and

[TABLE]

are the electron- and hole-excitation parts of the GF, respectively. $a^{\dagger}_{m}$ and $a_{m}$ are the creation and annihilation operators, respectively, of an electron at the $m$th spin orbital. $E_{\mathrm{gs}}^{N}$ is the ground-state energy and $E_{\lambda}^{N\pm 1}$ is the $\lambda$th energy eigenvalue of the $(N\pm 1)$-electron states.

[TABLE]

and

[TABLE]

are the transition matrix elements. The spectral function is defined via the GF as

[TABLE]

for a real $\omega$ with a small positive constant $\delta$ for ensuring causality.

It is clear from eqs. (3) and (4) that the calculation of GF requires not only the many-electron energy eigenvalues but also the transition matrix elements. Various approaches for obtaining many-electron energy eigenvalues on a quantum computer have been proposedHiggott et al. (2019); Jones et al. (2019); McClean et al. (2017); Santagati et al. (2018); Colless et al. (2018) and we can choose any alternative from them by comparing their precision and restriction from the viewpoints of algorithm and hardware. As for the transition matrix elements, however, there exists no established way for calculation of them on a quantum computer to our knowledge. We therefore propose a scheme for the construction of GF via statistical sampling and describe it below in detail. Our protocol is designed for obtaining the numerators on the RHSs in eqs. (3) and (4), provided that the denominators have been known.

II.2.2 Circuits for diagonal components

In a typical scheme for the construction of GFs on a classical computerNooijen and Snijders (1992, 1993, 1995); Kowalski et al. (2014); Bhaskaran-Nair et al. (2016), the equation-of-motion coupled-cluster (EOM-CC) approach is adopted to obtain the energy eigenvalues and the transition matrix elements for the $(N\pm 1)$-electron intermediate states. In the present case, one might think by looking at eq. (5) that $B_{\lambda mm^{\prime}}^{\mathrm{(e)}}$ can be easily calculated by preparing the qubit representations of $|\Psi_{\lambda}^{N+1}\rangle$ and $|\Psi_{m}^{\mathrm{(e)}}\rangle\equiv a_{m}^{\dagger}|\Psi_{\mathrm{gs}}^{N}\rangle,$ between which the inner product is calculated using the swap test or its versionsCincio et al. (2018); Higgott et al. (2019); Garcia-Escartin and Chamorro-Posada (2013) with phase factors. Such an approach is, however, difficult in fact. It is because the creation operator is not unitary and the norm of the electron-added state is not conserved in general, that is, $\langle\Psi_{m}^{\mathrm{(e)}}|\Psi_{m}^{\mathrm{(e)}}\rangle\neq 1$. This fact prevents one from preparing a specific electron-added state selectively since a quantum circuit can apply only unitary operations to qubits. This difficulty is similarly the case for the electron-removed (hole-added) state $|\Psi_{m}^{\mathrm{(h)}}\rangle\equiv a_{m}|\Psi_{\mathrm{gs}}^{N}\rangle$ . To circumvent this difficulty, we have to resort to another approach.

As explained in Introduction, the JWJordan and Wigner (1928) and BKBravyi and Kitaev (2002) transformations are often used for mapping a many-electron state to a many-qubit state. We do not distinguish between the kets as many-electron states and those as many-qubit states in what follows since such simplification will not cause confusion for the readers. By looking at the definitions of the transformations [see, e.g., eqs. (34), (39), and (40) in Ref.Seeley et al. (2012)], we can notice that for both transformations any pair of electronic creation and annihilation operators can be expressed by using two unitary operators $U_{0m}$ and $U_{1m}$ on qubits as

[TABLE]

and

[TABLE]

for a given $m$, regardless of the number of qubits comprising the quantum computer. Such a decomposition of electronic operators is in fact always possible since $a_{m}^{\dagger}+a_{m}$ and $a_{m}^{\dagger}-a_{m}$ are ensured to be unitary by the anti-commutation relation, reminding us of the Majorana fermions.Elliott and Franz (2015) We introduce a trick for state preparation by exploiting this fact. Specifically, we construct a circuit $\mathcal{C}_{m}$ equipped with an ancillary qubit $|q^{\mathrm{A}}\rangle$ by implementing the controlled operations of $U_{0m}$ and $U_{1m}$, as shown in Figure 1. The whole system consists of the ancilla and an arbitrary input register $|\psi\rangle$, whose state changes by undergoing the circuit as

[TABLE]

The action of the circuit to the whole system is easily confirmed to be unitary due to the anti-commutation relation between the electronic operators. The projective measurementNielsen and Chuang (2011) on the ancillary bit is represented by the two operators $\mathcal{P}_{q}=|q\rangle\langle q|\otimes I\,(q=0,1)$, for which $|q\rangle$ is observed with a probability $\langle\Phi_{m}|\mathcal{P}_{q}|\Phi_{m}\rangle$. The state of the whole system collapses immediately after the measurement as follows:

[TABLE]

This result implies that $\mathcal{C}_{m}$ allows us to prepare the two states $a_{m}^{\dagger}|\psi\rangle$ and $a_{m}|\psi\rangle$ probabilistically apart from their normalization constants. For the number $N_{\mathrm{meas}}$ of measurements on the ancilla, the probability distribution of counted outcomes for $|0\rangle$, or equivalently $|1\rangle$, is a binomial distribution. The probability distribution thus converges to the normal distribution for many repeated measurements and the error of normalization constant scales as $N_{\mathrm{meas}}^{-1/2}$. This circuit is used for obtaining the diagonal components of the GF, as explained later.

II.2.3 Circuits for off-diagonal components

For the $m$th and $m^{\prime}$th spin orbitals ($m\neq m^{\prime}$), we define the following auxiliary creation and annihilation operators

[TABLE]

and

[TABLE]

respectively, which are the Hermitian conjugate of each other. Unnormalized auxiliary $(N+1)$-electron states

[TABLE]

can have overlaps with the energy eigenstates as

[TABLE]

By solving eq. (16) for the off-diagonal component of $B_{\lambda}^{\mathrm{(e)}}$, we can calculate it from $D_{\lambda}^{\mathrm{(e)}\pm}$ and the diagonal components of $B_{\lambda}^{\mathrm{(e)}}$ as

[TABLE]

From the two expressions for both signs in eq. (17), we can obtain the off-diagonal component of $B_{\lambda}^{\mathrm{(e)}}$ only from $D_{\lambda}^{\mathrm{(e)}\pm}$ as

[TABLE]

For unnormalized auxiliary $(N-1)$-electron states

[TABLE]

the expression of $D_{\lambda mm^{\prime}}^{\mathrm{(h)}\pm}\equiv|\langle\Psi_{\lambda}^{N-1}|\Psi_{mm^{\prime}}^{\mathrm{(h)}\pm}\rangle|^{2}$ is the same as that in eq. (16) with $(\mathrm{e})$ replaced by $(\mathrm{h})$. This means that we can calculate the off-diagonal component of $B_{\lambda}^{\mathrm{(h)}}$ from $D_{\lambda}^{\mathrm{(h)}\pm}$ by using the same expression as eq. (18) with the replacement.

We construct a circuit $\mathcal{C}_{mm^{\prime}}$ equipped with two ancillary qubits $|q_{0}^{\mathrm{A}}\rangle$ and $|q_{1}^{\mathrm{A}}\rangle$ by implementing the controlled operations of $U_{0m},U_{1m},U_{0m^{\prime}}$, and $U_{1m^{\prime}}$, as shown in Figure 2. The whole system consists of the ancillae and an arbitrary input register $|\psi\rangle$, whose state changes by undergoing the circuit as

[TABLE]

The action of the circuit to the whole system is easily confirmed to be unitary due to the anti-commutation relation between the electronic operators. The projective measurement on the ancillary bits is represented by the four operators $\mathcal{P}_{qq^{\prime}}=|q\rangle\langle q|\otimes|q^{\prime}\rangle\langle q^{\prime}|\otimes I\,(q,q^{\prime}=0,1)$, for which $|q\rangle\otimes|q^{\prime}\rangle$ is observed with a probability $\langle\Phi_{mm^{\prime}}|\mathcal{P}_{qq^{\prime}}|\Phi_{mm^{\prime}}\rangle$. The state of the whole system collapses immediately after the measurement as follows:

[TABLE]

This result implies that $\mathcal{C}_{mm^{\prime}}$ allows us to prepare the four states $a_{mm^{\prime}}^{\pm\dagger}|\psi\rangle$ and $a_{m^{\prime}m}^{\pm}|\psi\rangle$ probabilistically apart from their normalization constants. For the number $N_{\mathrm{meas}}$ of measurements on the ancillae, the error of normalization constant scales as $N_{\mathrm{meas}}^{-1/2}$ similarly to the case for diagonal components. This circuit is used for obtaining the off-diagonal components of the GF, as explained below.

II.2.4 Transition matrices via statistical sampling

Given the results of measurement on the ancillary bit(s), we have the register $|\widetilde{\psi}\rangle$ representing the $N_{e}$-electron state with $N_{e}=N+1$ or $N-1$. Then we perform QPE by inputting $|\widetilde{\psi}\rangle$ to obtain the energy eigenvalue in the subspace spanned by the $N_{e}$-electron states. A QPE experiment inevitably suffers from probabilistic errors that depend on the number and the initial states of qubits.Nielsen and Chuang (2011); Harrow et al. (2009) Furthermore, the results are affected by the number of steps for the Suzuki–Trotter decomposition and the order of partial Hamiltonians. We assume for simplicity, however, that the QPE procedure is realized on a quantum computer with ideal precision. We will thus find the estimated value to be $E_{\lambda}^{N_{e}}$ with a probability $|\langle\Psi_{\lambda}^{N_{e}}|\widetilde{\psi}\rangle|^{2}$. Nielsen and Chuang (2011)

If we input $|\Psi_{\mathrm{gs}}^{N}\rangle$ to the diagonal circuit $\mathcal{C}_{m}$ in Fig. 1 and process the whole system the way described above, the energy eigenvalue $E_{\lambda}^{N+1}$ will be obtained with a probability [see eq. (12)]

[TABLE]

while $E_{\lambda}^{N-1}$ will be obtained with a probability [see eq. (11)]

[TABLE]

This means that we can get the diagonal components of transition matrices $B_{\lambda}^{(\mathrm{e})}$ and $B_{\lambda}^{(\mathrm{h})}$ via statistical sampling for a fixed $m$. It is easily confirmed that $\sum_{\lambda\in N-1}p_{m}(E_{\lambda}^{N-1})+\sum_{\lambda\in N+1}p_{m}(E_{\lambda}^{N+1})=1$ due to the completeness of $\{|\Psi_{\lambda}^{N+1}\rangle\}_{\lambda}$ for the $(N+1)$-electron states and that of $\{|\Psi_{\lambda}^{N-1}\rangle\}_{\lambda}$ for the $(N-1)$-electron states, as expected.

If we input $|\Psi_{\mathrm{gs}}^{N}\rangle$ to the off-diagonal circuit $\mathcal{C}_{mm^{\prime}}$ in Fig. 2 and process the whole system the way described above, the ancillary bits $|0\rangle\otimes|1\rangle$ or $|1\rangle\otimes|1\rangle$ will be observed and the energy eigenvalue $E_{\lambda}^{N+1}$ will be obtained with probabilities [see eqs. (22) and (24)]

[TABLE]

while the ancillary bits $|0\rangle\otimes|0\rangle$ or $|1\rangle\otimes|0\rangle$ will be observed and the energy eigenvalue $E_{\lambda}^{N-1}$ will be obtained with probabilities [see eqs. (21) and (23)]

[TABLE]

This means that we can get the off-diagonal components of transition matrices $B_{\lambda}^{(\mathrm{e})}$ and $B_{\lambda}^{(\mathrm{h})}$ from eq. (18) via statistical sampling for a fixed combination of $m$ and $m^{\prime}$. It is easily confirmed that $\sum_{\sigma=+,-}[\sum_{\lambda\in N-1}p_{mm^{\prime}}(\sigma,E_{\lambda}^{N-1})+\sum_{\lambda\in N+1}p_{mm^{\prime}}(\sigma,E_{\lambda}^{N+1})]=1,$ as expected.

We provide the pseudocodes in Appendix A for the calculation process of GF explained above.

II.3 Galitskii–Migdal formula

The Galitskii–Migdal (GM) formulaFetter and Walecka (2003) enables one to calculate the ground-state energy of an interacting electronic system solely from the time-ordered GF. It can be rewritten to a tractable form for representation using restricted Hartree–Fock (RHF) orbitals asDahlen et al. (2006); Caruso et al. (2013); Phillips et al. (2015)

[TABLE]

where the integrand on the RHS contains a convergence factor $e^{+i\omega 0}$, forcing us to pick up the poles of the GF for the states below the Fermi level. $E_{\mathrm{nucl}}$ is the nuclear-repulsion energy. For spatial HF orbitals $p$ and $p^{\prime}$, $h_{pp^{\prime}}$ is the matrix element of the one-electron operator $h(\bm{r})$, which is the sum of the kinetic-energy term and the ionic-potential term. $\varepsilon$ is the diagonal matrix whose components are the HF orbital energies.

[TABLE]

is the one-particle density matrixFetter and Walecka (2003); Stefanucci and van Leeuwen (2013) for spin $\sigma$. $\Sigma_{\mathrm{c}}$ is the self-energy obtained from the Dyson equation $\Sigma_{\mathrm{c}\sigma}[G]=G_{\mathrm{HF}\sigma}^{-1}-G_{\sigma}^{-1},$ where the HF GF is given solely by the orbital energies: $G_{\mathrm{HF}\sigma pp^{\prime}}(\omega)=\delta_{pp^{\prime}}(\omega-\varepsilon_{p})^{-1}.$ $\Sigma_{\mathrm{c}}$ is responsible for the correlation effects in the (interacting) GF, which are not taken into account in the HF solution. We can use the expression for $E_{\mathrm{GM}}$ as an energy functional for an arbitrary input GF. If we substitute the HF GF into eq. (29), the third term on the RHS vanishes and we get the well known expression for the HF total energy, $E_{\mathrm{HF}}=E_{\mathrm{GM}}[G_{\mathrm{HF}}]=E_{\mathrm{nucl}}+\sum_{\sigma}\mathrm{Tr}[(h+\varepsilon)\gamma_{\mathrm{HF}\sigma}]/2,$ where $\gamma_{\mathrm{HF}\sigma}$ is the HF density matrix. The total energy for the interacting case is thus written as $E_{\mathrm{GM}}[G]=E_{\mathrm{HF}}+\Delta E_{1}[G]+\Delta E_{2}[G],$ where the sum of

[TABLE]

and

[TABLE]

is the correlation energy. $\Delta E_{1}$ is interpreted as the energy correction coming from the variation in the occupancy of HF orbitals, while an interpretation for $\Delta E_{2}$ within the HF picture is difficult to draw. We should keep in mind that $E_{\mathrm{GM}}[G_{\mathrm{trial}}]$ calculated from the GF $G_{\mathrm{trial}}$ for a trial ground state $|\Psi_{\mathrm{trial}}\rangle$ via eqs. (3) and (4) can differ from the expected energy in general: $E_{\mathrm{GM}}[G_{\mathrm{trial}}]\neq\langle\Psi_{\mathrm{trial}}|\mathcal{H}|\Psi_{\mathrm{trial}}\rangle$ , since eqs. (3) and (4) use the fact that the true ground state is an eigenstate of $\mathcal{H}$. We use the expressions in eqs. (31) and (32), however, to examine quantitatively the accuracy of GFs calculated in the present study. It is because one of our purposes is to see how $E_{\mathrm{GM}}$ values for UCC GFs from statistical sampling approach the ideal values as the number of measurements increases.

III Computational details

We adopted STO-3G basis sets as the Cartesian Gaussian-type basis functionsHelgaker et al. (2000) for all the elements in our quantum chemistry calculations. The Coulomb integrals between the atomic orbitals were calculated efficiently.Fermann and Valeev (2003) We first performed RHF calculations to get the orthonormalized molecular orbitals in the target systems and calculated the two-electron integrals between them, from which we constructed the second-quantized electronic Hamiltonians. After that, we used JW transformation to get the Hamiltonians in qubit representation by using OpenFermionMcClean et al. (2017) to perform full configuration-interaction (FCI) and UCC calculations. The parameters in the UCC calculations were optimized by employing the constrained optimization by linear approximation (COBYLA) method.

Although our scheme for the calculation of GF assumes that the energy spectra of $(N\pm 1)$-electron states for a target system are already known (see Procedure 1), which can be obtained in various approaches for quantum computers,Higgott et al. (2019); Jones et al. (2019); McClean et al. (2017); Santagati et al. (2018); Colless et al. (2018) we simply use those obtained in (classical) FCI calculations for the $(N\pm 1)$-electron states in the present study. It is because the main purpose is to demonstrate succinctly the validity of our scheme for GFs using statistical sampling. Simulations of GFs by taking into account the restrictions on the accuracy of spectra of excited states imposed by hardware and/or specific algorithms should be performed in the future. Our calculations of GFs, including those simulated with statistical sampling, were performed by substituting the necessary quantities into the Lehmann representation, given by eqs. (3) and (4). We set $\delta$ in eq. (7) to $0.02$ a.u. for the spectral functions throughout the present study.

For numerical evaluation of the integrals in eqs. (30) and (32), we adopted rectangular contours on the complex plane so that they encircle all the poles on the negative real axis for the integrands.

IV Results and discussion

IV.1 LiH molecule

IV.1.1 UCC calculations

By fixing the bond length at $1.6$ Å in an LiH molecule, we performed an RHF calculation and obtained $E_{\mathrm{RHF}}=-213.9322$ eV and six spatial orbitals among which the two lowest ones were fully occupied. Therefore we adopted the RHF solution as the reference state $|\Psi_{\mathrm{ref}}\rangle=a^{\dagger}_{1\downarrow}a^{\dagger}_{1\uparrow}a^{\dagger}_{0\downarrow}a^{\dagger}_{0\uparrow}|\mathrm{vac}\rangle\propto X_{3}X_{2}X_{1}X_{0}|0\cdots 0\rangle$ in the JW representation, where $X_{j},Y_{j},$ and $Z_{j}$ for $j=0,\dots,11$ are the Pauli matrices acting on the $j$th qubit, for the subsequent simulations of quantum computation with twelve qubits for the STO-3G basis (twelve) functions. We tried two excitation operators $T_{1}(\theta_{1},\theta_{2})=\theta_{1}a^{\dagger}_{2\downarrow}a^{\dagger}_{2\uparrow}a_{1\downarrow}a_{1\uparrow}+\theta_{2}a^{\dagger}_{5\downarrow}a^{\dagger}_{5\uparrow}a_{1\downarrow}a_{1\uparrow}$ and $T_{2}(\theta_{1},\theta_{2})=\theta_{1}a^{\dagger}_{3\downarrow}a^{\dagger}_{3\uparrow}a_{1\downarrow}a_{1\uparrow}+\theta_{2}a^{\dagger}_{4\downarrow}a^{\dagger}_{4\uparrow}a_{1\downarrow}a_{1\uparrow},$ each of which excites the two electrons in the highest occupied molecular orbital (HOMO), composed mainly of the Li $2s$ orbital, to the unoccupied orbital. [See Fig. 3(a)] We rewrite each of $T_{1}$ and $T_{2}$ to Pauli tensors for the qubits and pick up only a single tensor from them for each parameter as an approximation similarly to Hempel et al.,Hempel et al. (2018) which is then substituted into eq. (1) to define the ansatz. The ansätze in this case thus read

[TABLE]

for $T_{1}$ and

[TABLE]

for $T_{2}$, where we have rescaled the real parameters. We constructed the circuits $\mathcal{C}^{\mathrm{LiH}}_{1}$ and $\mathcal{C}^{\mathrm{LiH}}_{2}$ that act as these unitary operators and optimized the parameters to obtain the UCC ground-state energies. $\mathcal{C}^{\mathrm{LiH}}_{1}$ actually operates only on the eight among the twelve qubits, as shown in Fig. 4. It is similarly the case with $\mathcal{C}^{\mathrm{LiH}}_{2}$. The optimized $U_{1}$ gave $E_{\mathrm{UCC1}}=-214.3323$ eV, closer to the FCI value $E_{\mathrm{FCI}}=-214.4889$ eV than the optimized $U_{2}$ did with $E_{\mathrm{UCC2}}=-213.9758$ eV.

IV.1.2 GFs exact within UCC

We calculated the GFs from the ground states of the FCI and optimized UCC solutions, as shown in Fig. 5 (a). The FCI spectra $A^{\mathrm{FCI}}(\omega)$ exhibit the weak satellite peaks, which are correlation effects and thus are absent in the HF spectra. Specifically, the weak peaks are seen for $-30<\omega<-15$ eV, $5<\omega<12.5$ eV, and $17.5$ eV $<\omega$. Although the major peaks, called the quasiparticle peaks, can be basically assigned to the individual HF orbitals, the two neighboring major peaks around $\omega=15$ are split due to the correlation effects on the HF orbital $5$. The satellite peaks are also seen in the UCC spectra $A^{\mathrm{UCC}}(\omega)$ for both $U_{1}$ and $U_{2}$. The quasiparticle peaks in the FCI spectra are closer to the Fermi level ($\omega=0$) than the HF orbital energies are, which is due to the well known fact that HF solutions overestimate energy gaps in general. The overall shapes of the FCI and UCC spectra look quite similar to each other despite the simple ansätze since an LiH molecule is a weakly correlated system. The locations of quasiparticle and satellite peaks in the UCC spectra for the optimized $U_{1}$ are closer to those in the FCI spectra than those for the optimized $U_{2}$, as expected.

IV.1.3 UCC GFs via statistical sampling

Hereafter we denote the ground state for the optimized $U_{1}$ simply by the UCC ground state $|\Psi_{\mathrm{gs}}^{N(\mathrm{UCC})}\rangle$. To simulate the scheme for obtaining GFs on a quantum computer proposed above, we calculated the transition matrix elements between $|\Psi_{\mathrm{gs}}^{N(\mathrm{UCC})}\rangle$ and the FCI energy eigenstates $|\Psi_{\lambda}^{N\pm 1(\mathrm{FCI})}\rangle$. We generated random numbers according to these values since they represent the probability distributions of the measurement results for the qubits. [See eqs. (25)-(28)] By building the histograms of the results of simulated measurements, we constructed the GF $G^{\mathrm{UCC-stat}}$ for the UCC ground state. We denote such construction of all the components of a GF by a single simulation of GF in what follows.

Typical spectral functions $A^{\mathrm{UCC-stat}}(\omega)$ simulated in this way are shown in Fig. 5(b). We can see that the quasiparticle peaks in $A^{\mathrm{UCC}}(\omega)$ are well reproduced by the statistical sampling even for the smaller $N_{\mathrm{meas}}$. For the satellite peaks, on the other hand, their shapes for the two values of $N_{\mathrm{meas}}$ can be quite different from each other. In particular, those near $\omega=-25$ and $20$ eV were not even detected for $N_{\mathrm{meas}}=10$ due to the too few measurements. These observations indicate that a number of measurements on a quantum computer have to be performed if one wants to capture the correlation effects accurately, just as PES experiments and their inverse have to be conducted many times for the rare physical processes.

Figure 6(a) shows the typical shapes of the traces of self-energies $\Sigma_{\mathrm{c}}$ calculated from $G^{\mathrm{UCC-stat}}$ with $N_{\mathrm{meas}}=1000$ and $8000$. We notice that the convergence of self-energy with respect to $N_{\mathrm{meas}}$ looks far from satisfaction even for $N_{\mathrm{meas}}=1000$, in contrast to the sampled GF. [See Fig. 5 (b)] This observation comes from the fact that the major contributions to the GF, nothing but the quasiparticle peaks, are already taken into account as the HF GF, while the presence of $\Sigma_{\mathrm{c}}$ results solely from the correlation effects.

IV.1.4 Correlation energy from GF

To examine the statistical behavior of $G^{\mathrm{UCC-stat}}$ quantitatively, we performed 100 simulations to obtain $G^{\mathrm{UCC-stat}}$ for each given value of $N_{\mathrm{meas}}$ and calculated the correlation energies by using the GM formula in eqs. (31) and (32). The results for $N_{\mathrm{meas}}=1000,2000,4000,8000,16000$, and $32000$ are shown in Fig. 6(b), where $\Delta E_{1}[G^{\mathrm{UCC-stat}}]$ and $\Delta E_{2}[G^{\mathrm{UCC-stat}}]$ scatter around the ideal values, $\Delta E_{1}[G^{\mathrm{UCC}}]$ and $\Delta E_{2}[G^{\mathrm{UCC}}]$, respectively. The deviations of the sampled values from the ideal values decrease as $N_{\mathrm{meas}}$ increases, as expected.

IV.2 H2O molecule

IV.2.1 UCC calculations

By fixing the O-H bond length at $0.96$ Å and the H-O-H bond angle at $104.5^{\circ}$ in an H2O molecule, we performed an RHF calculation and obtained $E_{\mathrm{RHF}}=-2039.8504$ eV and seven spatial orbitals among which the five lowest ones were fully occupied. Therefore we adopted the RHF solution as the reference state $|\Psi_{\mathrm{ref}}\rangle\propto X_{9}\cdots X_{0}|0\cdots 0\rangle$ in the JW representation, for the subsequent simulations of quantum computation with fourteen qubits for the STO-3G basis (fourteen) functions. We tried two excitation operators

[TABLE]

and

[TABLE]

each of which excites the electrons in the MOs near the Fermi level, composed mainly of the O $2p$ orbitals, to the unoccupied orbitals. [See Fig. 3(b)] We rewrite them to the qubit operators with approximations similarly to the case of an LiH molecule and introduced the ansätze

[TABLE]

for $T_{2}$ and

[TABLE]

for $T_{1}$, where we have rescaled the real parameters. We constructed the circuits $\mathcal{C}^{\mathrm{H}_{2}\mathrm{O}}_{1}$ and $\mathcal{C}^{\mathrm{H}_{2}\mathrm{O}}_{2}$ that act as these unitary operators and optimized the parameters to obtain the UCC ground-state energies. The optimized $U_{1}$ gave $E_{\mathrm{UCC1}}=-2040.4359$ eV, closer to the FCI value $E_{\mathrm{FCI}}=-2041.2013$ eV than the optimized $U_{2}$ did with $E_{\mathrm{UCC2}}=-2040.0492$ eV.

IV.2.2 GFs exact within UCC

We calculated the GFs from the ground states of the FCI and optimized UCC solutions, as shown in Fig. 7 (a). These three spectral functions admit analyses similar to those for an LiH molecule described above, since an H2O molecule is also a weakly correlated system.

IV.2.3 UCC GFs via statistical sampling

Hereafter we denote the ground state for the optimized $U_{1}$ simply by the UCC ground state. We performed simulations for obtaining GFs via statistical sampling in the same manner as in the case of an LiH molecule. Typical simulated spectral functions are shown in Fig. 7(b). We can see that the quasiparticle peaks in $A^{\mathrm{UCC}}(\omega)$ are well reproduced by the statistical sampling, while the sampled satellite peaks are not satisfactory. These results are similar to those in the LiH case.

Figure 8(a) shows the typical shapes of the traces of self-energies $\Sigma_{\mathrm{c}}$ calculated from $G^{\mathrm{UCC-stat}}$ with $N_{\mathrm{meas}}=32000$ and $64000$. The convergence of self-energy with respect to $N_{\mathrm{meas}}$ is found to be much slower than in the LiH case. This slow convergence propagates to that of the sampled correlation energies, as explained below.

IV.2.4 Correlation energy from GF

Similarly to the case of an LiH molecule, we performed 100 simulations to obtain $G^{\mathrm{UCC-stat}}$ for each given value of $N_{\mathrm{meas}}$ and calculated the correlation energies by using the GM formula, as shown in Fig. 8(b). Although the increase in $N_{\mathrm{meas}}$ leads to the convergence of sampled correlation energy as well as for an LiH molecule, the convergence for this case is much slower. $N_{\mathrm{meas}}=32000$ achieves the convergence of $\Delta E_{1}+\Delta E_{2}$ within about $0.2$ eV accuracy for an LiH molecule [see Fig. 6(b)], while the same $N_{\mathrm{meas}}$ only achieves an accuracy as large as $1.5$ eV for an H2O molecule. These observations reflect the generic fact that the increase in the number of electrons immediately means the rapid increase in the excitation channels, which forces us to perform measurements on a quantum computer much more times to reproduce the correct probability distribution. Although the accuracy achieved in our simulations is far from the chemical accuracy, 1 kcal/mol $\approx$ 43 meV, it seems that we are left with much room for improving the naïve scheme proposed in the present study. In particular, the pursuit of efficient construction of histograms leading to the suppression of the rapid increase in the necessary number of measurements is valuable in the future.

V Conclusions

We proposed a scheme for the construction of one-particle GF of an interacting electronic system via statistical sampling on a quantum computer. We were able to circumvent the restriction of unitarity of qubit operations by introducing the quantum circuits for probabilistic state preparation. The quantum circuits for the diagonal and off-diagonal components and the subsequent QPE allow us to calculate the GF straightforwardly from the histogram obtained via measurements on the qubits.

For demonstrating the validity of our scheme, we performed simulations of such construction of GFs for LiH and H2O molecules based on the UCC method by referring to the spectral functions exact within UCC. We found that the accurate reproduction of weaker satellite peaks requires more measurements to detect the small contributions to the spectra. We also examined the accuracy of sampling method by exploiting the GM formula to find that the increase in the number of electrons leads to the rapid increase in the excitation channels, which forces us to perform measurements many times to get a correct histogram.

We should keep in mind that our simulations were performed on the assumption that the many-electron energy eigenvalues of the target systems are known and the QPE experiments are conducted with no probabilistic error. The results in the present study thus indicate that the requirements of resources for the accurate description of correlation effects using a real quantum computer grow rapidly as the target systems become large, as long as we use the simple statistical sampling. Therefore we have to improve the scheme for obtaining GFs accurately by considering more realistic setups and simultaneously reducing costs in the future.

Acknowledgements.

This research was supported by MEXT as Exploratory Challenge on Post-K computer (Frontiers of Basic Science: Challenging the Limits) and Grants-in-Aid for Scientific Research (A) (Grant Numbers 18H03770) from JSPS (Japan Society for the Promotion of Science).

Appendix A pseudocodes for GF

Here we provide the pseudocodes for the calculation process of GF proposed in the present study. We assume that not only the energy of $N$-electron ground state but also the energy eigenvalues of $(N\pm 1)$-electron states have been obtained before entering the calculation process for GF. The main process, CalcGF, is given by procedure 1. CalcAmplsDiag in procedure 2 is called to calculate the diagonal components of transition matrices, while CalcAmplsOffDiag in procedure 3 is called to calculate the off-diagonal components. The latter calls CalcAmplsAux in procedure 4 to get $D^{(\mathrm{e,h})\pm}_{mm^{\prime}}$, from which the off-diagonal components $B^{(\mathrm{e,h})}_{mm^{\prime}}$ are calculated using eq. (18).