Probabilistic imaginary-time evolution by using forward and backward real-time evolution with a single ancilla: first-quantized eigensolver of quantum chemistry for ground states

Imaginary-time evolution (ITE) on a quantum computer is a promising formalism for obtaining the ground state of a quantum system. As a kind of it, the probabilistic ITE (PITE) takes advantage of measurements to implement the nonunitary operations. We propose a new approach of PITE which requires only a single ancillary qubit. Under a practical approximation, the circuit is constructed from the forward and backward real-time evolution (RTE) gates as black boxes, generated by the original many-qubit Hamiltonian. All the efficient unitary quantum algorithms for RTE proposed so far and those in the future can thus be transferred to ITE exactly as they are. Our approach can also be used for obtaining the Gibbs state at a finite temperature and the partition function. We apply the approach to several systems as illustrative examples to see its validity. We also discuss the application of our approach to quantum chemistry by focusing on the scaling of computational cost, leading to a novel framework denoted by first-quantized quantum eigensolver.


I. INTRODUCTION
When solving a practical computational problem for which quantum parallelism can outperform classical algorithms, e.g., quantum chemistry [1] and combination optimization [2,3], the problem is often rewritten according to the quantum dynamics governed by an appropriately defined many-qubit Hamiltonian. The problem is, at least computationally, reduced to finding the ground state of the target system. Imaginary-time evolution (ITE) on classical computers is a technique that is widely used to obtain the ground state of a quantum system. It is based on the imaginary-time Schrödinger equation [4], derived by replacing real time t with a purely imaginary value −iτ , wherein the real parameter τ is called the imaginary time. A wave function governed by this equation exhibits dynamics wherein excited states decay rapidly as the imaginary time proceeds.
The second type, QITE, solves the linear equation on a classical computer to determine the coefficients of Pauli tensors for unitary evolution. The equation is derived such that a unitarily evolved state approximates the exactly (nonunitarily) evolved state in the sense of the squared norm; this derivation is mathematically equivalent to McLachlan's variational principle for VITE. The matrix dimension for the linear equation tends to grow rapidly with an increase in the problem size, i.e., the number n of qubits involved in a given system and that of Pauli tensors for the approximate unitary evolution. More precisely, when one decides to expand the approximate unitary operator in terms of Pauli tensors of size D, called the domain size, there exist 4 D tensors (including the identity) for each combination of D qubits picked from n. [8] Naive adoption of such n C D combinations leads to a linear equation of dimension d = 4 D n C D . Although solving this equation by using conventional methods on a classical computer requires polynomial computational time with respect to d, the time grows rapidly as D and/or n increase. Its scaling can be regularized if certain approximations that consider the locality of partial Hamiltonians are adopted.
[21] A QITE circuit can be compressed in several ways. [21,22] The third type, PITE, exploits measurements to implement nonunitary operations. An input state undergoes the unitary operations, and the measurement is performed; then, the state collapses to the desired state with a certain probability. Liu et al. [15] analyzed their PITE circuit for the Grover's search algorithm [23,24], in which they introduced an additional qubit per pair of qubits for a nonunitary operation on the pair. The measurement-based generation of desired states is used for other purposes such as linear equations [25], Green's functions [26,27], linear-response functions [28], and ionization energies. [29] The three types of ITE summarized above are not necessarily exclusive to each other. In fact, a combination of QITE and PITE was proposed recently. [30] Although the implementation of PITE by using only one ancillary qubit is mathematically ensured to be possible[13, 28,31], it requires classical computational resources for the singular-value decomposition of the total many-qubit Hamiltonian or those for other methods having the same cost, which is typically estimated to be O((2 n ) 3 ) for an n-qubit system. In this study, we propose a generic nonvariational approach of PITE that uses only one ancillary qubit and that does not demand explicitly large classical cost for linear algebra. Within the firstorder of the time step, our approach enables one to construct the circuit from forward and backward real-time evolution (RTE) generated by the Hamiltonian without knowing its eigenvalues and eigenvectors. Thus, every unitary quantum algorithm for real-time dynamics can be transferred to the imaginary-time dynamics without modifications. Further, our approach can be used to obtain the Gibbs state at a finite temperature and partition function. We apply the approach to several systems to confirm its validity. In addition, we discuss the overall PITE procedure of quantum chemistry calculations for obtaining the ground states of molecules by focusing on the scaling of computational cost, which leads to a novel framework referred to as a first-quantized eigensolver (FQE).

II. METHODS
A. PITE for generic cases

Exact circuit
Let H be the Hamiltonian for an n-qubit system. The ITE operator for an imaginary-time step ∆τ is e −H∆τ . For an input state |ψ , we want to get the evolved state e −H∆τ |ψ up to a normalization constant. From the evolution operator and an adjustable real parameter m 0 , we define a nonunitary Hermitian operator M ≡ m 0 e −H∆τ , for which we impose conditions 0 < m 0 < 1 and m 0 = 1/ √ 2. These conditions allow us perform Taylor expansion safely to obtain Eq. (4) later. To establish the exact PITE, we define the following Hermitian operator formally for the n-qubit system: This definition with κ ≡ sgn The operators e ±iκΘ are unitary and hence they can be implemented as quantum gates for the n qubits.
We introduce an ancillary qubit to construct the circuit C PITE for the (n + 1)-qubit system, as shown in Fig. 1(a). The circuit uses the single-qubit gate The composite system undergoes the unitary operations as where the normalization constant is correctly taken into account.

Approximate circuit
It is known [31] that an arbitrary diagonalizable complex matrix can be written as a linear combination of at most two unitary matrices provided that the resources of classical computation for its singular-value decomposition is available. The probabilistic implementation of a linear combination of untaries [28] together with this fact ensures mathematically the possibility of probabilistic imaginary-time dynamics with a single ancillary qubit (see also Ref. [13]). Since the exact PITE described above involves no approximation and is applicable to an arbitrary ∆τ , it is mathematically desirable. The adoption of it, however, is hardly of practical use due to the operator Θ. Specifically, it is formidable to implement the exponentiated Θ as a sequence of elementary quantum gates for the multiple qubits without knowing its eigenvalues and eigenvectors. Therefore we have to circumvent such a difficulty.
As usual, we assume that the imaginary-time step is small. We can Taylor expand Θ, defined by Eq. (1), in ∆τ for such a case as κΘ = θ 0 −Hs 1 ∆τ +O(∆τ 2 ), where θ 0 ≡ κ arccos[(m 0 + 1 − m 2 0 )/ √ 2] and s 1 ≡ m 0 / 1 − m 2 0 . We can rewrite the anti-controlled e iκΘ and the subsequent controlled e −iκΘ in the exact circuit in Fig. 1(a) as where I 2 n is the identity for an n-qubit state and R z (θ) is the single-qubit z rotation. U RTE (∆t) ≡ e −iH∆t is the usual RTE operator with respect to the original Hamiltonian for a real-time step ∆t. By using Eq. (4), we can construct the circuit C PITE as shown in Fig. 1(b), which implements approximately the exact one, C PITE , within the first order of ∆τ . In fact, it is easily confirmed that the state of the (n + 1)-qubit system immediately before the measurement is consistent with the RHS of Eq. (3) within the first order of ∆τ . The forward U RTE (s 1 ∆τ ) and backward U RTE (s 1 ∆τ ) † RTE operators for the rescaled time step in the circuit are responsible for the imaginary-time dynamics. The approximate circuit found here is the main result of the present work. The circuit C PITE is suggestive of sophisticated development of quantum algorithms for ground states of quantum systems. The direct operations on the input register come solely as the RTE generated by the Hamiltonian, which means that we have pushed the intricacy about the nonunitary operation e −H∆τ into the two black boxes, or equivalently an oracle, in which the tractable unitary operations of the form e ±iH∆t are implemented. Every efficient unitary algorithm for the real-time dynamics can thus be transferred to the imaginary-time dynamics exactly as it is. This availability is favorable since the realtime dynamics in quantum computation is much more tractable than the imaginary-time dynamics in general.
Since the accuracy of evolution is worse for a larger ∆τ , we need to apply the approximate circuit for a small ∆τ as a Trotter step and perform the measurements repeatedly. We have to get the success states for all the measurements to reach the final state where all the excited states have diminished sufficiently. It is easily understood for an n steps -step evolution that the probability for obtaining the success state having survived all the n steps measurements is ψ|M 2nsteps |ψ , exponentially decaying as the measurement number increases.
The construction of our single-ancilla PITE circuit is possible primarily due to the Hermiticity of the generator for evolution. For the dynamics generated by a non-Hermitian operator, the probabilistic evolution is also possible if two ancillary qubits are available. For details, see Appendix B.

Comparison with other ITE approaches
The advantage and disadvantage of PITE compared with VITE and QITE have to be mentioned. First of all, PITE does not involve arbitrariness for the implementation of state evolution of qubits on a circuit. That is, PITE does not need an ansatz circuit unlike VITE and it does not need selection of Pauli tensors according to the locality of partial Hamiltonians unlike QITE. The expected energy in a PITE calculation is thus minimized in the full Hilbert space of qubits which encode the state of target system. On the other hand, the exponential decrease in success probability during the PITE steps is a deplorable drawback inherent to the PITE approach as long as we are not interested in the normalization constant. A promising alternative for alleviating this difficulty is the quantum amplitude amplification (QAA) technique, [32,33] proposed by Brassard and his coworkers, known as a generalization of the Grover's search algorithm. [23,24] In the original QAA, [33] the rotation angle within the subspace is determined by the quantum amplitude estimation (QAE), which is a kind of quantum phase estimation (QPE). [17] It is possible, however, to perform QAA without QAE by employing the maximum likelihood estimation, as demonstrated by Suzuki et al. [34] The examinations of PITE in combination with QAA will be interesting and important for practical applications of the PITE approach where the input state has only a small overlap with the ground state.

B. FQE framework for quantum chemistry
One of the most practical applications for finding the ground states of quantum systems is quantum chemistry calculations for molecules, where many electrons interact with each other in three-dimensional space. As a promising application of our PITE approach, we describe here the framework for obtaining the ground state wave function of a molecular system. This framework, which we call FQE, is based on the scheme proposed by Kassal et al. [35] for efficient implementation of real-time dynamics for quantum chemistry.

Single particle in one-dimensional space
We begin with the consideration of a single particle in one-dimensional space as the first step for establishing the FQE framework. We assume that the particle with a mass m moves within a range [0, L] on the x axis, whose first-quantized Hamiltonian is H =T +V (x). T ≡p 2 /(2m) is the kinetic-energy operator for the momentum operatorp. V is the external potential for the position operatorx. We assume that n qubits are provided us for representing the wave function in real space.
We want to implement the RTE operator e −iH∆t for a real-time step ∆t that acts on the qubits whose state encodes the wave function of the particle. We discretize the finite interval to define N ≡ 2 n equidistant points, that is, x (k) ≡ k∆x (k = 0, . . . , N −1) for the step ∆x ≡ L/N. We can thus assign each computational basis |k n to each of the position eigenstates:x|k n = x (k) |k n . The coefficient of each computational basis in an n-qubit state thus represents the wave function ψ(x) of the particle as As the canonical counterpart of the discretized positions, we define the N discrete momenta p ( s) ≡ s∆p ( s = −N/2, −N/2 + 1 . . . , N/2 − 1) for the step ∆p ≡ 2π/L in reciprocal space. Recalling that the quantum field theory [36,37] relates the creation and annihilation operators of momentum eigenstates with those of position eigenstates via Fourier transform, we should define in this case the momentum eigenstate as We have to notice, however, that these eigenstates of position and momentum do not transform between each other via the ordinary quantum Fourier transform (QFT). [17] That is, the minimum of position is 0, while that of momentum is −N ∆p/2. This discrepancy is taken into account appropriately by employing the centered QFT (CQFT) [38,39], which is defined as the application of QFT after σ x on the qubit corresponding to the highest bit for the index of an n-qubit computational basis: CQFT = QFT (σ x ⊗ I 2 n−1 ) . One can easily confirm that CQFT transforms the position eigenstate to the momentum eigenstate correctly: CQFT|k n = |p ( k) , where k ≡ k − N/2. The meaning of the tilde symbol for an integer is the same as this in what follows.
|p ( s) is also the eigenstate ofT belonging to the discrete kinetic energy E s ≡ s 2 (∆p) 2 /(2m). If we have a kinetic-phase gate U kin which acts on the computational basis diagonally as U kin (∆t)|j n ≡ exp(−iE j ∆t)|j n , it becomes under CQFT as where we used the completeness of the momentum eigenstates to get the last equality. Equation (7) means that the evolution coming from the kinetic operator can be implemented by using CQFT and the kinetic-phase gate. [39] For the details of the action of kinetic-evolution operator on the position eigenstate, see Appendix C.
As for the potential-evolution operator e −iV (x)∆t , we need a potential-phase gate U pot which acts on the position eigenstate diagonally as U pot (∆t)|k n ≡ exp(−iV (x (k) )∆t)|k n . There exist several alternatives for implementation of the kinetic-and potential-phase gates. While Kassal et al. in the original paper [35] adopted the phase kick-back technique [40], Ollitrault et al. [39] adopted the technique proposed by Benenti and Strini [41], which introduces to a computational basis the phase factor given as a polynomial of position. The extension of the latter technique for a piecewisely defined polynomial is possible [39] in combination with efficient comparators [42], which decide the interval the index of an input computational basis falls onto.
The kinetic energy and potential operators do not commute with each other in most systems. Therefore we have to consider application of them to the qubits sequentially by employing the first-order Suzuki-Trotter as e −iH∆t = e −iT ∆t e −iV (x)∆t + O(∆t 2 ) for the RTE. Given this approximation, we adopt the circuit C (1) PITE in Fig. 1(b) for PITE. Since the RTE operator in this case is U RTE (s 1 ∆τ ) = CQFT · U kin (s 1 ∆τ ) · CQFT † · U pot (s 1 ∆τ ), the neighboring CQFT and CQFT † in the forward and backward evolution fortunately cancel each other, leading to the circuit C (ST1) PITE shown in Fig. 2(a). The number of operations for each imaginary-time step is [17] c kin and c pot are those of the kineticand potential-phase gates, respectively. α and α ′ are constants independent of n, coming from the introduction of control bits. If we adopt the implementation of phase gates in Refs. [39,41], the operation number of the kinetic-phase gate is c kin = O(n 2 ). Similarly, if the potential is given as a single or a piecewisely defined pthorder polynomial, the operation number of the potentialphase gate is c pot = O(n p ). The operation number for the PITE step is thus c PITE = O(max(n 2 , n p )).
The analyses of operation numbers in the circuits provided in this paper should be mentioned here. If one aims to discuss strictly the operation number of a given circuit, an elementary gate set has to be defined so that each member is the unit for counting the operations. A typical elementary gate set consists of a small number of single-qubit gates and a CNOT gate. [17] In the present study, however, we do not need to introduce any specific elementary gate set for discussing the scaling of operation numbers with respect to the numbers of grid points and particles. It is because every many-qubit operation appearing in this study is a controlled single-qubit operation for which the number of control qubits is independent of the numbers of grid points and particles. Specifically, each QFT gate is implemented by using singly controlled rotations, CNOTs, and SWAPs [17]. The depth of each of these building blocks is independent of the grid points and the particle number. This is similarly the case for the phase gates comprising the entire PITE circuit since we implement them with the techniques in Refs. [39,41]

Interacting electrons in three-dimensional space
Let us move on to establishing the FQE framework. We discuss here the encoding of a many-electron wave function and the construction of PITE circuit for a quantum chemistry problem, based on which we estimate the scaling of computational cost as a function of problem size. The discussion below defines and characterizes the FQE framework. We ignore the spin degree of freedom for simplicity in what follows.
The extension of the scheme for the case of a single particle in one-dimensional space explained above to the case of interacting electrons in three-dimensional space is straightforward. Let us consider n el electrons confined to a cube as a simulation cell, each of whose edges has a length L. We assume that n d qubits are available for the degree of freedom in each direction for each electron. We discretize the simulation cell to define a regular mesh which has N d ≡ 2 n d points in each direction. The wave function for an electron can thus be encoded in 3n d qubits, whose coefficients of computational basis represent the wave function at N grid ≡ N 3 d equidistant grid points. The n el -electron wave function Ψ(r 1 , . . . , r n el ) is thus encoded in 3n d n el qubits as a linear combination of |k 1x , k 1y , k 1z ⊗ · · · ⊗ |k n el x , k n el y , k n el z , where |k µx , k µy , k µz is the position eigenstate of µth electron specified by the three integers. The required memory for storing the wave function on a classical computer is O((N grid ) n el ), scaling exponentially with respect to n el and n d . On the other hand, the required number of qubits for the storage is O(n el log N grid ), exhibiting polynomial scaling, that is the fundamental advantage of quantum simulation discussed by Feynman. [43] The extension of our formalism described below to a case in which the simulation cell is spanned by arbitrary three independent vectors will be straightforward. In addition, the numbers of qubits for encoding the wave functions can differ from each other for the three directions.
The first-quantized Hamiltonian for electrons in-volves typically interactions of the formV int = µ>µ ′ v int (r µ ,r µ ′ ), where v int is a pairwise interaction such as the Coulomb interaction. The implementation of RTE thus requires the interaction-phase gate U int that acts diagonally on the position eigenstate for a pair of electrons as U int (∆t)|r ⊗ |r ′ ≡ exp(−iv int (r, r ′ )∆t)|r ⊗ |r ′ . If v int is given as a simple or a piecewisely defined p ′ th-order polynomial, the evolution operator exp(−iV int ∆t) can be implemented with the operation number c int = O(n 2 el n p ′ d ) by using the techniques in Refs. [39,41] The whole FQE procedure for the interacting electrons is shown in Fig. 2(b). U ref prepares the reference state as the input to the first PITE step. We have to keep in mind that the reference state has to be generated carefully so that it is antisymmetric under arbitrary exchange of the electrons. [44] Berry et al. [42] proposed an efficient technique for antisymmetrization of a many-electron wave function. Such a process is needed in general when we work with the first-quantized formalism, in contrast to the second-quantized one. The kinetic-evolution operators for each electron can be implemented for the three directions independently, as shown in the figure. The operation numbers for the individual components in the PITE step are also shown. With possibly p ′ ≥ 2 in practice, the largest number c int of operations comes from U int [see Fig. 2(b)], which dominates the overall scaling of the computational cost c PITE of the whole PITE step for increase in the electron number.
The accuracy of quantum dynamics depends on the resolution of the wave function, or equivalently the step ∆x between neighboring grid points in each direction. Since the number of electrons is roughly proportional to the volume of simulation cell in a typical calculation, we have, for a fixed grid step, n d ∝ log(L/∆x) ∝ log(n 1/3 el /∆x). The required number of qubits for the many-electron wave function is thus estimated to be O(n el log(n 1/3 el /∆x)). The operation number in the whole PITE step in this case is c PITE = O(n 2 el (log(n 1/3 el /∆x)) p ′ ). If we adopted the second-quantized formalism like VQE instead of the first-quantized one, we had to work with the creation and annihilation operators of electrons for the molecular orbitals and hence the circuit in Fig. 2(b) could not be used. We have to instead construct the circuit in Fig. 1(b) naively for the second-quantized Hamiltonian. For n orbs localized orbitals at each atom in a target molecule, O(n orbs n el ) qubits are required for representing the many-electron state. Due to the numerous two-electron integrals in the Hamiltonian [45], c PITE in this case is at least O((n orbs n el ) 4 ). The firstquantized formalism is thus more favorable than the second-quantized one in order for a real quantum computer to conform to the coherence time. Furthermore, the classical computation of the two-electron integrals between molecular-orbitals is needed with cost O((n orbs n el ) 5 ) before the quantum computation starts.
The FQE framework is free from such an overhead and it is hence a promising candidate for quantum chemistry calculation of the ground states of large molecules.

C. PITE for finite-temperature states
The thermodynamic property of a system subject to a heat bath is described by the partition function Z(β) for an inverse temperature β. The thermal equilibrium is attained so that the Helmholtz free energy F is minimized. It is related with the partition function as e −βF (β) = Z(β). Various approaches for treating finitetemperature states on a quantum computer have been proposed. [7,8,12,[46][47][48][49][50][51][52][53] We demonstrate here that our PITE approach can also generate the Gibbs state and calculate the partition function.
The partition function is defined as Z(β) ≡ Tr sys exp(−βH), where the trace is for the n-qubit target system having the Hamiltonian H. We want to prepare the Gibbs state e −βH /Z(β). To this end, we use the maximally entangled 2n-qubit state generated by the U ME gate defined in Fig. 3(a). This gate entangles the input n pairs of initialized qubits to make them the Bell pairs. [17] We construct the circuit C Gibbs for probabilistic preparation of the Gibbs state, as shown in Fig. 3(b). This circuit contains U ME as a part of it to generate the maximally entangled state between the target system and the environment composed of the same number of qubits. We set the imaginary-time step for C PITE to ∆τ ≡ β/2.

The composite system undergoes the unitary operations as
which is the state immediately before the measurement on the ancilla. |j n (j − 0, . . . , 2 n − 1) is the computational basis of an n-qubit system. The reduced density operator [17] of the (n + 1)-qubit system consisting of the target system and the ancilla is obtained by tracing out the environmental qubits as The probability for obtaining |0 from the measurement is thus where P 0 is the projection operator onto the subspace spanned by the ancillary |0 state. The (n + 1)-qubit state for this outcome will be P 0 ρP 0 /P 0 = e −βH /Z(β) ⊗ |0 0|. This state is nothing but the direct product of the normalized Gibbs state and the ancillary |0 state. Since the partition function is defined as the summation of the Boltzmann factors over all the possible states of n qubits, the 2 n terms in it gives Z(β) = O(2 n ), leading to P 0 = O(1). [See Eq. (10)] This fact ensures that, for example, increasing n for the quantum simulation of a given real molecule to achieve higher resolution of the spatial discretization will converge to some nonzero P 0 unique to the molecule in the limit of n → ∞. It is clear that the partition function can be obtained simply from the ratio of the number of outcomes for the success state to the number of measurements, which immediately gives the free energy via the relation e −βF (β) = Z(β). In this sense, losing in this lottery is not waste of time. If we did not have the partition function, the free energy had to be calculated in another way, e.g. the thermodynamic relation F = E−S/β, where E is the internal energy and S is the von Neumann entropy. The calculation of S itself within the framework of PITE seems difficult, similarly to the case of VITE. [53] It is noted here that there exist quantum algorithms for obtaining the complex partition functions for analyses of quantum critical phenomena. [54][55][56] III. APPLICATIONS A. Two-level system Let us consider a two-level system for which we assign its ground state and excited state to single-qubit states |0 and |1 , respectively. The Hamiltonian of this system is H = ε gs |0 0| + ε ex |1 1|, where ε gs and ε ex are the energy eigenvalue of ground state and that of excited state, respectively. We express the input state for the kth PITE step (k = 0, 1, . . . ) as |ψ k = cos φ k |0 +sin φ k |1 by using a real mixing angle φ k . The operator M acts in this case as M|ψ k = m 0 (e −εgs∆τ cos φ k |0 + e −εex∆τ sin φ k |1 ).

For exact circuit
We first examine the case where the exact circuit C PITE in Fig. 1(a) is used for the PITE steps. By using (λ = gs, ex), the operator defined in Eq. (1) for the present system is written as Θ = θ gs |0 0| + θ ex |1 1|. The unitaries generated by this Hermitian operator are thus e ±iκΘ = e ±iκθ R z (κ∆θ), PITE as a special case of C (1) PITE in Fig. 1(b). This circuit implements PITE for a particle in one-dimensional real space within the first-order Suzuki-Trotter. U λ ≡ U λ (s1∆τ ) (λ = kin, pot), defined in the main text, is used in this figure. where θ ≡ (θ ex + θ gs )/2 and ∆θ ≡ θ ex − θ gs . They constitute the circuit C PITE , which is specifically depicted in Fig. 4(a). The normalized success state is obtained with the probability where ∆ε ≡ ε ex − ε gs is the excitation energy and w k ≡ tan 2 φ k is the relative weight of the excited state in the input state. It is clear from Eq. (13) that the mixing angle and the relative weight change as φ k −→ φ k+1 = arctan(α tan φ k ) and w k −→ w k+1 = α 2 w k , respectively, where α ≡ e −∆ε∆τ is the decay factor. We can then obtain w k = α 2k w 0 , indicating that the relative weight decays exponentially as the PITE procedure continues. The success probability in Eq. (14) is rewritten as a closed form p k = m 2 0 (e −2εgs∆τ + e −2εex∆τ α 2k w 0 )/(1 + α 2k w 0 ). It is easily proved that p k+1 > p k for an arbitrary k since α < 1. This means that the success probability at each measurement increases monotonically as the PITE procedure continues and it reaches the saturated value p ∞ ≡ m 2 0 e −2εgs∆τ , independent of the initial relative weight.

For approximate circuit
We next examine the case where the approximate circuit C (1) PITE in Fig. 1(b) is used. The RTE operator for the present system where ε ≡ (ε ex + ε gs )/2, constitutes the circuit, as depicted in Fig. 4(b). The success probability at the kth measurement is calculated as where w (1) k is the relative weight of excited state in this case. The normalized success state is where γ λ ≡ cos(θ 0 − ε λ s 1 ∆τ ) + sin(θ 0 − ε λ s 1 ∆τ ) (λ = gs, ex). By introducing the decay factor α ′ ≡ γ ex /γ gs for the approximate circuit, we obtain w (1) k = α ′2k w 0 for the exponential decay of the relative weight similarly to the case of exact circuit. The number of steps necessary . 3. (a) 2n-qubit gate UME for generating the maximally entangled state from |0 ⊗n ⊗ |0 ⊗n . (b) (2n + 1)-qubit circuit C Gibbs for probabilistic preparation of the Gibbs state of an n-qubit system. This circuit consists of CPITE in Fig. 1(a) for the target system and UME. If the measurement outcome is |0 , the target system has collapsed to the Gibbs state. The success probability gives the partition function.
The success probabilities p k for the exact circuit and p (1) k for the approximate circuit at each step as functions of k for ε gs = 0 and m 0 = 0.8 are plotted in Fig. 5(a). Although those probabilities at k = 0 are smaller for a larger ∆τ as expected, the order is reversed beyond for k > 3 and they increase monotonically toward the saturated value. Figure 5(b) is the plot of the probability P k ≡ k k ′ =0 p k ′ that all the measurements in the exact circuit by the end of kth step are successful. The probability P (1) k for the approximate circuit is also defined similarly. We can see that the differences between the probabilities at each k in terms of the circuits and ∆τ become smaller as the steps continue. Figure 5(c) is the plot of the relative weights w k and w (1) k . For both of the circuits, the weights decay more rapidly for a larger ∆τ , as explained above. The approximate circuit exhibits the smaller weight than the exact circuit does. Apart from this favorable behavior, which might be accidental in the present case, the asymptotic behavior found for this system can be qualitatively true for generic systems. It is because any system whose state is sufficiently close to the ground state can be treated approximately as a two-level system if there is no degeneracy.

B. Wave functions for an asymmetric double-well potential
We consider here a quantum mechanical particle confined to an asymmetric double-well potential in one- PITE PITE circuits for the two-level system, as special cases of those in Fig. 1. The input register in each of them consists of a single qubit. R (1) z ≡ Rz(−s1∆ε∆τ ) is used in this figure. The input state |ψ k to CPITE for the kth step undergoes the gate operations to form the state shown near the dashed line, immediately before the measurement. This state is entangled with the ancilla and contains the success state |ψ k+1 for the next step. It is similarly the case with |ψ dimensional space by using the approach described in Sect. II B.
The potential shape is drawn as the black curve in Fig. 6(a), whose expression is provided in Appendix D. We encoded the wave function of the particle of unit mass at the discretized points by using six qubits. The energy eigenstates of this system, also shown in the figure, are obtained by numerically diagonalizing the Hamiltonian matrix of dimension 2 6 = 64.
Given the double-well structure, it is reasonable to adopt an initial state for the PITE steps having two peaks at the locations of minima. By using the normalized Gaussian wave function |g(x c , σ) centered at x c with width σ, we defined the initial state |ψ (d) init ∝ |g(L/2+d/2, d/3) +|g(L/2−d/2, d/3) /2 with the length L of simulation cell and the distance d between the potential minima. We simulated the PITE procedure using m 0 ≡ 0.9 and ∆τ ≡ 0.1. Some of the success states |ψ k during the steps are plotted in Fig. 6(b), where the initial state converges correctly to the ground state. The probability p k for obtaining the success state at each step is plotted in Fig. 6(c), where it converges to the saturated value close to m 2 0 . The weight of each energy eigenstate in the input state at each step is plotted in the left panel of Fig. 6(d). Since the weight of ground state is rather high already in the initial state, the rapid convergence is achieved.
To see the effects of initial state on the success states during the steps, we also defined the initial state |ψ (s) init = |g(L/2 + d/2, d/3) having only a single peak at the lower minimum and simulated the PITE procedure. The weight of each energy eigenstate is plotted in the right For εgs = 0, m0 = 0.8, and some values of ∆τ , the success probabilities for the exact CPITE and approximate C (1) PITE circuits at each step k are plotted in (a). The probability P k that all the measurements in the exact circuit by the end of kth step are successful is plotted in (b). The similarly defined probability P (1) k for the approximate circuit is also plotted. (c) is the plot of the relative weights of excited states.
panel of Fig. 6(d). On the contrary to the case of |ψ (d) init , the convergence is rather slow. This result and that for a harmonic potential (see Appendix E) tell us that the rapid convergence to the ground state in a PITE calculation requires a good initial guess, which might be also the case for VITE and QITE calculations. When a molecular system is given and we try to decide an initial guess for PITE steps, a typical way for finding a good guess is to perform a cheap classical calculation such as that based on the Hartree-Fock theory [45] or the density functional theory. [57,58] The following two points should, however, be kept in mind. First, such cheap calculations may not give a good initial guess for a molecule where electronic correlation is strong. Second, the state prepa-ration for an adopted initial guess via gate operations become more complicated in general as we treat a larger molecule and/or the spatial resolution for the guess is higher. The second point will force us to find compromise between the resolution and the operation number in the state preparation.

IV. SUMMARY
In summary, we have proposed a generic approach of PITE which constructs the circuit from the RTE gates as black boxes within the first-order of imaginary-time step. We discussed the overall PITE procedure of quantum chemistry calculations and the scaling of computational cost to propose the FQE framework. As demonstrated in the simple examples, the PITE approach suffers from the inherent drawback, i.e., the probabilistic nature. If we overcome it to some extent with the aid of the QAA techniques, we will be able to concentrate on the development of efficient unitary algorithms for real-time dynamics. The examinations on the quality of existing and new QAA techniques adapted to PITE will thus be important in the future in order for this nonvariational approach to be a standard for practical quantum computation. We describe here the change in an input many-qubit state |ψ via the gate operations in the exact PITE circuit C PITE in Fig. 1(a).
Before examining the circuit, we provide an additional explanation of the definition of Θ. Since the inverse trigonometric function arccos for an operator is defined via Taylor series as usual and its value is within the range (0, π), Θ defined in Eq.
From double-peak initial state From single-peak initial state we can set the origin of energy in the ITE operator sufficiently close to the lowest energy eigenvalue, so that the Taylor expansion in H∆τ is justified and √ 1 − M 2 is real.
The composite state of the input register and the initialized ancilla changes via the gate operations as where the RHS of the equality gives the resultant state in Eq. (3).
We describe here the algorithm for probabilistic evolution of a quantum state |ψ(t) whose dynamics is governed by a differential equation of the form d|ψ(t) /dt = L|ψ(t) . The formal solution for a time step ∆t is |ψ(t + ∆t) = e L∆t |ψ(t) . If the generator L of the evolution is anti-Hermitian, the dynamics is unitary and thus easy to implement without any measurement. If the generator is Hermitian, C PITE circuit or its approximate version described above is applicable with an ancillary qubit. Therefore we assume it to be non-Hermitian, for which we demonstrate its probabilistic evolution can be implemented by using two ancillary qubits.
Since L + L † is Hermitian, we can define a nonunitary Hermitian operator where m 0 is a constant similar to the case for PITE. In addition, by using the fact that L − L † is anti-Hermitian, we define a non-Hermitian unitary operator We define the operator Θ and the constants κ and θ 0 similarly to the case of PITE. With them, we introduce two ancillary qubits to construct the circuit C L for the (n + 2)-qubit system, as shown in Fig. 7 which coincides with the correct evolved state e L∆t |ψ within the first order of ∆t. If the measurement outcome is |0 ⊗|0 , the success state e L∆t |ψ with a normalization constant will be obtained. It is easily understood that the success state via the steps will eventually fall onto the subspace corresponding the Jordan block having the lowest real part of eigenvalue among the blocks. It is noted here that there exists a quantum algorithm for solving an inhomogeneous linear equation proposed by Xin et al. [59] They expand the formal solution in time up to a finite order, for which the expansion coefficients are obtained probabilistically by introducing ancillary bits.
Appendix C: Kinetic propagator for one-dimensional dynamics From eqs. (6) and (7), the action of kinetic-evolution operator on the position eigenstate |k n (k = 0, . . . , N − 1) is explicitly calculated as where we have defined the kinetic propagator for ℓ = −N + 1, −N + 2, . . . , N − 1. This satisfies clearly the symmetry J(ℓ; λ) = J(−ℓ; λ) and the probability conservation N −1 ℓ=0 J(ℓ; λ) = 1. The propagator for five qubits is plotted in Fig. 8 as an example. We can understand easily that the periodic boundary condition for kinetic-evolution has entered naturally our formulation and that the evolution operator affects the amplitude of wave function mainly around a point at which the particle is located.
The matrix element of kinetic-energy operator is similarly calculated as  double-peak shape of the wave function has disappeared at the step k = 3 and the wave function at k = 10 already has a large overlap with the ground state. The probability p k for obtaining the success state at each step is plotted in Fig. 9(c), where it converges to the saturated value close to m 2 0 . The weight of each energy eigenstate in the input state at each step is plotted in the left panel of Fig. 9(d), where the weight of ground state increases monotonically while those of the second and third excited states decay rapidly.
We also prepared the superposition of the excited states with odd parity as the initial state |ψ (odd) init = (|φ ex1 + |φ ex3 + |φ ex5 )/ √ 3 and simulated the PITE procedure with ∆τ ≡ 0.1. The weight of each energy eigenstate during the steps are plotted in the right panel of Fig. 9(d), where only the first excited state survives the steps and the other two states decay to zero. This result indicates that the odd parity in the initial state was preserved correctly during the steps and the lowestenergy state |φ ex1 within the odd-parity subspace was obtained since the initial state had an overlap with it.
From asymmetric initial state From odd initial state FIG. 9. (a) Circles represent the exact ground state |φgs and the three lowest excited states |φexµ (µ = 1, 2, 3) at the discrete points encoded by six qubits for the particle in the harmonic potential V (x), shown as the black curve. (b) Red and purple circles represent |φgs and |ψinit , respectively. The other circles represent the success state |ψ k (k = 0, 1, . . . ) obtained immediately after the measurement at the PITE step k. (c) Probability for obtaining the success state at each PITE step. (d) Left panel shows the weight of each energy eigenstate in the input state at each step for the asymmetric initial state |ψinit , while the right panel shows that for the odd initial state |ψ (odd) init with ∆τ = 0.1. tum computer using quantum imaginary time evolution, Nature Physics 16, 205 (2020).