Adaptive Variational Quantum Imaginary Time Evolution Approach for Ground State Preparation

An adaptive variational quantum imaginary time evolution (AVQITE) approach is introduced that yields efficient representations of ground states for interacting Hamiltonians on near-term quantum computers. It is based on McLachlan's variational principle applied to imaginary time evolution of variational wave functions. The variational parameters evolve deterministically according to equations of motions that minimize the difference to the exact imaginary time evolution, which is quantified by the McLachlan distance. Rather than working with a fixed variational ansatz, where the McLachlan distance is constrained by the quality of the ansatz, the AVQITE method iteratively expands the ansatz along the dynamical path to keep the McLachlan distance below a chosen threshold. This ensures the state is able to follow the quantum imaginary time evolution path in the system Hilbert space rather than in a restricted variational manifold set by a predefined fixed ansatz. AVQITE is used to prepare ground states of H$_4$, H$_2$O and BeH$_2$ molecules, where it yields compact variational ans\"atze and ground state energies within chemical accuracy. Polynomial scaling of circuit depth with system size is demonstrated through a set of AVQITE calculations of quantum spin models. Finally, it is shown that quantum Lanczos calculations can also be naturally performed alongside AVQITE without additional quantum resource costs.


Introduction
Quantum computers promise to solve certain types of classically difficult problems more efficiently, with quantum simulation as an important example [1]. In the long term, given access to fault-tolerant quantum computers, adiabatic state preparation followed by quantum phase estimation may become the standard algorithm to determine the ground state energy of a quantum chemistry Hamiltonian [2,3,4]. The required circuit depth, however, is beyond capabilities of near-term noisy intermediate-scale quantum (NISQ) devices, making low-depth hybrid quantum-classical algorithm such as the variational quantum eigensolver (VQE) much more promising to achieve quantum advantage [3,4,5,6]. VQE takes the expectation value of a Hamiltonian, which is measured on a quantum device, as a cost function with a set of variational parameters that are optimized using classical algorithms. Excited-state calculations using VQE have also been proposed by modifying the cost function (e.g. by replacing the energy by the energy variance) or by low energy subspace expansion through linear response [5,7,8,9,10,11]. Meanwhile, quantum imaginary time evolution (QITE) has been developed as an alternative approach to prepare ground states on quantum computers [12]. QITE inherits the advantage of classical imaginary time evolution algorithm, which allows correlations to build faster than would be allowed by the Lieb-Robinson bound that governs real time evolution [13], and always converges to ground state [12]. The key idea of QITE is to represent the imaginary time propagator by unitary operators via least-square fitting [12]. The QITE circuit depth grows exponentially with the correlation domain sizes (being roughly the system's correlation length) and linearly with the number of imaginary time steps. For practical implementations, strategies to reduce circuit complexity, e.g., by utilizing symmetries and effectively combining unitaries [14,15,16,17], have been proposed. Meanwhile, the energy cost function of a fixed variational ansatz can be minimized following the variational QITE (VQITE) approach [18,19], which is a special case of VQE with quantum natural gradient optimization [20,21]. The VQITE method has the advantage of maintaining a fixed circuit depth along the imaginary-time path, but its accuracy is limited by the fidelity of the variational ansatz in representing the ground state. While various strategies to construct variational ansätze have been reported during the development of VQE [22,6,23,24,25,26], their accuracy can be system-dependent and the variational circuits can often be suboptimal [27,28]. One promising strategy is to perform VQE with an adaptively generated ansatz, where operators are drawn from a predefined operator pool and iteratively added to the ansatz during the calculation [29]. It was shown that adaptive VQE can yield highly accurate and compact variational ansätze for specific problems [29,30,31,32,33]. In this work, we develop an adaptive VQITE (AVQITE) method to perform quantum imaginary time evolution for efficient, high-fidelity ground-state preparation for interacting quantum systems. AVQITE bridges the QITE and VQITE approaches by evolving a quantum state in the system's Hilbert space towards the ground state, similar to QITE, yet with a circuit that grows sublinearly and saturates with imaginary time, leading to a final compact variational ansatz. The method generalizes the recently proposed adaptive variational quantum dynamics simulation (AVQDS) method [34] from real to imaginary time. Like AVQDS, the variational ansatz in AVQITE is automatically generated by choosing optimal multi-qubit Pauli rotation gates along the dynamical path to keep a measure of ansatz quality, the McLachlan distance [18,35], within a desired accuracy. We demonstrate the capabilities of AVQITE calculations by preparing the ground states of an H 4 chain and H 2 O and BeH 2 molecules at representative bond lengths with increasing electron correlation effects. We find total energies to chemical accuracy with compact ansätze similar to qubit-ADAPT-VQE results [30], yet without resorting to the complex optimization of parameters in a high-dimensional nonconvex energy landscape. Furthermore, Quantum Lanczos eigenvalue calculations can be carried out together with AVQITE. We additionally demonstrate the favorable polynomial system-size scaling of AVQITE calculations for local spin models. We envision AVQITE, with its compact variational circuits and avoidance of explicit high-dimensional optimization, as a viable way to efficiently prepare ground states of interacting fermion systems (e.g., molecules) on NISQ devices.

AVQITE Algorithm
The AVQITE algorithm generalizes the recently introduced AVQDS approach [34] from real to imaginary time evolution. By time evolving a quantum state in imaginary time, it yields the ground state of a system as the final state. While the derivation of AVQITE resembles that of AVQDS, there are a few key differences that we point out in the following.

Algorithm
The theory of VQITE has been developed in reference [18] and has been used to find ground-state energies of H 2 and LiH molecules on classical simulators [18,19]. Here, we review the VQITE formalism within the density matrix approach, which is insensitive to the global phase of the quantum state. Consider a system with HamiltonianĤ in a pure state |Ψ . The evolution of the density matrixρ ≡ |Ψ Ψ| under the imaginary-time propagator e −τĤ is governed by the Liouville-von Neumann-type equation [18,36] dρ which is a quadratic function of the time derivatives {θ µ ≡ ∂θ µ /∂τ }. The real symmetric N θ × N θ matrix M is specified as which is equivalent to the quantum Fisher information matrix [20,21,37]. The real vector V of dimension N θ is defined as The McLachlan variational principle amounts to minimizing the quadratic cost function L 2 with respect to {θ µ }. This leads to the following linear equations of motion: which has the same form as the parameter update rule arising within the quantum natural gradient approach [20,21]. The second term in equation (3) originates from the global phase of the wavefunction [19,18]. Below we adopt a purely real wave function ansatz in a pseudo-Trotter form (with odd number of Pauli-Y terms) for which the global phase contribution vanishes. Consequently, the density-matrix and wavefunction-based derivations reach exactly the same results. The optimal McLachlan distance L 2 of the variational ansatz Ψ[θ] given by VQITE is formulated in the exactly same form as real time variational quantum dynamics simulations (VQDS) [18,34], with exception of the definition of vector V in Equation (4) due to the different superoperator L[ρ] in the Liouville-von Neumann-type equation.

Flowchart
A typical VQITE calculation, which integrates the equation of motion (6) with a constant time step ∆τ according to the Euler method, is illustrated in the green charts of Figure 1. Given a fixed variational ansatz initialized to a reference state |Ψ 0 that is easily prepared on a quantum computer, the energy expectation value is first measured. The convergence criterion, such as the energy difference between two  consecutive steps, is checked. If the convergence condition is not satisfied, the real symmetric matrix M in Equation (3)  is a product of parametrized native gates of the real device, e.g., single-qubit rotation gates plus two-qubit entangling gates [22]. Alternatively, at a higher algorithmic level,Û[θ] can be expressed as a product of N θ multi-qubit rotation gates in a pseudo-Trotter form: whereÂ µ are Hermitian operators. The unitary coupled cluster ansatz and its variants [6,23,24], the Hamiltonian variational ansatz [25], and the ansatz in the quantum approximate optimization algorithm (QAOA) all belong to this category [26]. In AVQITE, we adopt a variational ansatz in the above form (8), and allow the number ofÂ-operators to be dynamically expanded along the imaginary-time evolution path to maintain high accuracy in representing the evolving quantum state, as shown in the blue module of Figure 1.
The initial steps of an AVQITE calculation are the same as those of a VQITE calculation, with the exception of great flexibility in the choice of initial ansatz |Ψ[θ] . While the simulation accuracy is tied to the initial ansatz in VQITE calculations, an AVQITE calculation monitors the quality of |Ψ[θ] in representing the evolving quantum state and adaptively expands the form of |Ψ[θ] to maintain a fixed accuracy. In fact, AVQITE calculations can simply take any |Ψ 0 as the initial ansatz. For convenience, we use a product state as our initial state in all our calculations below, which is easily prepared on a quantum processor unit (QPU). In the first iteration where no variational parameters are present, the McLachlan distance L 2 in Equation (7) is determined by the energy variance var[Ĥ] in state |Ψ 0 , which is generally larger than the threshold L 2 cut , since |Ψ 0 is not an eigenstate ofĤ. As a result, a predefined operator pool is scanned and an operatorÂ ν is chosen to construct a unitary e −iθ Â ν to be appended to the variational ansatz |Ψ[θ] , which produces minimal McLachlan distance, as illustrated in Figure 1. The additional parameter θ associated with the new operator is always initialized to zero to keep the imaginary time evolution of the variational state continuous. Nevertheless, the McLachlan distance can still be reduced by addition of e −iθ Â ν , because it involves derivatives of the ansatz. The adaptive procedure of selectively appending a new unitary to the ansatz continues until the updated McLachlan distance satisfies L 2 < L 2 cut . The expanded variational parameter vector θ is subsequently updated at the imaginary time step as in VQITE, and new iterations proceed until energy convergence is reached.

Important Technical Details
The accuracy of AVQITE calculations is controlled by the McLachlan distance threshold L 2 cut , while the ability to reach a compact final ansatz |Ψ[θ] is tied to the operator pool. For quantum chemistry calculations, different ways to efficiently construct operator pools and pool completeness conditions have been extensively discussed in the context of adaptive approaches to VQE [29,30,31]. The fermionic operator pool proposed in reference [29] is composed of the single excitation operators and double excitation operators with respect to a Hartree-Fock (HF) reference state. When translated to a qubit representation using encoding methods such as Jordan-Wigner (JW) mapping [38], a single excitation operator can result in a weighted sum of two Pauli strings, and six Pauli strings for a double excitation. Here a Pauli string is defined as a product of single qubit Pauli operators, namely X, Y and Z. Alternatively, a qubit operator pool can also be constructed directly with rudimentary Pauli strings [30,31]. Compared with fermionic operator pools, ADAPT-VQE calculations with qubit operator pools generate variational ansätze with considerably shallower circuits at the price of more variational parameters. As the explicit complex nonconvex optimization problem in the high-dimensional parameter space of VQE is avoided in AVQITE, the qubit operator pool is therefore very appealing and adopted in the following AVQITE calculations. In the AVQITE calculations, we construct a qubit operator pool by choosing all the Pauli strings present in the fermionic single and double excitation operators of the unitary coupled cluster (UCCSD) ansatz [6,23]. The parity mapping is used to transform the fermionic excitation operators to qubit operators [39,40], as fermion to qubit mappings other than JW have not been studied before in the context of operator pool construction [30]. Since the fermionic excitation operators are real, every Pauli string in the operator pool contains an odd number of Pauli Y operators and is therefore antisymmetric. The unitaries in equation (8) are thus real and an initially real wave function (i.e. with real coefficients) remains real when evolving along the imaginary time path. Consequently, the expres- vanishes for any µ. Therefore, the second term of M in Eq. (3) which originates from the global phase of the wavefunction vanishes, and the density-matrix and wavefunction-based approaches lead to exactly the same results. To evolve the quantum state |Ψ[θ] along the imaginary time path by updating parameters as the imaginary time step ∆τ needs to be small enough to maintain high numerical accuracy, while being sufficiently large for fast convergence. In practical AVQITE calculations of molecules (with energy measured in Hartree atomic units) and local spin models (with energy measured in units of the spin-spin in-teraction J), we find ∆τ = 0.1 works well, leading to a ground-state solution with relatively fewer steps.
To stabilize the inversion of the matrix M against potential proximity to singularity in numerical simulations, we adopt the Tikhonov regularization approach by adding a small diagonal element ξ = 10 −6 to M . Alternatively, the matrix inversion problem can be avoided by solving the linear equation (6) using optimization techniques [41].

Implementation Strategies and Measurement Costs on Real Devices
The implementation of AVQITE amounts to measuring the symmetric matrix M in Eq. (3), the vector V in Eq. (4) and the scalars Ĥ θ , Ĥ 2 θ on a quantum device, all of which have been discussed in the context of VQITE in reference [6,42]. Specifically, the energy Ĥ θ can be obtained as a weighted sum of expectation values of the Pauli strings in the Hamiltonian, a procedure often termed "Hamiltonian averaging" [6]: Here h i andÔ i are the coefficient and Pauli string of the ith Hamiltonian term in the qubit representation. The expectation value Ĥ 2 θ can be directly measured similarly by replacingĤ withĤ 2 . The circuit implementation to measure V , which is essentially the energy gradient, has also been discussed in the context of VQE optimization [43]. The proposed indirect measurement circuit introduces an ancillary qubit with a Hadamard-type test. Two controlled-unitary gates, specifically, controlled multi-qubit Pauli gates, need to be implemented. Alternatively, general strategies to replace indirect measurements by direct measurements have also been proposed [44]. It is especially appealing to use the parametershift rule to evaluate the energy gradient V [45,46,47,48], as it only requires some additional measurements of Hamiltonian expectation values at shifted parameters, without the necessity of introducing new circuits. For the matrix element where we have used that the global phase term vanishes with the pseudo-Trotter ansatz (8), the diagonal elements can be simplified to This can be measured with the Hamiltonian averaging method for the Hermitian operatorÂ 2 µ . BecauseÂ µ is a single Pauli string for the qubit operator pools adopted here,Â 2 µ is an identity operator and M µµ = 2. The off-diagonal elements of M can be simplified to M µν = 2 Re Ψ µ [θ]|Â µÛµ+1,ν−1 [θ]Â ν |Ψ ν−1 [θ] for µ < ν, which can be measured using generalized Hadamard test circuit [18,42], or direct measurement circuits [44]. Generally for an N -qubit system with HamiltonianĤ composed of N H Pauli strings, and a parameterized ansatz |Ψ[θ] of N θ parameters with an operator pool of dimension N p , the upper bound for the number of distinct measurement circuits N m for AVQITE calculations is given by N H + N 2 H + N θ (N θ − 1)/2 + N p N θ , where N p N θ comes from the operator selection step. It consists of direct measurements of H,Ĥ 2 , and Hadamard tests to measure matrix M . Here we assume the gradient vector V is evaluated using the parameter shift rule [45,46,47,48]. Assuming a polynomial scaling of N H ∝ N h , N p ∝ N p , and N θ ∝ N q , the leading order for distinct measurement circuits becomes N m ∝ N max(2h,2q,p+q) . For the local spin chain model calculations with local operator pool (see Sec. 4), both N H and N p scale linearly with system size N . Therefore, the number of measurement circuits scales as N m ∝ N max (2,1+q) . Quantum chemistry calculations (see Sec. 3) are usually more challenging, with N H ∝ N 4 and N p ∝ N 4 for the adopted UCCSD pool, which leads to N m ∝ N max(8,2q,q+4) . The N 8 scaling for direct measurement ofĤ 2 and the N q+4 scaling for operator screening can be limiting factors. Remarkably, linear-scaling (p = 1) operator pools have also been proposed in the context of quantum chemistry calculations [30]. (Nevertheless, we leave the study of AVQITE calculations with different pools for future work.) Furthermore, the direct measurement circuits for Ĥ 2 θ can be replaced with the quantum power method for more NISQ-friendly scaling [49], where theĤ 2 is approximated as a linear combination of unitaries. Alternatively, the evaluation of Ĥ 2 θ can be completely avoided by modifying the McLachlan distance criterion. We observe that while it is necessary to have Ĥ 2 θ to evaluate the McLachlan distance L 2 , the estimate of the decrease of L 2 upon appending a new unitary to the ansatz is independent of Ĥ 2 θ . Therefore, the criterion L 2 < L 2 cut can be replaced with a Ĥ 2 θ -free criterion based on the maximal reduction of L 2 due to an added unitary generated by an operator in the operator pool, i.e., max Aν ∈P However, this new criterion implies that screening through all operators in the pool P becomes necessary at every time step. By adopting the above way to circumvent the measurement of Ĥ 2 θ , the main measurement cost for AVQITE calculation at each time step is tied to measuring matrix M and gradient V . The measurement cost analysis in this case closely follows Ref. [50] for generic metric-aware variational quantum algorithms. For simplicity, we consider a uniform distribution of measurements among the N 2 θ + N θ N p elements of M and N θ + N p elements of V , where the N p -dependent parts come from the operator screening step in the adaptive ansatz expansion procedure. To reduce the uncertainty in evaluatingθ = M −1 V due to quantum noise to a precision , i.e., N θ i=1 Var[θ i ] = 2 , the total number M of measurements is upper bounded as: Compared with the measurement cost for a fixed ansatz (see Eqs. (3) and (4) in Ref. [50]), a factor of (1 + Np N θ ) is introduced due to the additional N θ N p measurements for M and N p measurements for V in the operator selection step. The first term in the above equation is associated with measurements of |V | max is the largest absolute value of all elements of V . Because the generator of a multi-qubit Pauli rotation gate is a single Pauli string in the current calculations, the constant factor f F obeys f F ≤ 2 [50]. The constant factor f V depends on the Hamiltonian measurement circuits, and is bounded as 1 ≤ f V ≤ N H . In the above analysis, the indirect Hadamard test circuits are used to measure the gradient V . The measurement cost for estimating the quantum Fisher information matrix M relative to that for V is shown to be asymptotically negligible with increasing number of iterations, due to the diminishing gradient |V | max [50].

Quantum Lanczos Calculation
The focus of the QITE approach lies on the final quantum state, which converges to the ground state for large enough time. Within QITE, all previous quantum states along the path are discarded. In contrast, the quantum Lanczos (QL) method was developed to achieve a more efficient calculations of ground and excited state energies by exploiting information contained in all quantum states along the path [12]. The essence of QL is to diagonalize the Hamiltonian within the Krylov subspace spanned by a subset of imaginary time states {|Ψ n = C n e −n∆τĤ |Ψ 0 } with even time step indices n = 0, 2, . . . N and normalization constants C n . In the classical Lanczos algorithm, the reduced Hamiltonian in the Krylov subspace is brought to a tridiagonal form by sequentially applying the Hamiltonian on orthonormalized Krylov basis vectors [51]. In contrast, QL directly represents the reduced HamiltonianĤ in the normalized basis {|Ψ n } characterized by an overlap matrix S nn = CnC n C 2 (n+n )/2 . As the normalization coefficient C (n+n )/2 is used, only the even-n imaginary-time states (for which (n + n )/2 is an integer) are chosen for the construction of the Krylov subspace. Accordingly, the dense Hamiltonian matrix elements can be evaluated as H nn = S nn E (n+n )/2 , where the expectation value E n ≡ H nn = Ψ n |Ĥ |Ψ n . The normalization coefficient C n can be calculated recursively as C −2

AVQITE calculations of molecules
The ab initio nonrelativistic molecular electron Hamiltonian is given bŷ where the one-electron core part h pq = drφ * p (r)(T + V ion )φ q (r), and the two-electron Coulomb integral is obtained as Here p, q, r, s are composite indices for atom and orbital, and σ is the spin index. T , V ion and V ee are the kinetic energy, ionic potential operator and Coulomb interaction operator, respectively. The standard STO-3G minimal basis set is adopted for the basis orbital functions {φ(r)}. In the following AVQITE calculations, the PySCF quantum chemistry package is first used to generate the molecular Hamiltonian (13) and produce the restricted Hartree-Fock(HF) solution [53]. A basis transformation from atomic orbitals to molecular orbitals is performed to the Hamiltonian (13) for the convenience of preparation of the initial HF state as a tensor product state on the quantum computer. The parity transformation is applied to get the qubit representation of the molecular Hamiltonian, where two qubits are tapered due to the Z 2 symmetry from total electron number and total spin conservation [39,54,55]. All the molecular calculations reported below, including AVQITE and qubit-ADAPT-VQE, are performed in this representation.
In Figure 2 we illustrate a detailed numerical AVQITE calculation of an H 4 chain molecule at a uniform bond length R = 1.5Å, where significant electron correlation effects are present. Starting with the initial Hartree-Fock state |Ψ 0 , the total energy monotonically decreases and converges toward the exact result with increasing imaginary time τ = N ∆τ at ∆τ = 0.1, as shown by the blue circles in Figure 2(a). The associated error, which is defined as the difference between the AVQITE energy and the exact result from a full configuration interaction (FCI) calculation E − E Exact , is shown in the inset of Figure 2(b) on a log scale. The AVQITE calculation reaches chemical accuracy of 1 kcal/mol after τ = 4.7. The quantum state fidelity [1] f ≡ | Ψ[θ]| Ψ Exact | 2 , i.e., the squared overlap between the ansatz state and the exact ground state, goes beyond 99.9%. The number of variational parameters N θ , or equivalently the number of multi-qubit Pauli rotation gates, increases at several imaginary time steps and flattens at 24. For comparison, the qubit-ADAPT-VQE calculation of H 4 using the same operator pool generates a final ansatz of 30 variational parameters upon reaching chemical accuracy. Specific information about the McLachlan distance L 2 is plotted in Figure 2(c), where L 2 is reduced below the threshold L 2 cut = 5 × 10 −4 by adaptively expanding the variational ansatz whenever the initial value L 2 > L 2 cut at any time step. Together with L 2 , the energy variance var θ [Ĥ] along the imaginary path is shown to be almost linear on a semi-log scale, implying an exponential convergence. In fact, the AVQITE total energies E(τ ) can also be fitted with an exponential function, E(τ ) = E ∞ + e −aτ , with E ∞ = −1.9957 Ha in the infinite time step limit, which is within chemical accuracy with an error of 0.2 kcal/mol. Alternatively, one can fit a function E(var θ [Ĥ]) = E 0 + bvar θ [Ĥ] due to the approximately linear relation between E and var θ [Ĥ] in the numerical results, where E 0 = −1.9959 Ha in the zero-variance limit is also very accurate with an error 0.1 kcal/mol. The quantum Lanczos results, which are conveniently calculated along with AVQITE at no additional quantum resource cost, are plotted as orange squares in Figure 2(a) for the total energy, with the error shown in the inset of panel (b). The QL calculation converges relatively faster than AVQITE, since the reduced Lanczos Hamiltonian encodes information beyond the latest quantum state. Because the QL method is generally susceptible to numerical inaccuracies and the way to evaluate the Hamiltonian and overlap matrix in the Krylov subspace used in the QL method accumulates sizable errors, excited states cannot be directly accessed in the current QL calculations. To demonstrate the general applicability of the AVQITE approach in generating compact ground state ansätze, we perform AVQITE calculations for H 4 chains at bond length R HH = 0.8Å, 1.5Å, 2.4Å, H 2 O at R OH = 0.8Å, 1.5Å, 2.4Å, and BeH 2 at R BeH = 0.8Å, 1.4Å, and 3.6Å, as shown in Figure 3. This benchmark set covers a variety of directional covalent bonding and electron correlation effects, with atomic states of spin singlet, doublet and triplet in the dissociation limit. In Figure 3(a) we show the AVQITE energies of H 4 at three bond lengths colored in red, green and blue, which agree with the FCI results within chemical accuracy. The detailed error convergence, E − E Exact , as a function of imaginary time step τ = N ∆τ is shown in panel (d). With increasing bond length or correlation energy, the critical imaginary time τ c , which is the time where chemical accuracy is reached, generally increases. Here the correlation energy is defined as the difference between the HF and FCI energies, which corresponds to the initial AVQITE energy error as the AVQITE ansatz starts with the HF state. At R HH = 2.4Å proximate to the atomic limit, the critical time step increases significantly to τ c ≈ 110. However, as observed earlier in numerical step-merged QITE calculations [14], one can choose a bigger step size ∆τ for molecules at larger bond length to reduce the number of imaginary time steps N . In this case, ∆τ is increased from 0.1 to 0.5 to speed up the convergence, although sizable energy fluctuations are present in the initial time steps. The number of controlled-NOT gates (CNOTs) generally increases in the initial steps and levels off for τ above 3, as shown in panel (g). The maximal number of CNOTs reaches 250 for H 4 at R HH = 2.4Å. Here the number of CNOTs is estimated by the rule that each multi-qubit rotation gate e −iθσ with Pauli stringσ of length p needs 2(p − 1) CNOTs (assuming all-to-all connectivity) [1]. Moving on to the AVQITE calculations of H 2 O and BeH 2 at representative bond lengths, the error convergence behavior generally remains similar to the H 4 calculations, with final energies within chemical accuracy. Likewise, the number of CNOTs in the AVQITE state preparation circuit generally increases initially and levels off for τ > 3. The total number of CNOTs remains close to or under 400. Remarkably, the positive correlation between the number of CNOTs in the AVQITE ansatz and correlation energy, which seems to exist in the calculations of H 4 , does not apply to H 2 O and BeH 2 . The AVQITE state preparation circuits in the test-set are generally over one order of magnitude shallower than the original UCCSD ansatz. In panels (g-i), we show that the number of variational parameters N θ of the AVQITE ansatz is generally quite close to that of the qubit-ADAPT-VQE approach. The H 4 chain at R HH = 2.4Å marked by a blue cross symbol in Fig. 3(j) represents an interesting example, where the qubit-ADAPT-VQE becomes trapped in a local minimum with energy 15 mHa higher than the groundstate energy. The qubit-ADAPT-VQE method is implemented following references [27,30]. The Broyden-Fletcher-Goldfarb-Shannon (BFGS) algorithm is used to optimize the variational ansatz in a new generation with an additional parameterized unitary at the optimal solution of the previous generation. This implies that the energy cost function is optimized along a continuous parameter path, which could lead to a local minimum with vanishing gradients in a generally nonconvex high-dimensional energy landscape, as exemplified here. This is related with the general barren plateau problem where the optimization of random circuits is trapped in a flat cost function region in the parameter space [56]. Various techniques to address the barren plateau problem have been proposed, including the block-identity initialization strategy [57], and entanglement devised mitigation techniques [58]. In fact, a modified qubit-ADAPT-VQE calculation, which successively adds unitaries in accordance with the AVQITE calculation and optimizes the ansatz initialized at the AVQITE parameter values already reaches the accuracy of 0.1mHa with the first 31 unitaries out the total 33 unitaries in the AVQITE ansatz. This implies a proper combination of qubit-ADAPT-VQE and AVQITE can be mutually beneficial.  4 System-size scaling of AVQITE circuit complexity Quantum resources for AVQITE calculations have been analysed in Sec. 2.2.3. A measure of NISQ circuit complexity is the number of CNOTs N cx , which is related to the number of parameters N θ associated with multi-qubit Pauli rotation gates. It is assumed that N θ grows polynomially as N q with system size N , where the exact order q is tied to ground state complexity of the system. Because of the limitations of classical simulations, it is difficult to show directly the system-size scaling of N cx in quantum chemistry calculations. Instead, we apply AVQITE to prepare ground states of local spin chain models, where the scaling of N cx is more easily accessible. We find that N cx AVQITE scales quadratically or linearly, depending on how close the system is to a quantum critical point. This scaling is similar to VQE using the well-established Hamiltonian variational ansatz (HVA) for local spin models [25,59,60]. We thus conjecture that a favorable polynomial order of system-size scaling for the adaptive preparation of ground states carries over to complex quantum chemistry problems, which should be addressed in future work.
To establish the scaling of N cx with system size, we consider the mixed-field Ising model (MFIM): with Z 0 ≡ Z N for periodic boundary conditions. Energy is measured in units where J = 1. When setting the longitudinal field h z = 0, the MFIM reduces to the integrable transverse-field Ising model (TFIM). The competition between ferromagnetic and paramagnetic phases is controlled by the transverse field h x of the model. When h x = 1, the TFIM is at a quantum critical point. We apply AVQITE to prepare the ground state of the TFIM at the quantum critical point (h x , h z ) = (1, 0), and of the MFIM at (h x , h z ) = (1, 0.5), with 4 ≤ N ≤ 10. The following complete operator pool is adopted for the calculation [61,11]: which is composed of only local one-qubit and two-qubit Pauli strings with a single Y operator. Figure 4 shows the system-size dependence of N cx of the AVQITE circuit which prepares the ground state with fidelity f > 99.9%. N cx grows as N 2 for preparation of the quantum critical state of the TFIM. The finite longitudinal field h z = 0.5 in MFIM breaks the integrability of the TFIM. However, it also drives the system away from the quantum critical point. Consequently, the ground state of the MFIM becomes less entangled, and can be prepared by an AVQITE circuit with N cx = 2N . The above circuit analysis for N cx assumes that the qubit connectivity is cyclic, like the spin-site configuration of the model, such that each two-qubit rotation gate requires 2 CNOT gates. The favorable scaling of N cx in preparing the ground state of local spin models can also be obtained using HVA, as shown in Fig. 4. Perfect fidelity is obtained with an N/2-layer ansatz for a TFIM with even system size N [59,60]: The reference state is chosen to be the ground state of H 1 : |Ψ 0 = ⊗ N i=1 |+ , with |+ = (|↑ + |↓ )/ √ 2. The CNOT gates in the N/2-layer HVA circuit can be simply counted as N cx = 2N × N/2 = N 2 . For the MFIM, we find that a constant two-layer HVA ansatz of the following form can prepare the ground state with similar fidelity f > 99.9%: Here we defineĤ 3 ≡ − N i=1Ẑ i , and choose the ground state ofĤ 3 , |Ψ 0 = ⊗ N i=1 |↑ , as the reference state. The number of CNOT gates can be similarly estimated as N cx = 2N × 2 = 4N . For both local spin models considered, the ground state preparation circuits for AVQITE and HVA have similar complexity as measured by N cx . Nevertheless, they approach the ground state along different paths. AVQITE follows the QITE path starting with the reference state, and the imaginary time evolution approach is guaranteed to converge to the ground state [12]. Compared with the VQE-HVA approach, the operator pool screening procedure to adaptively expand the ansatz introduces substantial overhead for AVQITE calculations. In contrast, proper parameter initialization is important for VQE-HVA calculations. In fact, VQE-HVA calculations must not start from state |Ψ 0 , as otherwise the gradient vanishes. This can be easily achieved by initializing the parameters to nonzero values. For random parameter initialization, a finite fraction of VQE calculations on the TFIM with the N/2-layer HVA ansatz fail to converge [60].

Conclusion
The adaptive variational quantum imaginary time evolution (AVQITE) approach is developed by generalizing the recently proposed adaptive quantum dynamics simulation (AVQDS) method [34] from real time to imaginary time. We present benchmark quantum chemistry calculations on H 4 , H 2 O and BeH 2 with different chemical bonding character and atomic spin multiplicity upon dissociation limit. They demonstrate the general applicability of AVQITE to finding accurate and compact variational ansätze for interacting many-electron models. The key advantage of AVQITE over adaptive VQE approaches is that it bypasses the complicated nonconvex optimization problem in high-dimensional parameter space by performing energy cost function minimization along a single imaginary time axis. AVQITE can also be applied to reduced "low energy" models in quantum chemistry such as those that arise within the complete active space self-consistent field method (CASSCF) [62]. More generally, it can be used as an impurity solver for quantum embedding approaches like the rotationally invariant Gutzwiller embedding and density-matrix embedding methods [63,64,65,66,67,68,69], extending its applicability to larger molecules and solid state materials. With the favorable polynomial scaling of the AVQITE calculations which automatically generates compact variational ansätze of high accuracy, we envision feasible followup applications of this method on NISQ devices to calculate ground-state properties of spin models and molecules. Excited states are accessible if AVQITE is used with cost functions related to the energy variance, such as in the folded spectrum method [5,6,11,70,71].