A Comparison of Three Ways to Measure Time-Dependent Densities With Quantum Simulators

Quantum algorithms are touted as a way around some classically intractable problems such as the simulation of quantum mechanics. At the end of all quantum algorithms is a quantum measurement whereby classical data is extracted and utilized. In fact, many of the modern hybrid-classical approaches are essentially quantum measurements of states with short quantum circuit descriptions. Here, we compare and examine three methods of extracting the time-dependent one-particle probability density from a quantum simulation: direct Z-measurement, Bayesian phase estimation, and harmonic inversion. We have tested these methods in the context of the potential inversion problem of time-dependent density functional theory. Our test results suggest that direct measurement is the preferable method. We also highlight areas where the other two methods may be useful and report on tests using Rigetti's quantum virtual device. This study provides a starting point for imminent applications of quantum computing.


INTRODUCTION
The real time simulation of quantum systems on a classical computer is a difficult problem even for a supercomputer due to the fact that the Hilbert space grows exponentially with the system size [1]. A universal quantum computer is believed to be the solution of the difficulty, where it is known that a wide class of physical systems can be simulated efficiently on a quantum computer [1][2][3][4][5][6]. But running a practically meaningful quantum algorithm may require a large amount of qubits, e.g., factoring 2,048 bit RSA integers may take up to 20 millions qubits [7], which is far beyond the capacity of the current best 53-qubit quantum computer [8]. So the current quantum technology works best when paired with classical algorithms. We have been studying the application of such a hybrid algorithm in quantum chemistry. The primary example is the time-dependent density functional theory (TDDFT) [9]. To utilize quantum technology in classical algorithms, quantum measurement are necessary, here we measured the density operator on Rigetti's quantum device and then utilized the density to perform the potential inversion within the framework of TDDFT.
Density functional theory (DFT) is a powerful tool in modeling condensed matter systems [10]. In the framework of DFT, a non-interacting system with a self-consistently determined potential is constructed to replace the interacting system. The additional potential term in the non-interacting system is known as the Kohn-Sham potential. Such a system with non-interacting particles is the Kohn-Sham(K-S) system. In the K-S system, the calculation of an exchange-correlation term is required. However, the exact form of the exchange-correlation potential is not yet known. This term is usually obtained with some approximation methods [11][12][13][14], machine learning methods [15,16]. In fact, the utilization of a quantum computer can help generate an accurate exchange-correlation potential. This idea is mentioned in the article [17], where a hybrid method of generating exchangecorrelation potential for classical DFT calculation is proposed.
The time dependent counterpart of DFT, time-dependent density functional theory (TDDFT) is widely used in finding the dynamics of the system when a time-dependent potential is present. Similar to DFT, TDDFT uses the time dependent K-S system where a time-dependent K-S potential is required. We call the task of constructing such a K-S potential when given the timeevolution of the on-site probability density, the K-S potential inversion problem. In article [18], a scheme of solving the K-S potential inversion problem utilizing a quantum computer was proposed. We have recently returned to this proposal with improved numerical methods for inverting the potential [19]. To obtain the K-S potential, we need to get the density of the time evolved many-particle system using a quantum computer.
In this paper, we will present three different methods of measuring the density operator on a quantum computer and compare the performance of the methods.
An outline for the remainder of the article is as follows: first, we discuss the phase estimation approach to measurement. Then we describe the circuit implementation for measuring the on-site fermionic density. Qubit descriptions for the fermionic operator are explained in the next part followed by the illustration of a twoelectron test. Finally, three schemes for extracting the density are tested numerically and compared.

Phase Estimation
Quantum phase estimation [20,21] plays an important role in the quantum algorithm zoo [22], it is a key sub module of many quantum algorithms [23][24][25]. It is also an important procedure to measure the on-site density operator in our work.
We'll next describe the general picture of doing the measurement of an arbitrary operator and how quantum phase estimation plays a role in our work. To implement the measurement of an arbitrary observable, we will consider the circuit as shown in Figure 1. The circuit has two parts, the part before the dashed line is for evolving the initial state at time t under a fixed fermionic Hamiltonian of chemical interests.
The system of the most chemical interests is the interacting electron system. The Hamiltonian of a many-body interacting system is given by where V ext (r i ) is the external potential energy consists of the interaction between the electrons and the external field. The second quantized form of the above many-body Hamiltonian is given by where the fermionic operators {a p , a † p } satisfy a † q a p + a p a † q = δ pq , a p a q = −a q a p and a † p a † q = −a † q a † p . Given the basis set {χ p (r)}, the coefficients h pq , h pqrs are given by The latter half is a phase estimation circuit where U O (τ ) = e −iOτ where O is the observable to be measured. For a general state |ψ , we can expand it in the eigenspace of the operator O. To be precise, for a general state |ψ = k c k |k , where O k and |k are the eigenvalue and eigenvector of the operator O. Thus, the probability of measuring zero on the top register is given by where c k (t) = k|U(t)|ψ . In this article, we only consider O = n j = a † j a j in order to measure the local on-site density at site j. This is because the inverse potential is determined by the on-site density and its first and second order derivatives [18,19].
The eigenvalues of n j = a † j a j are 0 and 1, so the wave function after the unitary evolution U(t) is given by |ψ(t) = c 0 (t) |ψ n j =0 + c 1 (t) |ψ n j =1 . Thus, the expectation value of the density is given by,

Qubit Encoding
To implement the evolution and phase estimation algorithm on a quantum computer, we need to encode the Hamiltonian into qubits. A standard way is to use Jordan-Wigner (JW) transformation, which encodes a fermionic system of M orbitals into M qubits.
With the above transformation, the fermionic Hamiltonian can be encoded into qubit representation. Thus the Hamiltonian can be written as H = i h i , where all the h i 's are tensor product of Pauli operators. There is not an easy way to construct arbitrary unitary operators on a quantum computer [26]. A pragmatic way to simulate the propagator U(t) = e −iHt is applying the Trotter decomposition.
Each individual term in the decomposition above can be simulated efficiently on a quantum computer [27].

RESULTS AND DISCUSSION
We tested our methods on two different Hamiltonians, one is the 4-orbital HeH + model, the other is the 8-orbital HeH + . For the 4-orbital model, the basis set used to examine the HeH + molecule is that given in [28] which results in four spin orbitals. The interatomic distance is 1.401 Bohr. The basis functions are orthogonalized and then transformed such that the onebody Coulomb matrix is diagonal. This transformation was chosen so that a corresponding scalar time-dependent Kohn-Sham potential could be calculated using for this system using the method of [19].
For the larger case, we used HeH + at the same geometry but in the 6-31G basis set [29]. This resulted in twice the number of basis functions as the minimal example. The integrals in the 6-31 basis were computed using Pis4 [30].
The basis functions are orthogonalized and then transformed such that the one-body Coulomb matrix is diagonal. This transformation was chosen so that a corresponding timedependent Kohn-Sham potential could be calculated using for this system using the method of [19]. In both models, the initial state at t = 0 places two electrons in the first two modes of opposite spin. This state is obtained by employing two X-gates to prepare |ψ(0) = |1100 in the 4-orbital model and |ψ(0) = |11000000 in the 8-orbital model.
Using Rigetti's quantum virtual machine [31], we then evolve the system under its Hamiltonian for times less than three atomic units. The propagation is implemented via the first-order Trotterization with Trotter step equal to three. To reduce the Trotter error in evolution, either a shorter Trotter step or a higher order Trotter approximation must be used [32]. This means more quantum gates are needed, making it harder to be implemented on a near term device. Additional sources of error are associated with finite sampling from the binomial distribution and the error associated with the inference steps. To make the virtual machine slightly closer to a real quantum computer, in all methods below, measurement noise was added into the system, giving 1% probability of flipping the qubit. It should be noted that the quantum noise found on the actual device was much higher, so we will not present the results from the actual quantum device.
In Figures 2A,C

Method 1: Z-Basis Measurement
In the first method, we rely on the fact that the Jordan-Wigner transformation of the on-site density operator has a simple form a † p a p = (1 − Z p )/2. Thus, we can directly measure the local density operator by measuring Z p without passing in the phase estimation circuit after the dotted line in Figure 1.
By repeating the measurement at each time step in the time range 0 ≤ t ≤ 3, we obtain the expectation value of the density. The results based on 15 equally spaced time-points with 3,000 measurements at each fixed time are shown in Figures 2A,D. The exact time evolution of the density is also shown in the figure for comparison along with error bars of 2σ reflective of the N = 3, 000 sample variance of the binomial distribution.
The simplicity of this measurement approach reduces the classical runtime to the lowest of the three methods compared, and the convergence of the error bars is faster than the Bayesian measurement discussed later.

Method 2: Harmonic Inversion
Harmonic inversion is a technique of extracting the amplitudes A j , frequencies f j , phases φ j , and exponential decay constants α j out of a signal, which is evenly sampled [33,34]. The signal reconstructed from harmonic inversion has the same form as the probability P(0|τ ) except for the decaying term which is negligible when the decoherence is not considered. By comparing the form of the reconstructed signal with the probability, we can obtain the density from the reconstructed signal.
The results of density measurement through harmonic inversion are shown in Figures 2B,E. Each point in Figures 2B,E was computed through harmonic inversion using the HarmInv package [35]. Because the local density operator a † p a p only has eigenvalues zero and one, the measurement outcome has a simple form where A 0 (t) = 1 2 (2 − |c 1 (t)| 2 ), A 1 (t) = |c 1 (t)| 2 /4, and f = 1/2π.

Method 3: Bayesian Inference
Bayesian inference can be used to estimate the density as well.
As a powerful tool of making inferences, Bayesian inference has wide applications. We applied Bayesian inference to infer the unknown parameters in a quantum system which, in our case, is the on-site density. The density estimation was implemented via sequential Monte Carlo (SMC) [36]. This method requires the most communication between the classical and quantum processors since the SMC suggests each τ -point based on the previous outcomes. The Bayesian experimental design is based on the implementation found in the QInfer package [37]. Bayesian inference gives the probability distribution of a parameter over the parameter space. The final decision is made according to the posterior probability P(θ |d 1 , d 2 . . . d N ), where θ is the parameter we want to estimate, d i 's are the outcome of each measurement.
In the present application, θ ≡ n j (t) .
Recall the Bayesian rule, the posterior probability is updated by carrying out experiments sequentially, where P(θ ) is the prior probability, P(d i |θ ) is the likelihood function.
The likelihood function contains the information about the parameters before conducting any experiments. Since we know nothing before the experiment, we can initialize the prior with a uniform distribution over the parameter space. For the phase estimation circuit of Figure 1, the likelihood function is given by where U O (τ ) = exp(−iτ a † j a j ) and d = 0 or 1. Note, when d = 0 we recover Equation (5).
With this we can rewrite the likelihood function as This can be compared with Equation (5) in the case that d = 0.
The results of Bayesian inference are shown in Figures 2C,F. Bayesian inference has good performance within a wide range of the time domain except at the boundary of the estimate domain e.g., when the density is one or zero. This is based on numerical evidence since the majority of the points at or near the boundary of the estimation domain needed to be discarded when cleaning the data as discussed below.
Unlike harmonic inversion, τ in the phase estimation circuit is not required to be evenly spaced. Another advantage of Bayesian inference is that we do not need to know the exact form of the function to be estimated a priori. Bayesian inference could also be applied to estimate more general parameters.

Comparison
To quantify the accuracy of these density extraction methods, we employ the L 1 norm to measure the deviation from the Trotter solution. For discrete data points, the deviation is given by the loss function on the density at the first site: L = N k |ñ 1 (t k ) − n 1 (t k ))|/N, whereñ 1 (t k ) is the outcome of the measurement at time t k , n 1 (t k ) is the solution of the Trotterized Hamiltonian. Figure 3 shows how the loss function scales with the number of trials for each of the three approaches. The convergence rate for determining the bias of a coin would be 0.5 but here additional measurement error has been introduced into the model which prevents L = 0 situation even with an infinite number of samples. Further, in our implementation, the Bayesian and harmonic inversion techniques sometimes reported anomalously poor estimates of the density at a given time. A single fluctuation of this type along the time trace of the density entirely dominates the loss function. For the sake of comparison, we did not include the data points that are 5σ away from the exact solution in all three methods. This led to more stable results when the number of trials is small. Another benefit of filtering the data is that for the Bayesian inference, estimates close to the boundary of the domain are subject to large fluctuations giving poor estimates. So we can exclude the wrong data points by setting a 5σ window. Although the discarding procedure is ad hoc and requires knowing the exact answer, we have tested our data at various levels of cutoff finding that at any fixed cutoff harmonic inversion had the most points discarded and consistently displayed marginally faster convergence rates. Figure 3A shows the scaling of loss function of the 4-orbital model. The slope of the fitting lines are −0.4961, −0.3747, and −0.4103, respectively. Points 5σ away from the exact density under the Trotter approximation are not used for calculating the loss function. This resulted in 73.87, 87.26, and 90.07% of points used in the plotted data, respectively. Figure 3B shows the scaling of loss function of the 8-orbital model. The slope of the fitting lines are −0.4346, −0.4595, and −0.4184, respectively. Points 5σ away from the exact density under the Trotter approximation are not used for calculating the loss function. This resulted in 80.53, 84.40, and 93.20% of points used in the plotted data, respectively.
Harmonic inversion measures 40 times more than the other two methods, so the actual data and fitting line should be shifted to the right by 40 times the number of measurements showing in the figure.
Regardless of the possible improvement in convergence, it should be reminded that the harmonic inversion technique uses many quantum computer queries to estimate P(0|τ , t) at variable τ before inferring the density at a fixed time t. In comparing the three methods, all require time evolution of the system wave function to time t. In the harmonic inversion and Bayesian estimation techniques, additional gates are needed for the τ propagation under the observable for density. The difference between queries in harmonic inversion and Bayesian inference is the selection of the τ parameter in U O (τ ).
While the convergence rates are all approximately the same, it is clear that the Z-basis measurement has the best performance in terms of the number of queries of the quantum computer. In the case considered here, the direct Z measurements are convenient for the Jordan-Wigner encoding. In other circumstances with different fermion-to-spin transforms, the direct measurement technique may not be as fruitful. For existing and near-term quantum devices, the constraints of low circuit depth suggests direct measurement of the Z operators as the best path forward when using a Jordan-Winger transformed qubit Hamiltonian.
The runtime of these three methods also varies. Since direct Zmeasurements are the simplest from an inference point of view, the classical computation time is also the least. Bayesian inference requires many steps for the sequential Monte Carlo to converge [36]. Consequently, this method used the longest amount of classical computational time. Although harmonic inversion uses 40 times more measurement per time-point, it is interesting to note that it only took an intermediate amount of classical processing time.

CONCLUSIONS
We have tested three different methods of measuring the on-site density operator for a toy model inspired by TDDFT. We were able to conclude that direct Z measurements obtains the best estimates of the on-site density for a given number of quantum computer queries. This is based on the use of the Jordan-Wigner transform and simulated measurement noise. Of course, we could have considered other fermion-to-spin transforms which lead to different encodings of the a † i a i . For improving our noise models, we can do no better than testing our circuits on current and future quantum devices. We tested our circuits on Rigetti's quantum device but found that the loss function depends heavily on which qubits are used as well as the permutation of qubit labels within the circuit. Time evolution under the full Hamiltonian did not return any signal even when using only one first-order Trotter step. We therefore resorted to using a truncated Hamiltonian which included the one-body Hamiltonian and only the Coulomb-like (h ijji ) terms of the two-body Hamiltonian. After encoding and exponentiation, this Hamiltonian results in 66 universal gates and compiled non-deterministically using the PyQuil package PyQuil to approximately 200 allowable gates on the Rigetti device. Due to decoherence, only a weak signal was present where amplitudes recovered were between three and twenty percent of the exact solution. The recovered amplitude depended mostly on qubit selection but also changed run-to-run. The frequency and sinusoidal shape of the signal was recovered more reliably. In our present study, the eigenenergies were not interesting but we suspect that problems that depend on the frequencies may be more successfully calculated on the current Rigetti device.
We plan to continue our inquiry into the TDDFT potential inversion problem using existing and forthcoming quantum technology. Tasks that avoid QMA-hard [38] state preparation problems will continue to be of interest to those looking for new applied areas of quantum computation.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.