Linear-response functions of molecules on a quantum computer: Charge and spin responses and optical absorption

We propose a scheme for the construction of linear-response functions of an interacting electronic system via quantum phase estimation and statistical sampling on a quantum computer. By using the unitary decomposition of electronic operators for avoiding the difﬁculty due to their nonunitarity, we provide the circuits equipped with ancillae for probabilistic preparation of qubit states on which the necessary nonunitary operators have acted. We perform simulations of such construction of the charge and spin response functions and photoabsorption cross sections for C 2 and N 2 molecules by comparing with the results from full conﬁguration-interaction calculations. It is found that the accurate detection of subtle structures coming from the weak poles in the response functions requires a large number of measurements.


I. INTRODUCTION
Since the information carrier of a programmable quantum computer is a set of qubits that exploits the principle of superposition, essentially parallel algorithms can exist and perform computation for classically formidable problems [1,2].Quantum chemistry [3] is believed to be one of the most suitable research fields for quantum computation since its problem setting is quantum mechanical by definition.Indeed, a quantum computer can treat a many-electron state composed of lots of Slater determinants as it is in a sense that the electronic state is encoded as a superposition of qubit states via an appropriately chosen map such as the Jordan-Wigner (JW) [4] and Bravyi-Kitaev (BK) [5] transformations.
The quantity which a quantum chemistry calculation is asked to first provide is the total energy of a target system [6].One of the most widespread methods for obtaining the total energy is the variational quantum eigensolver (VQE), in which a trial many-electron state is prepared via a parametrized quantum circuit.The parameters are optimized iteratively with the aid of a classical computer aiming at the ground state.This approach was first realized [7] by using a quantum photonic device, after which the realizations by superconducting [8,9] and ion trap [10] quantum computers have been reported.There exist algorithms for obtaining the energy spectra of excited states [11][12][13][14][15][16].
Not only academic interest but also industrial demands for accurate explanations and predictions of material properties Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.make it an urgent task to develop methodologies on quantum computers for various electronic properties other than energy levels.The one-particle Green's functions (GFs) are important in correlated electronic systems [17,18] since they are often used in explanation of spectra measured in photoemission and its inverse experiments.Recently we proposed [19] a method for the GFs [20,21] via statistical sampling by employing quantum phase estimation (QPE) [1].Endo et al. [22], on the other hand, proposed a method for the GFs by focusing on noisy intermediate-scale quantum (NISQ) devices.
The charge and spin response functions, formulated in the linear-response theory [23,24], describe the leading contributions to the electric and magnetic excitations when perturbation fields are applied to a target system.Since the response functions are the fundamental building blocks in constructing the elaborated methods for correlated electrons such as GW theory [23,25], the accurate calculation of them is needed.In addition, there exist response functions directly related to measurable quantities such as dielectric constants, electric conductivities, and magnetic susceptibilities.
Given the recent rapid development of fabrication techniques for quantum hardware and the growing demands for quantum computation in material science, it is worth making tools for analyses on correlation effects.In this study, we propose a scheme for the construction of the response functions of an interacting electronic system via statistical sampling on a quantum computer.We track the changes in the state of qubits analytically when they undergo gate operations and measurements to derive the exact probability distributions of the outcomes.For examining the validity of our scheme, we perform the full configuration-interaction (FCI) calculations for diatomic molecules, from which we simulate the measurements on qubits to get the response functions by generating random numbers based on the exact probability distributions.
Circuit C O for the probabilistic preparation of a state on which a linear combination O of 2 n unitaries has acted.H in the circuit represents the Hadamard gate.The desired state is prepared according to a measurement outcome of n ancillary qubits.This circuit contains C (n) as a partial circuit defined in Fig. 2.This paper is organized as follows.In Sec.II, we explain the basic ideas and our scheme in detail by providing the quantum circuits for obtaining the response functions via statistical sampling.In Sec.III, we describe the computational details for our simulations on a classical computer.In Sec.IV, we show the simulation results for C 2 and N 2 molecules.In Sec.V, we provide the conclusions.

A. Circuit for a linear combination of unitaries
The calculation of a physical quantities of interest quite often involves the evaluation of matrix elements of various electronic operators.Such an operator is, however, not necessarily unitary and it prohibits one from implementing it straightforwardly as logic gates for qubits.We describe a workaround for this difficulty by providing a circuit for probabilistic state preparation.
For 2 n unitaries U k , where k is a bit string of length n, we want to apply an arbitrary linear combination of them, to an input register |ψ .c k ≡ |c k | exp(iφ k ) is the complex coefficient.We construct a circuit C O equipped with n ancillary qubits, as depicted in Fig. 1, containing a partial circuit C (n)  defined recursively in Fig. 2. If the measurement outcome for Recursive definitions of circuits C ( j) λ for (a) j = 0 and (b) j > 0. λ is empty or a bit string.exp(iφ λ ) is an identity operator multiplied by the phase factor.R y (θ ) = e −iθσy/2 is a rotation for an angle θ .These circuits are used in Fig. 1.
the ancillae is |0 ⊗n , the state of the whole system collapses to |0 ⊗n ⊗ O|ψ up to a normalization constant.The proof of this fact is provided in Appendix A.
The state preparation techniques for the response functions described below are special cases of the one introduced above with some target-specific modifications.

Linear-response functions
We work with the n orbs orthonormalized spatial oneelectron orbitals for each spin direction in a target N-electron system, which we assume is in the many-electron ground state | gs at zero temperature.The formalism described below can be easily extended to a system with multiple ground states and/or a nonzero temperature.
The response functions in terms of Hermitian operators O and O in time domain is defined as [23] where O(t ) is the operator in the Heisenberg picture and the expectation value is for the ground state.We assume that the system has been in the equilibrium when the perturbation is turned on.Since the response function defined in Eq. ( 2) depends on time only via t − t in this case, its expression in frequency domain is written as for a real frequency ω.
is the Lehmann summation over the energy eigenstates for a complex frequency z.E N gs is the ground-state energy and E λ is the energy eigenvalue of the λth excited state | λ .
is the transition matrix element, which satisfies clearly L λOO = (L λO O ) * .The positive infinitesimal constant δ appears in Eq. ( 3) due to the retarded nature of the response function, rendering all the poles immediately below the real axis.It is clear that the real part of χ OO (ω) is even with respect to ω, while the imaginary part is odd.

Charge and spin responses
By using the creation a † pσ and annihilation a pσ operators for the pth spatial orbital of σ spin |φ pσ (σ = α, β ), the electron number operator is given by n pσ = a † pσ a pσ .The spin operator is given by s p = σ,σ a † pσ (σ el σ σ /2)a pσ , where σ el is the Pauli matrix for the electronic state.
The orbital-wise response functions involving the charge and spin of the individual spin orbitals are obtained by putting O pn ≡ σ n pσ and O p j ≡ s p j ( j = x, y, z), respectively, into Eq.(4).To rewrite the expression for transition matrix elements in Eq. ( 5) into a more tractable form, we define the charge-charge transition matrix element, for the λth energy eigenstate | λ of the N-electron states.We define similarly the spin-spin one, S λp j,p j ≡ gs |s p j | λ λ |s p j | gs , (7) for j, j = x, y, z and the spin-charge one, M λp j,p σ ≡ gs |s p j | λ λ |n p σ | gs ≡ M * λp σ ,p j .(8) From these matrix elements, we define the following Hermitian matrices L λ : L λp j,p j ≡ S λp j,p j , (10) which are used for Eq. ( 5).

Generic one-body operators
A generic one-body operator is given in the form, where m and m are the composite indices of spatial orbitals and spins.O mm is the matrix element between the oneelectron orbitals.Our scheme is actually applicable to such a generic case for O and O , for which we have to obtain the transition matrix elements, which satisfy clearly B λm 1 m 2 ,m 3 m 4 = (B λm 4 m 3 ,m 2 m 1 ) * .From them, we sum up the contributions to Eq. ( 5) as We calculate the electric polarizabilities of molecules in the present study as examples for the generic scheme.The contribution from electrons to the electric dipole of a molecule is d = − m,m d mm a † m a m , where the negative sign on the righthand side comes from the electron charge.d mm ≡ φ m |r|φ m is the matrix element of position operator.The linear electric polarizability tensor α j j ( j, j = x, y, z) [26,27] is defined as the first derivative of the expected electric dipole with respect to an external electric field, obtained via the response function as α j j (ω) = −χ d j d j (ω), from which the photoabsorption cross section [28] is given by This cross section is a measurable quantity in optical absorption experiments.

Unitary decomposition of electronic operators
Although there are alternatives for mapping the electronic operators of a target system to the qubit ones such as JW [4] and BK [5] transformations, we do not distinguish between an electronic operator and its corresponding qubit representation in what follows since no confusion will occur for the readers.For each combination of a spatial orbital p and a spin σ , we perform the Majorana fermion-like [29] transformation for the qubits [19]: and which are unitary regardless of the adopted qubit representation thanks to the anticommutation relation between the electronic operators and can thus be implemented as logic gates in the quantum computer.This means that we can prepare at least probabilistically an electronic state on which an arbitrary product of the creation and annihilation operators has acted, similarly to the case for GFs [19].
In what follows, we assume that the many-electron ground state | gs is already known and can be prepared on a quantum computer.

C. Charge-charge responses
Let us first consider the determination of the charge-charge transition matrices N λ .

Circuits for diagonal components
From the unitary operators in Eqs. ( 16) and ( 17) for a combination of a spatial orbital p and a spin σ , we define U (p)  κκ σ ≡ U κ pσ U κ pσ for κ, κ = 0, 1, which are also unitary.With them, the electron number operator is written as while the whole number operator is written as We construct a circuit C pσ equipped with two ancillary qubits |q A 0 and |q A 1 by implementing the controlled operations of U (p)  κκ σ , as depicted in Fig. 3.The whole system consists of the ancillae and an arbitrary input register |ψ .Its state changes by undergoing the circuit as pσ,p σ |ψ and n ± pσ,p σ |ψ from an arbitrary input state |ψ and three ancillary qubits.Z (±π/4) = diag(1, e ±iπ/4 ) is a phase gate.The partial circuits defined in Fig. 3 are contained as the controlled subroutines.
since (a † pσ ) 2 and a 2 pσ vanish due to the Fermi statistics.The projective measurement [1] on |q A 0 is represented by the two operators P q = I ⊗ |q q| ⊗ I (q = 0, 1).The state of the whole system collapses immediately after the measurement as follows:

Circuits for off-diagonal components
For mutually different combinations (p, σ ) and (p , σ ) of spatial orbitals and spins, we define the four non-Hermitian auxiliary operators, and Un-normalized auxiliary states with the energy eigenstates, from which the charge-charge transition matrix elements in Eq. ( 6) can be calculated as We construct a circuit C pσ,p σ equipped with three ancillary qubits by using the controlled operations of the partial circuits U (p)  σ and U (p ) σ , defined in Fig. 3, as depicted in Fig. 4. The whole system consists of the ancillae and an arbitrary input register |ψ .Its state changes by undergoing the circuit as The projective measurement on |q A 2 and |q A 1 is represented by the four operators P qq = |q q| ⊗ |q q | ⊗ I ⊗ I (q, q = 0, 1).The two outcomes among the possible four are of our interest, immediately after which the whole system collapses as follows:

Transition matrices via statistical sampling
Given the result of a measurement on the ancillary bits for a diagonal or an off-diagonal component, we have the register | ψ different from the input N-electron state.Then we perform QPE for the Hamiltonian H by inputting | ψ to obtain one of the energy eigenvalues in the Hilbert subspace for the N-electron states.A QPE experiment inevitably suffers from probabilistic errors that depend on the number of qubits and the various parameters for the Suzuki-Trotter decomposition of H.We assume for simplicity, however, that the QPE procedure is realized on a quantum computer with ideal precision as well as in our previous study [19].We will thus find the estimated value to be E N λ with a probability | λ | ψ | 2 [1].Since the probabilistic state preparation and obtaining the transition matrix elements are dominated mathematically by multinomial distributions, the errors in the response function scale as N −1/2 meas for N meas repeated measurements.If we input | gs to the diagonal circuit C pσ in Fig. 3 and observe the ancillary bit |q A 0 = 0 for QPE, the energy eigenvalue E N λ will be obtained with a probability [see Eq. ( 21)], This means that we can get the diagonal components of transition matrices N λ via statistical sampling for a fixed combination of p and σ .
If we input | gs to the off-diagonal circuit C pσ,p m in Fig. 4

and observe the ancillary bits |q
1 for QPE, the energy eigenvalue E N λ will be obtained with probabilities [see Eqs. ( 27) and ( 28)], This means that we can get the off-diagonal components of N λ from Eq. ( 25) via statistical sampling for a fixed combination of p, p , σ , and σ .
|ψ FIG. 5. Spin-spin diagonal circuit C p j ( j = x, y) for probabilistic preparation of s p j |ψ and s p j |ψ from an arbitrary input state |ψ and an ancillary qubit.We define the partial circuit U (p)  j by enclosing it with dashed lines.

D. Spin-spin responses
Let us next consider the determination of the spin-spin transition matrix elements S λp j,p j for j, j = x, y.Since the operators of the z component of spin and the electron number are related via s pz = (n pα − n pβ )/2, the matrix elements S λpz,p z can be calculated from the diagonal and off-diagonal components of N λ , already obtained in Sec.II C. S λp j,p z can be calculated from M λ which will be explained in Sec.II E.

Circuits for diagonal components
From the unitary operators in Eqs. ( 16) and ( 17) for a combination of a spatial orbital p and a spin σ , we define the following four unitary operators: With them, the spin operators for the x and y directions are written as for j = x, y.We define for later convenience.We construct a circuit C p j equipped with an ancillary qubit by implementing the controlled operations of U (p)  0 j and U (p) 1 j , as depicted in Fig. 5.The whole system consists of the ancilla and an arbitrary input register |ψ .Its state changes by undergoing the circuit as The projective measurement on |q A is represented by the two operators P q = |q q| ⊗ I (q = 0, 1).The state of the whole system collapses immediately after the measurement as follows: FIG. 6. Spin-spin off-diagonal circuit C p j,p j ( j, j = x, y) for probabilistic preparation of s ± p j,p j |ψ and s ± p j,p j |ψ from an arbitrary input state |ψ and two ancillary qubits.The partial circuits defined in Fig. 5 are contained as the controlled subroutines.

Circuits for off-diagonal components
For mutually different combinations (p, j) and (p , j ) of spatial orbitals and spin components j, j = x, y, we define the four non-Hermitian auxiliary operators, and Un-normalized auxiliary states | ± p j,p j ≡ s ± p j,p j | gs can have overlaps T ± λp j,p j ≡ | λ | ± p j,p j | 2 with the energy eigenstates, from which the spin-spin transition matrix elements in Eq. ( 7) can be calculated as S λp j,p j = e −iπ/4 (T + λp j,p j − T − λp j,p j ) + e iπ/4 (T + λp j ,p j − T − λp j ,p j ).
We construct a circuit C p j,p j equipped with two ancillary qubits by using the controlled operations of the partial circuits U (p) j and U (p ) j , defined in Fig. 5, as depicted in Fig. 6.The whole system consists of the ancillae and an arbitrary input register |ψ .Its state changes by undergoing the circuit as The projective measurement on |q A 1 and |q A 0 is represented by the four operators P qq = |q q| ⊗ |q q | ⊗ I (q, q = 0, 1).The two outcomes among the possible four are of our interest, immediately after which the whole system collapses as follows: x, y and σ = α, β ) for probabilistic preparation of v ± p j,p σ |ψ and other states from an arbitrary input state |ψ and three ancillary qubits.The partial circuits defined in Figs. 3 and 5 are contained as the controlled subroutines.

Transition matrices via statistical sampling
We can get the transition matrices S λ via statistical sampling similarly to the charge-charge ones.If we input | gs to the diagonal circuit C p j in Fig. 5 followed by a measurement and QPE for H, the energy eigenvalue E λ will be obtained with a probability 4S λp j,p j .[See Eqs. ( 7) and (34).]If we use the off-diagonal circuit C p j,p j in Fig. 6, on the other hand, the energy eigenvalue E λ will be obtained with probabilities 4T ± λp j,p j depending on the measurement outcome.[See Eqs. ( 40) and ( 41).]The off-diagonal components of transition matrices are then calculated from Eq. ( 38).

E. Spin-charge responses
Having found ways to determine the charge-charge and spin-spin contributions, let us consider the determination of the spin-charge transition matrices M λp j,p σ for j = x, y and σ = α, β.Those involving the z component of spin, M λpz,p σ , can be calculated from the diagonal and off-diagonal components of N λ , already obtained in Sec.II C.

Circuits for off-diagonal components
For combinations (p, j) and (p , σ ) of spatial orbitals and spin components j = x, y with σ = α, β, we define the two non-Hermitian auxiliary operators, Un-normalized auxiliary states with the energy eigenstates, from which the spin-charge transition matrix elements in Eq. ( 8) can be calculated as We construct a circuit C ± p j,p σ equipped with three ancillary qubits by using the controlled operations of the partial circuits U (p) j in Fig. 5 and U (p ) σ in Fig. 3, as depicted in Fig. 7.The whole system consists of the ancillae and an arbitrary input register |ψ .Its state changes by undergoing the circuit as The projective measurement on |q A 2 and |q A 0 is represented by the four operators P qq = |q q| ⊗ I ⊗ |q q | ⊗ I (q, q = 0, 1).Only one of the four possible outcomes is of our interest, immediately after which the whole system collapses as follows:

Transition matrices via statistical sampling
We can get the transition matrices M λ via statistical sampling similarly to the charge-charge and spin-spin ones.If we input | gs to the off-diagonal circuit C ± p j,p σ in Fig. 7 followed by a measurement and QPE for H, the energy eigenvalue E λ will be obtained with a probability T ± λp j,p σ .[See Eq. ( 45).]The off-diagonal components of transition matrices are then calculated from Eq. (43).
We provide the pseudocodes in Appendix B for the calculation procedures of response functions explained above.

F. Generic cases
Here we describe briefly the scheme for the response function involving generic one-body operators given as Eq. ( 12).

Circuits for diagonal components
For a combination of spin orbitals m 1 and m 2 , we define the four unitary operators U (m 1 m 2 ) κκ ≡ U κm 1 U κ m 2 (κ, κ = 0, 1).We construct a circuit C m 1 m 2 equipped with two ancillary qubits by implementing the controlled operations for an arbitrary input register |ψ , as depicted in Fig. 8.It is easily confirmed that the whole state changes by undergoing the circuit as 46) When the target ground state is input, the projective measurement on the two ancillae leads to the probabilistic preparation and other states from an arbitrary input state |ψ and three ancillary qubits.The partial circuits defined in Fig. 8 are contained as the controlled subroutines. of a † m 1 a m 2 | gs , for which the subsequent QPE process gives the "diagonal" matrix element B λm 2 m 1 ,m 1 m 2 [see Eq. ( 13)] via statistical sampling.
The charge-charge diagonal circuit C pσ in Fig. 3 is a special case of the generic circuit described here.

Circuits for off-diagonal components
For mutually different combinations (m 1 , m 2 ) and (m 3 , m 4 ) of spin orbitals, we define the two non-Hermitian auxiliary operators, with the energy eigenstates, from which the "off-diagonal" matrix element is calculated as We construct a circuit C m 2 m 1 ,m 3 m 4 equipped with three ancillary qubits by using the controlled operations of the partial circuits defined in Fig. 8, as depicted in Fig. 9, for an arbitrary input register |ψ .It is easily confirmed that the whole state changes by undergoing the circuit as When the target ground state is input, the projective measurement on the three ancillae leads to the probabilistic preparation of | ± m 2 m 1 ,m 3 m 4 .This circuit and the similarly constructed C m 3 m 4 ,m 2 m 1 allows one to calculate B λm 1 m 2 ,m 3 m 4 from Eq. (48) via statistical sampling.After obtaining the necessary matrix elements, they are put into Eq.( 14) to calculate the Lehmann summation.
The charge-charge off-diagonal circuit C pσ,p σ in Fig. 4 is a special case of the generic circuit described here.

III. COMPUTATIONAL DETAILS
As stated in Introduction, we simulated the measurements on qubits for obtaining the response functions of molecules by generating random numbers based on the analytically derived exact probability distributions.The computational details are described here.We adopted STO-6G basis sets as the Cartesian Gaussiantype basis functions [30] for all the elements in our quantum chemistry calculations.The two-electron integrals between the atomic orbitals (AOs) were calculated efficiently [31].We first performed restricted Hartree-Fock (RHF) calculations to get the orthonormalized molecular orbitals (MOs) in the target systems and calculated the two-electron integrals between them, from which we constructed the second-quantized electronic Hamiltonians.
In the FCI calculations for the large target Hilbert subspaces, we performed exact diagonalization of the electronic (not in qubit representation) Hamiltonians by using the Arnoldi method [32].We can take the z axis as the quantization axis for spins without loss of generality since our calculations are nonrelativistic.
We calculated the FCI response functions simply by substituting the necessary quantities into Eq.( 3).For the simulations of response functions from statistical sampling, we generated random numbers according to the matrix elements between the FCI energy eigenstates to mimic the measurements on ancillae and the ideal QPE procedures.We set δ in Eq. ( 3) to 0.01 atomic unit throughout the present study.

A. C 2 molecule
We used the experimental bond length of 1.242 Å [33] for a C 2 molecule in the RHF calculation and obtained the total energy E RHF = −2045.2939eV.This system contains six electrons per spin direction which occupy the lowest six MOs, as shown in Fig. 10(a).We found via the subsequent FCI calculation with E FCI = −2052.6918eV that the major electronic configuration in the nondegenerate many-electron   for N meas =5000 σ FCI-stat   for N meas =40000 σ(ω) (a.u.) 12. The photoabsorption cross section σ FCI exact within the FCI solution for a C 2 molecule.The simulated ones σ FCI−stat for the numbers of measurements N meas = 5000 and 40 000 are also shown.ground state, denoted by X 1 + g in spectroscopic notation [34], is the same as in the RHF solution.The ground state was found via the exact diagonalization of the Hilbert subspace for n α = n β = 6, from which we obtained the lowest 2000 among the 44 100 energy eigenvalues.
We calculated the response functions χ FCI exact within the FCI solution, from which the components χ FCI 1πn,1πn and χ FCI 1πz,3σ z are plotted in Fig. 11.We also performed simulations of statistical sampling for the construction of response function χ FCI−stat based on our scheme and plotted those for N meas = 10 000 and 40 000.It is seen that χ FCI 1πn,1πn in Fig. 11(a), which involves only the HOMO, is well reproduced by χ FCI−stat 1πn,1πn with N meas = 10 000.On the other hand, χ FCI 1πz,3σ z in Fig. 11(b) is not accurately reproduced by χ FCI−stat 1πn,1πn with N meas = 10 000.These results mean that the response involving a weak excitation channel requires a large number of measurements for its accurate reproduction, just like the situation for the GFs [19].We calculated the photoabsorption cross section σ FCI from the FCI solution, as shown in Fig. 12.By simulating the process for a generic one-body operator to obtain the matrix elements in Eq. ( 13), we also calculated the cross section σ FCI−stat for N meas = 5000 and 40 000 and plotted it in the figure.It is seen that those from sampling can take negative values despite the original definition in Eq. ( 15), which is ensured to be non-negative.It is due to our naïve implementation of Eq. (48) by using random numbers without considering any symmetry.The major peaks in σ FCI were well reproduced by σ FCI−stat with the smaller N meas , while the detailed structure between them required the larger N meas for the accurate reproduction.Since symmetry-adapted construction of transition matrix elements should reduce the necessary total number of measurements and physically appropriate results, such techniques will be useful as well as in VQE calculations [35].In a more practical setup where the QPE procedure is not ideal, the resolution and the reliability of two-qubit gates in the QPE procedure will have significant impacts on the results.The adequate number of measurements should thus be discussed by considering specific setups.

B. N 2 molecule
We used the experimental bond length of 1.098 Å [36] for an N 2 molecule in the RHF calculation and obtained the total energy E RHF = −2953.5952eV.This system contains seven electrons per spin direction which occupy the lowest seven MOs, as shown in Fig. 10(b).We found via the subsequent FCI calculation with E FCI = −2957.9124eV that the major electronic configuration in the nondegenerate many-electron ground state is X 1 + g [34], the same as in the RHF solution.The ground state was found via the exact diagonalization of the Hilbert subspace for n α = n β = 7, from which we obtained the lowest 2000 among the 14 400 energy eigenvalues.
We calculated the response functions from the FCI solution, from which the components χ FCI 3σ n,3σ n and χ FCI 3σ n,1π * n are plotted in Fig. 13.We also performed simulations for the construction of χ FCI−stat and plotted those for N meas = 10 000 and 40 000.It is seen that χ FCI 3σ n,3σ n in Fig. 13(a), which involves only the HOMO, is not well reproduced by χ FCI−stat 3σ n,3σ n with N meas = 10 000, in contrast to the C 2 molecule case.Since the strength of χ FCI 3σ n,1π * n in Fig. 13(b) is similar to that of χ FCI 3σ n,3σ n , N meas = 10 000 is not sufficient for the accurate reproduction of correct values as well.We found that χ FCI 3σ z,1π * z (not shown) is much weaker and even N meas = 40 000 is insufficient for obtaining good χ FCI−stat 3σ z,1π * z .We calculated the photoabsorption cross section σ FCI from the FCI solution, as shown in Fig. 14.Those from simulations of measurements σ FCI−stat for N meas = 5000 and 40 000 are also shown in the figure.

V. CONCLUSIONS
We proposed a scheme for the construction of linearresponse functions of an interacting electronic system via QPE and statistical sampling on a quantum computer.By using the unitary decomposition of electronic operators for avoiding the difficulty due to their nonunitarity, we provided the circuits equipped with at most three ancillae for probabilistic preparation of qubit states on which the necessary nonunitary operators have acted.We performed simulations of such construction of the response functions for C 2 and N 2 molecules by comparing with the accurate ones based on the FCI calculations.It was found that the accurate detection of subtle structures coming from the weak poles in the response functions requires a large number of measurements.
Since the unitary decomposition of electronic operators is applicable regardless of the adopted qubit representation, an electronic state on which an arbitrary product of the creation and annihilation operators has acted can be prepared at least probabilistically.The quality of results for our approach comes mainly along with that of QPE procedure, implying that the feasibility of our approach may be promising for small systems which have been already treated on realized quantum computers such as an H 2 molecule.The approach described in this study enables one to access not only the response functions and GFs but also various physical quantities on a quantum computer.Invention and enrichment of tools for such properties will enhance the practical use of quantum computers for material simulations., . . ., θ (1)  1•••1 .For an arbitrary input state |ψ and the linear combination O of unitaries in Eq. ( 1), O|ψ up to a normalization constant can be prepared by setting the parameters to appropriate values, as demonstrated below.We can assume that the coefficients in Eq. ( 1) are positive real values without loss of generality since each of the phase factors can be absorbed into the unitary coupled to it.
We use the notation | n ≡ |0 ⊗n ⊗ |ψ , H n ≡ H ⊗n ⊗ I, and R y (−2θ )|q ≡ |q, θ for q = 0, 1.The action of circuit on the initial state | n can be tracked by referring to the recursive definition in Fig. 2. Specifically, we find For each bit string k of length n, U k appears only once on the right-hand side of Eq. (A1).If the outcome of a projective measurement on the n ancillae is |0 ⊗n , the whole state collapses immediately after the measurement as follows: We derived the expression above by using |0, θ = cos θ |0 − sin θ |1 and |1, θ = sin θ |0 + cos θ |1 .
In order for the resultant state to be proportional to O|ψ , we first determine the 2 n−1 parameters θ (1)  0•••00 , . . ., θ (1)  1•••11 (each having a bit string of length n − 1) such that tan θ (1)  0 All the remaining parameters can also be determined in this way, so that the circuit C O allows one to prepare the desired state O|ψ probabilistically.It is noted that the probability ψ|O † O|ψ is determined not by the artifact (the number n of ancillae the user has introduced), but by the physics (the quantum state O|ψ itself) regardless of n.This fact also means that the probability of "success" in the state preparation gives the information about the normalization constant, which can be as important as the desired state.In this sense, measurement results other than |0 ⊗n are neither failures nor waste of time.
The construction procedure described above should be understood as a starting point for a generic nonunitary operator.There can be room for making the circuit more efficient by considering the characteristics of the operator and/or the states you want.For example, the circuit in Fig. 3 does not contain a rotation gate and only one of the two ancillae needs to be measured.It works thanks to the Fermi statistics.Another interesting point is seen in the circuit in Fig. 4, which can prepare the two important states n + pσ,p σ |ψ and n − pσ,p σ |ψ .[See Eqs. ( 27) and ( 28).]

Charge-charge contributions
The following procedures calculate the matrix elements in the manner described in Sec.II C.

a. Diagonal components
Procedure 2 Calculation of diagonal components of transition matrices.

Spin-spin contributions
The following procedures calculate the matrix elements in the manner described in Sec.II D.

Spin-charge contributions
The following procedures calculate the matrix elements in the manner described in Sec.II E. T ν λp j,p σ + = 1 8: T ν p j,p σ * = 1/N meas 9: return T ν p j,p σ

FIG. 10 .
FIG. 10.Schematic illustration of RHF orbitals and their electronic occupancies for (a) a C 2 molecule and (b) an N 2 molecule.The descriptions beside the energy levels represent the orbital characters.Those with asterisks are the antibonding orbitals.The 1σ and 1σ * MOs, coming from the 1s AOs of the constituent atoms, are not shown.The 1π and 1π * MOs come mainly from the π bonding of 2p AOs, while the 3σ and 3σ * MOs from the σ bonding of 2p AOs.The nondegenerate many-electron ground state for each system consists mainly of the same electronic configuration as the RHF one.

χ
FIG. 11.The response function (a) χ FCI1π n,1π n and (b) χ FCI 1π z,3σ z exact within the FCI solution for a C 2 molecule.The simulated ones χ FCI−stat for for the numbers of measurements N meas = 10 000 and 40 000 are also shown.

χFIG. 14 .
FIG. 13.The response function (a) χ FCI 3σ n,3σ n and (b) χ FCI 3σ n,1π * n exact within the FCI solution for an N 2 molecule.The simulated ones χ FCI−stat for the numbers of measurements N meas = 10 000 and 40 000 are also shown.

11 FIG. 15 .
FIG. 15. Circuit C O for n = 2 as a special case of that in Fig. 1.
Generic diagonal circuit C m 1 m 2 for probabilistic preparation of a † m 1 a m 2 |ψ and other states from an arbitrary input state |ψ and two ancillary qubits.We define the partial circuit U (m 1 m 2 ) by enclosing it with dashed lines.
Procedure 3 Calculation of off-diagonal components of transition matrices.