Boosting quantum amplitude exponentially in variational quantum algorithms

Quantum computers promise exponential speedup over classical counterparts in solving certain tasks [1]. When fault-tolerant general-purpose quantum computers become available, adiabatic state preparation and quantum phase estimation may become the standard quantum routines for determining the ground-state energy of sophisticated physical Hamiltonians [2–5]. However, such schemes are very costly in terms of required overhead and hence are not suitable for the current era of noisy intermediate-scale quantum (NISQ) hardware [6–10]. This limitation of quantum computers today shifts central attention towards low-depth hybrid quantum–classical algorithms, known as NISQ algorithms [11–14]. The variational quantum eigensolver (VQE) [15, 16] serves as a prototypical example, as an algorithm that computes the expectation value of a Hamiltonian, which is measured on a quantum machine, resulting in a cost function with a set of variational parameters, which are optimized using classical computers. The process is repeated until the cost function reaches its local minimum.

On the other hand, the variational quantum simulator [17] has been proposed for hybrid quantum–classical simulations of quantum dynamics based on the McLachlan's variational principle [18, 19], including quantum imaginary time evolution (QITE) to prepare ground states [18–20]. Here, we introduce the 'quantum iterative power algorithm' inspired by the variational quantum simulator to provide an accelerated method to the general problem of global optimization with near term quantum computers.

Global optimization is central to many important problems in science and engineering, from back-propagation in machine learning [21] and molecular geometry optimization/protein structure prediction [22, 23] to route planning and control of drone/unmanned aerial vehicles [24]. However, the brute force approach of considering each possible element of a search space often becomes computationally intractable. For example, identification of the optimal configuration of a protein faces Levinthal's paradox [25]—that the native configuration must be identified out of about 10³⁰⁰ possibilities. This has inspired a broad array of both classical [26] and quantum computing [27] optimizers. Recently, we have shown that tensor trains [28, 29] (also known as matrix product states) provide a way to vastly reduce the computational cost of exploring low-rank optimization cost functions, and have employed the approach to introduce an optimization algorithm that deterministically explores the full search space in data-compressed form, the tensor-train 'iterative power algorithm (IPA)' [30].

We recognize the strategy of tensor-train IPA can be implemented on quantum computers to enable global optimization of an even broader class of optimization problems. In tensor-train IPA, the optimization cost function of interest is taken to be a potential energy surface. A density is initialized in the potential energy surface, and an oracle is iteratively applied in a sifting approach akin to imaginary time propagation (with infinite mass) to localize the density as a delta function at the global minimum position. The expectation value of position then gives the location of the global minimum. Tensor-train IPA represents the density and potential energy surface as tensor trains to avoid calculation of the cost function everywhere in search space, which is efficient for representation of problems amenable to low-rank representations, such as prime factorization or molecular geometry optimization [30]. However, the tensor-train strategy faces the roadblock that highly-correlated systems cannot be efficiently represented in low-rank tensor-train format. In contrast, quantum computers excel in the simulation of highly-correlated systems, as the coupling or entanglement between qubits is limited only by the choice of ansatz [31].

The quantum iterative power algorithm (QIPA) takes advantage of the high degree of entanglement possible on quantum computers with a hybrid variational scheme. In standard variational approaches such as the VQE [15, 16], classical optimizers are used to determine the parameters of a quantum circuit, which are used to prepare trial wavefunctions measured to obtain expectation values. Analogously, the variational quantum simulator [17] evolves the parameters that define the time-evolved wavefunction by using a classical computer that integrates the Euler–Lagrange equation obtained from the Schrödinger equation with the McLachlan's variational principle. Parameters required by the Euler–Lagrange equation are obtained with a quantum circuit with a small number of quantum operations. QIPA generalizes the variational quantum approach to evolve an initial wavefunction such that the corresponding probability density (modulus squared of the wavefunction) becomes localized at the global minimum of a given cost function (see figure 1). As in IPA, the propagator of QIPA is not limited to the imaginary time quantum propagator enabling the use of other propagators that are maximal at the minimum of the cost function.

Zoom In Zoom Out Reset image size

Figure 1 and algorithm 1 show the overall workflow of the QIPA. First, we select parameters $\theta_1,\theta_2,\dots,\theta_{\mathcal{N}_\theta}$ corresponding to the initial state. Second, we use a quantum coprocessor to perform the Hadamard tests, resulting in ${A}_{k,m}$ and C_k of $\sum_m {A}_{k,m}\dot{\theta}_m = C_k$. There are a number of ways to get an approximate solution to the linear equation. Here, we used the conjugate gradient method without needing to invert the matrix A to get a solution $\dot{\boldsymbol{\theta}}$, using a classical computer. Next, we update the parameters θ by the Euler method. Having updated the parameters, the process is iterated until convergence to obtain parameters θ corresponding to a distribution function localized at the global minimum.

Algorithm 1. Variational QIPA.
Require: Hamiltonian $\hat{H}$ and initial state $\vert{\psi(0)}\rangle = \vert{0}\rangle^{\otimes\mathcal{N}}$
1: Start with an ansatz $\vert{\phi(\boldsymbol{\theta}_0)}\rangle = U(\boldsymbol{\theta}_0)\vert{\bar{0}}\rangle$ at time τ = 0;
2: Evaluate Hadamard tests to form the matrix $\textbf{A}(\tau)$ and the vector $\textbf{C}(\tau)$ (quantum computer subroutine);
3: Compute an approximate solution ξτ of $\textbf{A}(\tau)\dot{\boldsymbol{\theta}}(\tau) = \textbf{C}(\tau)$ via the conjugate gradient method(classical computer subroutine);
4: Update the parameter as $\boldsymbol{\theta}(\tau+\delta\tau)\leftarrow\boldsymbol{\theta}(\tau)+\xi_\tau\delta\tau$ and set $\tau\leftarrow\tau+\delta\tau$;
5: Repeat steps 2–4 until $\tau = \tau_{\textrm{total}}$ or the convergence criteria is met;

We are interested in a particular case of global optimization involving the search of the ground state of a Hamiltonian $\hat{H}$: a problem that is typically solved by the imaginary time propagation. QIPA can solve the same problem analogously by using $\hat{H}$ in the normalized oracle function $f(\hat{H};\tau)$ that acts on the initial wavefunction $\vert{\psi(0)}\rangle$, as follows (onwards setting $\hbar = 1$):

$$\begin{align} \vert{\psi\left(\tau\right)}\rangle & = f\left(\hat{H};\tau\right)\vert{\psi\left(0\right)}\rangle\nonumber\\ & = \frac{{U}\left(\tau\right)\vert{\psi\left(0\right)}\rangle}{\sqrt{\langle {U}\left(\tau\right)\psi\left(0\right)\vert {U}\left(\tau\right)\psi\left(0\right)\rangle}}, \end{align} \tag{ 1 }$$

where ${U}(\tau)$ is an arbitrary oracle function with maximum at the global minimum position of the potential energy surface M, or here the Hamiltonian $\hat{H}$. In the following, we show that oracles defined by concatenated exponential functions,

$$\begin{align} {{U(\tau)}} = U_n(\tau) = \beta_n(-\hat{H}\tau)) = e^{b_n\beta_{n-1}(-\hat{H}\tau)}, \end{align} \tag{ 2 }$$

with $n\unicode{x2A7E}1$ the number of concatenated exponentials, $\beta_0(y) = y,\,\beta_1(y) = e^{b_1 y},\beta_2(y) = e^{b_2e^{b_1y}},\ldots$ and real constants $b_1,\dots,b_n\neq 0$, provide effective QIPAs based on a generalization of the McLachlan's variational principle (SI appendix A.2.). For example, the oracle defined by the double-exponential $U_2(\tau) = e^{e^{-\tau\hat{H}}}$ is obtained by setting n = 2 and $b_2 = b_1 = 1$.

We remark that the choice of U₁ corresponds to the standard QITE, which is widely used in quantum Monte Carlo algorithms. References [18, 19] show that one can perform imaginary time evolution [20] with unitary gates defined by equation (1) with n = 1 that evolve the initial state according to the Wick-rotated Schrödinger equation: ${\partial\vert{\psi(\tau)}\rangle}/{\partial\tau} = -\left(\hat{H}-E_1(\tau)\right)\vert{\psi(\tau)}\rangle,$ where $E_1(\tau) = \langle{\psi(\tau)}\vert\hat{H}\vert{\psi(\tau)}\rangle$. Here, we introduce a family of near-term quantum algorithms defined by β_n with $n\unicode{x2A7E} 1$ that evolve the initial state according to the generalized Wick-like-rotated Schrödinger equation (SI appendix A.3.):

$$\begin{align} \frac{\partial}{\partial\tau}\vert{\psi\left(\tau\right)}\rangle & = -\prod_{k = 1}^{n}b_k\left(\hat{H}\exp\left(\hat{S}_{n-1}\right) -\textrm{Re}\left\langle\hat{H}\exp\left(\hat{S}_{n-1}\right)\psi\left(\tau\right)\mid\psi\left(\tau\right)\right\rangle\right)\vert{\psi\left(\tau\right)}\rangle,\nonumber \end{align} \tag{ 3 }$$

where $\hat{S}_{n-1} = \sum_{k = 1}^{n-1}b_k\beta_{k-1}(-\hat{H}\tau)$.

With the choice n = 2 and $b_2 = b_1 = 1$ we arrive at a double-exponential function and the following Wick-like-rotated Schrödinger equation:

$$\begin{align} \frac{\partial}{\partial\tau}\vert{\psi\left(\tau\right)}\rangle = -\left(\hat{H}e^{-\hat{H}\tau}-E_2\left(\tau\right)\right)\vert{\psi\left(\tau\right)}\rangle, \end{align} \tag{ 4 }$$

with $E_2(\tau) = \langle\psi(\tau)\mid\hat{H}e^{-\hat{H}\tau}\mid\psi(\tau)\rangle$. According to the McLachlan's variational principle, when we constrain the equation of motion as such

$$\begin{align} \delta\left\|\left(\partial/\partial\tau+\left[\hat{H}e^{-\hat{H}\tau}-E_2\left(\tau\right)\right]\right)\vert{\psi\left(\tau\right)}\rangle\right\|^2 = 0, \end{align} \tag{ 5 }$$

the result is equivalent in finding a solution of the linear equation: $\sum_mA_{k,m}\dot{\theta}_m = C_k,$ where the entries of the symmetric and positive semi-definite matrix A and the right-hand side C can be computed on a quantum computer by deploying the Hadamard tests. The parameters θ are updated with $\dot{\boldsymbol{\theta}}$ for a short time step $\delta\tau\gt0$ according to the Euler method as $\boldsymbol{\theta}(\tau+\delta\tau)\approx\boldsymbol{\theta}(\tau)+\dot{\boldsymbol{\theta}}(\tau)\delta\tau$. The underlying assumption is that we can approximate $\vert{\psi(\tau)}\rangle$ by $\vert{\phi(\boldsymbol{\theta}(\tau))}\rangle = U(\theta_1(\tau))U(\theta_2(\tau))\cdots U(\theta_{\mathcal{N}_\theta}(\tau))\vert{\bar{0}}\rangle$, where $\vert{\bar{0}}\rangle = \vert{0}\rangle^{\otimes\mathcal{N}}$ and $U(\theta_1(\tau)),\dots,U(\theta_{\mathcal{N}_\theta}(\tau))$ are parameterized quantum circuits (PQCs), with $\boldsymbol{\theta} = (\theta_1,\dots,\theta_{\mathcal{N}_\theta})$ the corresponding real-valued parameter vector.

2.1. Exponential amplitude amplification

2.2. Resource estimate and error analysis

2.3. Numerical experiments

Figure 2 shows that one can efficiently search for the ground-state energy of hydrogen molecule across various bond stretching distances with QIPA with fewer time steps than QITE. Results are shown in figure 2(a) comparing QIPA and QITE. For consistency and fair comparison, we use the same time step for each bond distance for both QIPA and QITE runs. The error difference between the exact energy obtained from full configuration interaction (CI) calculations and the QIPA and QITE results can be seen in figure 2(b). QIPA features less error for all bond distances considered even though the accuracy of both QITE and QIPA deteriorates as the system goes into the highly correlated regime. Alternative ansatz circuits could be considered for increased accuracy.

Zoom In Zoom Out Reset image size

As shown in figure 2(c), both QITE and QIPA successfully minimize the energy of a given transmon circuit. Moreover, QIPA requires significantly fewer iterations than QITE for global optimization of the energy. These results constitute the successful implementation of two forms of imaginary-time-like evolution to perform quantum computer-aided design (QCAD), i.e. to think in a 'reverse-engineering way' to find out optimal quantum system parameters, in this case being classical circuit parameters of the superconducting circuit, given one desired energy spectrum of the quantum system. Previous methods [33–35] used VQE approach whereas ours is in the form of variational quantum simulator.

According to figures 3(a) and (b), we find QIPA identifies the prime factors of biprimes with fewer iterations than QITE for all integers factored. The number of required iterations varies depending on the biprime factored and the ansatz, with speedups shown here of up to 50%. In addition, QIPA successfully identifies both factors of biprimes. As expected, QIPA succeeds for stable integration with a small, converged time step, for which the twin solutions are readily identified as the two components of the final wavefunction with the largest and equal amplitude. Moreover, QIPA succeeds for unstable integration with a larger time step, for which the twin solutions are identified as the two largest components of the final wavefunction with unequal amplitude, as depicted in figure 3(c). To our knowledge this ability to identify multiple prime factors is unique among published quantum computing prime factorization results.

Zoom In Zoom Out Reset image size

In summary, we have presented a family of generalized imaginary-time-like near-term quantum algorithms which we coin the 'quantum iterative power algorithm,' inspired by its classical counterpart. (Plural 'algorithms' is used since depending on the choice of oracle function, the performance and behavior will differ. However, they all fall under the same family.) We have analyzed its convergence rate. One caveat is that since the proposed algorithm relies heavily on the ansatz circuit used, its convergence rate is difficult to discern in the generic case. We have also determined QIPA's estimated resource count as well as analytical error analysis, and demonstrated it can outperform the QITE while it reduces the number of required iterations, at the cost of a moderate increase in the number of gates. We note that even when the initial quantum state has an exponentially small overlap with the final target state, QIPA needs only a polynomial number of steps to reach convergence (SI appendix D). This is particularly important when starting with an initial state defined by a uniform superposition, or a low-rank reference state for a highly correlated system [32]. Furthermore, we have used the three numerical case studies—quantum chemical molecular simulation of the hydrogen molecule for various bond dissociation distances, quantum computer-aided design of a superconducting transmon, and finding optimal solutions for double prime factorization—to highlight how QIPA outperforms QITE. The simulation accuracy could be further improved by integrating with existing quantum error mitigation methods [17, 36–41].

We would like to point out that an additional consideration for such an algorithm besides the choice of the ansatz circuit is setting the right parameter for the time step δτ at each evolution step. A number of proposals [42, 43] suggest the use of an adaptive time step to overcome such an issue. However, there exists an opportunity to develop a systematic way to adjust the time step δτ according to gradient descent, rather than with a heuristic argument on how to adjust δτ. A major drawback with quantum imaginary-time-like evolution algorithm is that it involves constructing the matrix elements of A with NISQ hardware. Since we are working with noisy quantum hardware, any large fluctuation in the matrix elements would result in suboptimal angles θ . Future work will explore conditioning the matrix based on errors associated with performing the Hadamard tests.

T H K, C S, and A A-G acknowledge funding from Dr Anders G Frøseth. A A-G also acknowledges support from the Canada 150 Research Chairs Program, the Canada Industrial Research Chair Program, and from Google, Inc. in the form of a Google Focused Award. M B S acknowledges financial support from the Yale Quantum Institute Postdoctoral Fellowship. V S B acknowledges support from the Center for Quantum Dynamics on Modular Quantum Devices (NSF Grant No. CHE-2124511) and computational resources from NERSC and the Yale Center for Research Computing High Performance Computing. We acknowledge fruitful discussions with Jakob Kottmann, Abhinav Anand, Sabre Kais, Eitan Geva, Raja Selvarajan, Suguru Endo, Xiao Yuan, Ellen Mulvihill, and Zi-jian Zhang. Numerical simulations of quantum circuit evaluations for the $\textrm{H}_2$ and transmon quantum systems were performed on the Niagara supercomputer at the SciNet HPC Consortium [44, 45]. SciNet is funded by the Canada Foundation for Innovation, the Government of Ontario, Ontario Research Fund-Research Excellence, and the University of Toronto. The Tequila package [46] has been used to implement all the quantum circuit operations, using Qulacs [47] as a quantum backend.

The data that support the findings of this study are openly available at the following URL/DOI: https://github.com/aspuru-guzik-group/QIPA/.

This section outlines a quantum version of a generalized IPA method [30], the so-called quantum iterative power algorithm (QIPA) inspired by the imaginary time propagation (appendix A.1) and the variational quantum simulator [17]. We introduce a series of oracles defined by concatenated exponential functions which evolve the initial state according to a Wick-like-rotated Schrödinger equation obtained from a generalized McLachlan's variational principle. Then, we focus on the particular case of a double-exponential function that provides an efficient implementation of QIPA. For completeness, we include a section introducing the QITE method obtained from the McLachlan's variational principle.

A.1. Quantum imaginary time evolution

A.1.1. McLachlan variational principle approach

Let us consider a many-body Hermitian Hamiltonian $\hat{H}$. ¹² Given an initial state $\vert{\psi(0)}\rangle$, non-unitary quantum imaginary time evolution (QITE) is defined by

$$\begin{align} \vert{\psi\left(\tau\right)}\rangle = \frac{e^{-\hat{H}\tau}\vert{\psi\left(0\right)}\rangle}{\sqrt{\langle{\psi\left(0\right)}\vert e^{-2\hat{H}\tau}\vert{\psi\left(0\right)}\rangle}}. \end{align} \tag{ A.1 }$$

Note that the denominator is a normalization factor. Since a key to our proposal is the realization that there is nothing except the requirement that there be a continuous, integrable, strictly positive oracle $\hat{U}$ that is maximized at the global minima of $\hat{H}$ to prevent us from assuming a particular form of the oracle; we present QITE more generally in terms of non-unitary time evolution according to an oracle, as follows:

$$\begin{align} \vert{\psi\left(\tau\right)}\rangle = \frac{\beta\left(-\hat{H}\tau\right)\vert{\psi\left(0\right)}\rangle}{\sqrt{\langle{\psi\left(0\right)}\vert\beta\left(-2\hat{H}\tau\right)\vert{\psi\left(0\right)}\rangle}}, \end{align} \tag{ A.2 }$$

where β could be any strictly increasing positive function. When $\beta(y) = e^y$, we recover imaginary time evolution, which corresponds directly to the Wick rotation $(\tau = -it)$. With that choice of $\beta(y)$, we obtain that the above quantum state satisfies the Wick-like rotated Schrödinger equation:

$$\begin{align} \frac{\partial\vert{\psi\left(\tau\right)}\rangle}{\partial\tau} = -\left(\hat{H}-E_1\left(\tau\right)\right)\vert{\psi\left(\tau\right)}\rangle, \end{align} \tag{ A.3 }$$

where $E_1(\tau) = \langle{\psi(\tau)}\vert\hat{H}\vert{\psi(\tau)}\rangle$. Even though $\vert{\psi(\tau)}\rangle$ is a valid wavefunction that can be represented in a quantum computer, the non-unitary time evolution cannot be straightforwardly mapped to a quantum circuit based solely on unitary gates. Here, the McLachlan's variational principle [18, 19, 48, 49] comes to the rescue and demands that

$$\begin{align} \delta\left\|\left(\partial/\partial\tau+\left[\hat{H}-E_1\left(\tau\right)\right]\right)\vert{\psi\left(\tau\right)}\rangle\right\| = 0, \end{align} \tag{ A.4 }$$

with $\|\cdot\|$ representing the L² norm and δ the functional derivative. In the following, we intend to simulate the action of non-unitary dynamics, (A.3), on a quantum computer via McLachlan's variational principle.

In variational quantum simulations, instead of directly encoding the quantum state $\vert{\psi(\tau)}\rangle$ at time τ, we approximate it with a PQC $\vert{\psi(\tau)}\rangle\approx\vert{\phi(\boldsymbol{\theta}(\tau))}\rangle$ with a real-valued parameter vector $\boldsymbol{\theta}(\tau) = (\theta_1 (\tau),\theta_2 (\tau),\dots,\theta_{\mathcal{N}_\theta}(\tau))$. We assume that physically relevant quantum states span a restricted region of the full Hilbert space [50] for a given time interval, such that the trial state parameterized by θ is sufficient to prepare a desired quantum state by applying a sequence of parameterized unitary gates $U(\boldsymbol{\theta}) = U_{\mathcal{N}_\theta}(\theta_{\mathcal{N}_\theta})\cdots U_k(\theta_k)\cdots U_1(\theta_1)$ to the initial state $\vert{\bar{0}}\rangle = \vert{0\cdots 0}\rangle$. Thus, we have $\vert{\phi(\boldsymbol{\theta})}\rangle = U(\boldsymbol{\theta})\vert{\bar{0}}\rangle$, where $U(\boldsymbol{\theta})$ is referred to as the ansatz, and $U_k (\theta_k)$ is the kth unitary gate controlled by classical parameter θ_k . Here, we are only concerned with single- or two-qubit gates, which is sufficient for universal quantum computing.

According to McLachlan's variational principle, we require $\partial\|(\partial/\partial\tau+(\hat{H}-E_\tau))\vert{\phi(\boldsymbol{\theta}(\tau))}\rangle\|/\partial\dot{\theta}_k = 0$. We have

$$\begin{align} & \delta\|\left(\partial/\partial\tau+\left[H-E_1\left(\tau\right)\right]\right)\vert{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\rangle\|^2 \nonumber\\ & \quad = \sum_{m,n}\frac{\partial\langle{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\vert}{\partial\theta_m}\frac{\partial\vert{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\rangle}{\partial\theta_n}\dot{\theta}_m\dot{\theta}_n+\langle{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\vert^2\left(\hat{H}-E_1\left(\tau\right)\right)^2\vert{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\rangle\nonumber\\ &\qquad+\sum_m\frac{\partial\langle{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\vert}{\partial\theta_m}\left(\hat{H}-E_1\left(\tau\right)\right)\vert{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\rangle\dot{\theta}_m\dots\nonumber\\ &\qquad+\sum_m\langle{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\vert\left(\hat{H}-E_1\left(\tau\right)\right)\frac{\partial\vert{\phi(\boldsymbol{\theta}(\tau))}\rangle}{\partial\theta_m}\dot{\theta}_m. \end{align} \tag{ A.5 }$$

By differentiating with respect to $\dot{\theta}_k$, we obtain

$$\begin{align} &\partial\|(\partial/\partial\tau+\left[\hat{H}-E_1(\tau)\right])\vert{\phi(\boldsymbol{\theta}(\tau))}\rangle\|^2/\partial\dot{\theta}_k \nonumber\\ & \quad = \sum_m\left(\frac{\partial\langle{\phi(\boldsymbol{\theta}(\tau))}\vert}{\partial\theta_k}\frac{\partial\vert{\phi(\boldsymbol{\theta}(\tau))}\rangle}{\partial\theta_m}+\frac{\partial\langle{\phi(\boldsymbol{\theta}(\tau))}\vert}{\partial\theta_m}\frac{\partial\vert{\phi(\boldsymbol{\theta}(\tau))}\rangle}{\partial\theta_k}\right)\dot{\theta}_m\nonumber\\ & \qquad+\frac{\partial\langle{\phi(\boldsymbol{\theta}(\tau)))}\vert}{\partial\theta_k}\hat{H}\vert{\phi(\boldsymbol{\theta}(\tau))}\rangle+\langle{\phi(\boldsymbol{\theta}(\tau))}\vert\hat{H}\frac{\partial\vert{\phi(\boldsymbol{\theta}(\tau))}\rangle}{\partial\theta_k}, \end{align} \tag{ A.6 }$$

where we use $\langle{\phi(\boldsymbol{\theta}(\tau))}\rangle{\phi(\boldsymbol{\theta}(\tau))} = 1$. Finally, we conclude that

$$\begin{align} \partial\|\left(\partial/\partial\tau+\left[\hat{H}-E_1\left(\tau\right)\right]\right)\vert{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\rangle\|/\partial\dot{\theta}_k = 0 \end{align} \tag{ A.7 }$$

is equivalent to the following linear equation

$$\begin{align} \sum_{m}A_{k,m}\dot{\theta}_m = C_k \end{align} \tag{ A.8 }$$

with

$$\begin{align} A_{k,m}& = \textrm{Re}\left(\frac{\partial\langle{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\vert}{\partial\theta_k}\frac{\partial\vert{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\rangle}{\partial\theta_m}\right) \quad\textrm{and}\quad \nonumber\\ C_k & = -\textrm{Re}\left(\frac{\partial\langle{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\vert}{\partial\theta_k}\hat{H}\vert{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\rangle\right). \end{align} \tag{ A.9 }$$

We note that

$$\begin{align} \frac{\partial\langle{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\vert}{\partial\theta_k}\hat{H}\vert{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\rangle = \sum_\beta\lambda_\beta\frac{\partial\langle{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\vert}{\partial\theta_k}h_\beta\vert{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\rangle, \end{align} \tag{ A.10 }$$

where the h_β are Pauli matrices and λ_β are corresponding coefficients. Hence, QITE reduces to solving the linear equation $\textbf{A}\dot{\boldsymbol\theta} = \textbf{C}$ for $\dot{\boldsymbol\theta}$ which could be accomplished by inversion of the matrix A, as follows $\dot{\boldsymbol\theta} = \textbf{A}^{-1}\textbf{C}$. In our numerical simulations, however, we solve the linear equation via the conjugate gradient method [51] using a subroutine from the SciPy library with a tolerance for convergence of 10⁻⁶. That approach by-passes the need of inverting the matrix A.

A.1.2. Quantum circuit evaluations of A and C

We efficiently evaluate the components of A and C following [18, 52, 53] by implementing the Hadamard test with an additional ancilla qubit. Recall that $\vert{\phi(\boldsymbol{\theta})}\rangle = U_{\mathcal{N}_\theta}(\theta_{\mathcal{N}_\theta})\cdots U_1 (\theta_1)\vert{\bar{0}}\rangle$ and that in the variational ansatz circuit we are only concerned with single- and two-qubit unitary gates $U_n(\theta_n)$, namely rotational or controlled rotational gates. The required derivatives of the ansatz wavefunction are then determined as follows.

Suppose $U_n(\theta_n)$ is a single-qubit rotational gate $R_{\theta_n}^{Z} = e^{-i\theta_n\hat{Z}/2}$, with derivative $\partial U_n(\theta_n)/\partial\theta_n =$ $-(i/2)\times\hat{Z}e^{-i\theta_n\hat{Z}/2}$. If $U_n(\theta_n)$ is a two-qubit controlled rotational gate $\vert 0\rangle\langle 0\vert\otimes\hat{I}+\vert 1\rangle\langle 1\vert\otimes R^{Z} _{\theta_n}$, the derivative is given by $\partial U_n(\theta_n)/\partial\theta_n = \vert 1\rangle\langle 1\vert\otimes\partial R_{\theta_n}^{Z}/\partial\theta_n = (-i/2)\times\vert 1\rangle\langle 1\vert\otimes\hat{Z} e^{-i\theta_n\hat{Z}/2}$. In our numerical simulations, we do not consider parameterized two-qubit gates for simplicity. Using the notation

$$\begin{align} \bar{U}_{n} & = U_{\mathcal{N}_\theta}\left(\theta_{\mathcal{N}_\theta}\right)\cdots U_{n+1}\left(\theta_{n+1}\right)U_n \left(\theta_n\right)g_{n}U_{n-1}\left(\theta_{n-1}\right)\cdots U_2\left(\theta_2\right)U_1\left(\theta_1\right), \end{align} \tag{ A.11 }$$

we conclude that

$$\begin{align} \bar{U}_{n} & = U_{\mathcal{N}_\theta}\left(\theta_{\mathcal{N}_\theta}\right)\cdots U_{n+1}\left(\theta_{n+1}\right)\left(2i\times\partial U_n\left(\theta_n\right)/\partial\theta_n\right) U_{n-1}\left(\theta_{n-1}\right)\cdots U_2\left(\theta_2\right)U_1\left(\theta_1\right). \end{align} \tag{ A.12 }$$

This implies

$$\begin{align} \frac{\partial\vert{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\rangle}{\partial\theta_n} = \left(-i/2\right)\bar{U}_{n}\vert{\bar{0}}\rangle. \end{align} \tag{ A.13 }$$

Consequently, we have

$$\begin{align} A_{k,m} = \textrm{Re}\left(\frac{\partial\langle{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\vert}{\partial\theta_k}\frac{\partial\vert{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\rangle}{\partial\theta_m}\right) = \frac{1}{4}\textrm{Re}\left(\langle{\bar{0}}\vert\bar{U}^\dagger_{k}\bar{U}_{m}\vert{\bar{0}}\rangle\right) \end{align} \tag{ A.14 }$$

and

$$\begin{align} C_k & = \textrm{Re}\left(-\frac{\partial\langle{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\vert}{\partial\theta_k}\hat{H}\vert{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\rangle\right)\nonumber \\ & = -\textrm{Re}\left(\sum_\beta\lambda_\beta\frac{\partial\langle{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\vert}{\partial\theta_k}\hat{h}_\beta\vert{\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)}\rangle\right)\nonumber\\ & = -\frac{1}{2}\textrm{Re}\left(i\sum_{\beta}\lambda_\beta\langle{\bar{0}}\vert\bar{U}_{k}^\dagger\hat{h}_\beta U\vert{\bar{0}}\rangle\right). \end{align} \tag{ A.15 }$$

Since we are evaluating $\textrm{Re}(\langle{\bar{0}}\vert\bar{U}^\dagger_{k}\bar{U}_{m}\vert{\bar{0}}\rangle)$ and $\textrm{Re}(\langle{\bar{0}}\vert\bar{U}^\dagger_{k}\hat{h}_\beta {U}\vert{\bar{0}}\rangle)$, one can implement them on a quantum computer by carrying out the Hadamard tests shown in figure A1.

Zoom In Zoom Out Reset image size

A.1.3. Parameter update

We make use of the Euler method (the first order Taylor series expansion) to update the variational parameters as

$$\begin{align} \boldsymbol{\theta}\left(\tau+\delta\tau\right) \simeq\boldsymbol{\theta}\left(\tau\right)+\dot{\boldsymbol{\theta}}\left(\tau\right)\delta\tau \approx\boldsymbol{\theta}\left(\tau\right)+\xi_\tau\delta\tau, \end{align} \tag{ A.16 }$$

where ξ_τ is the numerical solution to $\textbf{A}(\tau)\dot{\boldsymbol{\theta}}(\tau) = \textbf{C}(\tau)$. One needs to repeat this procedure $\mathcal{N}_T = \tau_{\textrm{total}}/\delta\tau$ times to simulate imaginary-time-like evolution. The difference between the above parameter update and the gradient descent method is that the latter uses

$$\begin{align} \boldsymbol{\theta}\left(\tau+\delta\tau\right) \simeq\boldsymbol{\theta}\left(\tau\right)+\textbf{C}\left(\tau\right)\delta\tau, \end{align} \tag{ A.17 }$$

which only considers information about the average energy at each time step without taking into account of information about the ansatz circuit itself.

As mentioned in the main text, the main ingredient of QIPA is the choice of a suitable cooling function or oracle to quickly reach the optimal solution. Here, we show that oracles defined by the concatenated exponential functions introduced by equation (2) (main text) evolve the initial state according to the generalized Wick-like rotated Schrödinger equation introduced by equation (3) (main text). We apply the chain rule for derivatives to equation (2) (main text) to obtain

$$\begin{align} \frac{\partial}{\partial y}\beta_n(y) = \beta_n(y)b_n\beta^{^{\prime}}_{n-1}(y) & = \beta_n(y)b_n\beta_{n-1}(y)b_{n-1}\beta^{^{\prime}}_{n-2}(y)\nonumber\\ & = \beta_n(y)b_n\beta_{n-1}(y)b_{n-1}\beta_{n-2}(y)\cdots b_1\beta^{^{\prime}}_0(y)\nonumber\\ & = b_n\cdots b_{1}\exp(b_n\beta_{n-1}(y)+\dots+b_1\beta_0(y))\nonumber\\ & = \prod_{k = 1}^n b_k\exp\left(\sum_{k = 1}^n b_k\beta_{k-1}(y)\right). \end{align} \tag{ B.1 }$$

For a given initial wave function $\vert{\psi(0)}\rangle$, let us introduce the auxiliary functions

$$\begin{align} g_{1,n}(\tau) & = U_n(\tau)\vert{\psi(0)}\rangle,\nonumber\\ g_{2,n}(\tau) & = \sqrt{\langle U_n(\tau)\psi(0)\mid U_n(\tau)\psi(0)\rangle}. \end{align} \tag{ B.2 }$$

Note that the normalized time-evolved wavefunction $\vert{\psi(\tau)}\rangle$ can be written as

$$\begin{align} \vert{\psi(\tau)}\rangle & = \frac{g_{1,n}(\tau)}{g_{2,n}(\tau)} = \frac{U_n(\tau)\vert{\psi_0}\rangle}{\sqrt{\langle U_n(\tau)\psi(0)\mid U_n(\tau)\psi(0)\rangle}}. \end{align} \tag{ B.3 }$$

For the derivatives of the above functions we obtain

$$\begin{align} \frac{\partial}{\partial\tau}g_{1,n}(\tau) & = \frac{\partial}{\partial\tau}U_n(\tau)\vert{\psi(0)}\rangle\nonumber\\ & = -\hat{H}U_n(\tau)b_n\beta^{^{\prime}}_{n-1}(-\hat{H}\tau)\vert{\psi(0)}\rangle = -\hat{H} b_n\beta^{^{\prime}}_{n-1}(-\hat{H}\tau)g_{1,n}(\tau),\nonumber \end{align}$$

$$\begin{align} \frac{\partial}{\partial\tau}g_{2,n}(\tau) & = \frac{1}{2}(\langle U_n(\tau)\psi(0)\mid U_n(\tau)\psi(0)\rangle)^{-1/2}\bigg(\left\langle\frac{\partial}{\partial\tau}U_n(\tau)\psi(0)\mid U_n(\tau)\psi(0)\right\rangle\nonumber\\ & \quad + \left\langle U_n(\tau)\psi(0)\left|\frac{\partial}{\partial\tau}U_n(\tau)\psi(0)\right\rangle\right.\bigg)\nonumber\\ & = -\frac{1}{2}g_{2,n}(\tau)^{-1}a_n(\langle\hat{H}\beta^{^{\prime}}_{n-1}(-\hat{H}\tau)U_n(\tau)\psi(0)\mid U_n(\tau)\psi(0)\rangle \nonumber\\ & \quad +\langle U_n(\tau)\psi(0)\mid\hat{H}\beta^{^{\prime}}_{n-1}(-\hat{H}\tau)|U_n(\tau)\psi(0)\rangle),\nonumber \end{align}$$

$$\begin{align} \frac{\partial}{\partial\tau}\vert{\psi(\tau)}\rangle & = \frac{\partial g_{1,n}(\tau)/\partial\tau}{g_{2,n}(\tau)}-\frac{g_{1,n}(\tau)\partial g_{2,n}(\tau)/\partial\tau}{g_{2,n}(\tau)^2}\nonumber\\ & = -\,\hat{H} b_n\beta^{^{\prime}}_{n-1}(-\hat{H}\tau)\vert{\psi(\tau)}\rangle+\frac{b_n}{2}(\langle\hat{H}\beta^{^{\prime}}_{n-1}(-\hat{H}\tau)\psi(\tau)\mid\psi(\tau)\rangle\nonumber\\ & \quad +\langle\psi(\tau)\mid\hat{H}\beta^{^{\prime}}_{n-1}(-\hat{H}\tau)\psi(\tau)\rangle\vert{\psi(\tau)}\rangle) \nonumber\\ & = -b_n(\hat{H}\beta^{^{\prime}}_{n-1}(-\hat{H}\tau)-\textrm{Re}\langle\hat{H}\beta^{^{\prime}}_{n-1}(-\hat{H}\tau)\psi(\tau)\mid\psi(\tau)\rangle)\vert{\psi(\tau)}\rangle\nonumber\\ & = -\,\prod_{k = 1}^{n}a_k\left(\hat{H}\exp\left(\sum_{k = 1}^{n-1} b_k\beta_{k-1}(-\hat{H}\tau)\right)\right.\nonumber\\ & \quad \left.-\textrm{Re}\left\langle\hat{H}\exp\left(\sum_{k = 1}^{n-1} b_k\beta_{k-1}(-\hat{H}\tau)\right)\psi(\tau)\mid\psi(\tau)\right\rangle\right)\vert{\psi(\tau)}\rangle, \end{align} \tag{ B.4 }$$

which is identical to equation (3) (main text). As previously mentioned, the standard QITE is recovered for the choice $n = 1,\,b_1 = 1$ (standard exponential function, cf equation (A.3), as follows:

$$\begin{align} \frac{\partial}{\partial\tau}\vert{\psi(\tau)}\rangle = -(\hat{H}-E_1(\tau))\vert{\psi(\tau)}\rangle,\quad E_1(\tau) = \langle\psi(\tau)\mid\hat{H}\mid\psi(\tau)\rangle. \end{align} \tag{ B.5 }$$

Moreover, for the choice $n = 2,\,b_2 = 1,\,b_1 = 1$ (double-exponential function) we obtain

$$\begin{align} \frac{\partial}{\partial\tau}\vert{\psi(\tau)}\rangle & = -(\hat{H}e^{-\hat{H}\tau}-E_2(\tau))\vert{\psi(\tau)}\rangle,\nonumber\\ E_2(\tau) & = \langle\psi(\tau)\mid\hat{H}e^{-\hat{H}\tau}\mid\psi(\tau)\rangle.\nonumber \end{align}$$

We see that the resulting equation for double-exponential QIPA takes a similar form to the original Wick-rotated equation of motion, but with more rapid convergence to the ground state of $\hat{H}$.

The McLachlan's variational principle is equivalent to the following minimization problem:

$$\begin{align} \delta\left\|\left(\partial/\partial\tau+\left[\hat{H}e^{-\hat{H}\tau}-E_2\left(\tau\right)\right]\right)\vert{\psi\left(\tau\right)}\rangle\right\|^2 = 0 .\end{align} \tag{ C.1 }$$

Since θ is real, we have

$$\begin{align} &\left\|\left(\partial/\partial\tau+\left[\hat{H}e^{-\hat{H}\tau}-E_2\left(\tau\right)\right]\right)|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\rangle\right\|^2\nonumber\\ & \quad = \sum_{m,n}\frac{\partial\langle\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)|}{\partial\theta_m}\frac{\partial|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\rangle}{\partial\theta_n}\dot{\theta}_m\dot{\theta}_n+\langle\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)|^2\left[\hat{H}e^{-\hat{H}\tau}-E_2\left(\tau\right)\right]^2|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\rangle\dots\nonumber\\ & \qquad +\sum_m\frac{\partial\langle\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)|}{\partial\theta_m}\left[\hat{H}e^{-\hat{H}\tau}-E_2\left(\tau\right)\right]|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\rangle\dot{\theta}_m \nonumber\\ & \qquad +\sum_m\langle\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)|\left[\hat{H}e^{-\hat{H}\tau}-E_2\left(\tau\right)\right]\frac{\partial|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\rangle}{\partial\theta_m}\dot{\theta}_m. \end{align} \tag{ C.2 }$$

Differentiating equation (C.2) with respect to $\dot{\theta}_k$, we obtain

$$\begin{align} &\partial\left\|\left(\partial/\partial\tau+\left[\hat{H}e^{-\hat{H}\tau}-E_2\left(\tau\right)\right]\right)|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\rangle\right\|^2/\partial\dot{\theta}_k\nonumber\\ &\quad = \sum_m\left(\frac{\partial\langle\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)|}{\partial\theta_k}\frac{\partial|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\rangle}{\partial\theta_m}+\frac{\partial\langle\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)|}{\partial\theta_m}\frac{\partial|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\rangle}{\partial\theta_k}\right)\dot{\theta}_m\dots\nonumber\\ &\qquad+\frac{\partial\langle\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)|}{\partial\theta_k}\left[\hat{H}e^{-\hat{H}\tau}-E_2\left(\tau\right)\right]|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\rangle\nonumber\\ &\qquad+\langle\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)|\left[\hat{H}e^{-\hat{H}\tau}-E_2\left(\tau\right)\right]\frac{\partial|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\rangle}{\partial\theta_k}. \end{align} \tag{ C.3 }$$

Consequently, since $\langle\phi(\boldsymbol{\theta}(\tau))\mid\phi(\boldsymbol{\theta}(\tau))\rangle = 1$, we conclude that

$$\begin{align} \left\|\left(\partial/\partial\tau+\left[\hat{H}e^{-\hat{H}\tau}-E_2\left(\tau\right)\right]\right)|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\rangle\right\| = 0 \end{align} \tag{ C.4 }$$

is equivalent to the following linear equation

$$\begin{align} \sum_{m} A_{k,m}\dot{\theta}_m = C_k \end{align} \tag{ C.5 }$$

with

$$\begin{align} {A}_{k,m} & = \textrm{Re}\left(\frac{\partial\langle\phi(\boldsymbol{\theta}(\tau)|}{\partial\theta_k}\frac{\partial|\phi(\boldsymbol{\theta}(\tau))\rangle}{\partial\theta_m}\right) \quad\textrm{and},\nonumber\\ {C}_k & = -\textrm{Re}\left(\frac{\partial\langle\phi(\boldsymbol{\theta}(\tau))|}{\partial\theta_k}\hat{H}e^{-\hat{H}\tau}|\phi(\boldsymbol{\theta}(\tau))\rangle\right). \end{align} \tag{ C.6 }$$

Therefore, as in the case of the standard exponential oracle, the McLachlan's principle reduces to a linear system of equations.

C.1. Quantum circuit evaluations of A and C

The matrix A and the vector C can be obtained from the circuits shown in figure C1. As seen from the above equation (C.6), the main difference between QIPA and QITE, in terms of comparing the use of double-exponential and exponential cooling functions, is the presence of $e^{-\hat{H}\tau}$ in C_k . Here, we approximate the exponential by its Taylor series expansion to the second order, as follows:

$$\begin{align} {C}_k & = \textrm{Re}\left(\langle{\bar{0}}\vert\bar{U}_{k}^\dagger\hat{H}U\vert{\bar{0}}\rangle\right)-\tau\textrm{Re}\left(\langle{\bar{0}}\vert\bar{U}_{k}^\dagger\hat{H}^2U\vert{\bar{0}}\rangle\right)\nonumber\\ & \quad +\frac{\tau^2}{2}\textrm{Re}\left(\langle{\bar{0}}\vert\bar{U}_{k}^\dagger\hat{H}^3 U\vert{\bar{0}}\rangle\right)+\mathcal{O}\left(\hat{H}^4\right)\,. \end{align} \tag{ C.7 }$$

Zoom In Zoom Out Reset image size

In [30] we presented a formal proof of convergence for IPA applied to discrete optimization problems with a single global minimum and showed that the number of IPA steps necessary to amplify the amplitude of the global minimum to a value higher than $1/2 = 50\%$ scales logarithmically with the number $n\unicode{x2A7E} 2$ of possible states. More precisely, assuming that the initial wavefunction is given by the uniform superposition over all elements of the finite search space and that the unique maximum of the oracle is given by $\lambda_1\gt0$, while all other values are given by the same value $0\lt\lambda_2\lt\lambda_1$ (here λ can be taken as eigenvalues of the problem Hamiltonian), it was shown that after $k\unicode{x2A7E} 1$ IPA iteration steps, the ratio between the minimum and the maximum amplitude of the evolved density function ρ_k (modulus squared of the wavefunction) is given by

$$\begin{align} \frac{\rho_{k,\min}}{\rho_{k,\max}} = \left(\frac{\lambda_2}{\lambda_1}\right)^k. \end{align} \tag{ D.1 }$$

In the following we compare the ratio $\mu = \lambda_2/\lambda_1$ for the exponential oracle $U_1(\tau) = e^{-\tau\hat H}$ and the double-exponential oracle $U_2(\tau) = e^{e^{-\tau\hat H}}$.

Let us assume that $\hat H = M$ (quantum propagator with infinite mass) and that the potential energy surface M attains its unique global minimum at position $x^*$. The global maximum of the oracle function $U_1(\tau)$ is then given by

$$\begin{align} \lambda_{1,U_1} = e^{-\tau M\left(x^*\right)}, \end{align} \tag{ D.2 }$$

where the index U₁ is used to indicate that we consider the exponential oracle. Since $\lambda_{1,U_1}$ is the unique global maximum and all other values of the oracle are equal, it also follows that

$$\begin{align} \lambda_{2,U_1} = e^{-\tau M\left(x_1\right)}, \end{align} \tag{ D.3 }$$

where $x_1\ne x^*$ is one of the remaining points of the search space. Moreover, for the double-exponential oracle $U_2(\tau)$ we obtain

$$\begin{align} \lambda_{1,U_2} = e^{e^{-\tau M\left(x^*\right)}} = e^{\lambda_{1,U_1}} \quad\textrm{and}\quad \lambda_{2,U_2} = e^{e^{-\tau M\left(x_1\right)}} = e^{\lambda_{2,U_1}}, \end{align} \tag{ D.4 }$$

where the index U₂ is used to indicate that we are now considering the double-exponential. Hence we get

$$\begin{align} \mu_\mathrm{QIPA} = \frac{e^{\lambda_{2,U_1}}}{e^{\lambda_{1,U_1}}} = e^{\lambda_{2,U_1}-\lambda_{1,U_1}} = e^{-\left(\lambda_{1,U_1}-\lambda_{2,U_1}\right)}, \end{align} \tag{ D.5 }$$

which shows that the ratio of the double-exponential oracle decays exponentially in the gap $\lambda_{1,U_1}-\lambda_{2,U_1}$ of the exponential oracle. The analysis presented in [30] therefore yields that the minimum number of required QIPA iterations (with double-exponential oracle) is given by

$$\begin{align} k_\mathrm{QIPA} \unicode{x2A7E}-\frac{\log\left(n-1\right)}{\log\left(\mu_\mathrm{QIPA}\right)} = \frac{\log\left(n-1\right)}{\lambda_{1,U_1}-\lambda_{2,U_1}}, \end{align} \tag{ D.6 }$$

which scales logarithmically with the size of the search space n and inversely with the difference $\lambda_{1,U_1}-\lambda_{2,U_1}$. Recall that this must be compared with the minimum number of required IPA iterations (exponential oracle), given by

$$\begin{align} k_\mathrm{QITE} \unicode{x2A7E}\frac{\log\left(n-1\right)}{\log\left(\lambda_{1,U_1}/\lambda_{2,U_1}\right)}, \end{align} \tag{ D.7 }$$

which is significantly larger if $\lambda_{1,U_1}/\lambda_{2,U_1}\approx 1$ and $\lambda_{1,U_1}-\lambda_{2,U_1}\gg 1$.

Concerning the error of our proposed QIPA as compared to QITE, the difference lies in equation (C.6). To be precise,

$$\begin{align} {C}_{k,\mathrm{QIPA}} = -\textrm{Re}\left(\frac{\partial\langle\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)|}{\partial\theta_k}\hat{H}e^{-\hat{H}\tau}|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\rangle\right), \end{align} \tag{ E.1 }$$

while

$$\begin{align} {C}_{k,\mathrm{QITE}} = -\textrm{Re}\left(\frac{\partial\langle\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)|}{\partial\theta_k}\hat{H}|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\rangle\right). \end{align} \tag{ E.2 }$$

In that sense, the main difference lies in $e^{-\hat{H}\tau}|\phi(\boldsymbol{\theta}(\tau))\rangle$ and $|\phi(\boldsymbol{\theta}(\tau))\rangle$. Since we are taking the Taylor expansion of the exponential term in $C_{k,\mathrm{QIPA}}$, the error can be calculated by the Euclidean distance between the exact imaginary time-evolved state and the approximate one as

$$\begin{align} \epsilon = ||e^{-\hat{H}\delta\tau}\vert{\Psi\left(t\right)}\rangle-\left(1-\delta\tau\hat{H}+\left(\delta\tau\right)^2\hat{H}^2 /2+\mathcal{O}\left(\delta\tau^3\right)\right)\vert{\Psi\left(t\right)}\rangle||, \end{align} \tag{ E.3 }$$

with $||\vert{\psi}\rangle|| = \sqrt{\langle{\psi}\rangle{\psi}}$ being the L² norm. By expanding the norm and using the triangle inequality, we arrive at

$$\begin{align} \epsilon = \sqrt{\Delta^2\delta\tau^2+\mathcal{O}\left(\delta\tau^3\right)} \unicode{x2A7D}\Delta\delta\tau+\mathcal{O}\left(\delta\tau^{3/2}\right), \end{align} \tag{ E.4 }$$

where

$$\begin{align} \Delta^2 & = \langle{\Psi\left(t\right)}\vert\left( \left(1+e^{-\hat{H}\delta\tau}\right)^2/\left(\delta\tau^2\right)+2\left(e^{-\hat{H}\delta\tau}-1\right)\hat{H}/\delta\tau -\left(e^{-\hat{H}\delta\tau}-2\right)\hat{H}^2 \right)\vert{\Psi\left(t\right)}\rangle. \end{align} \tag{ E.5 }$$

The time-evolved wavefunction can be expressed as

$$\begin{align} \left|\psi\left(\tau\right)\right\rangle = \frac{g_{1,n}\left(\tau\right)}{g_{2,n}\left(\tau\right)}, \end{align} \tag{ F.1 }$$

where

$$\begin{align} g_{1,n}\left(\tau\right) & = U_{n}\left(\tau\right)\left|\psi\left(0\right)\right\rangle, \end{align} \tag{ F.2 }$$

$$\begin{align} g_{2,n}\left(\tau\right) & = \sqrt{\left\langle U_{n}\left(\tau\right)\psi\left(0\right)|U_{n}\left(\tau\right)\psi\left(0\right)\right\rangle}. \end{align} \tag{ F.3 }$$

The evolution of the wavefunction is then

$$\begin{align} \frac{\partial\left|\psi\left(\tau\right)\right\rangle }{\partial\tau} & = \frac{\partial g_{1,n}\left(\tau\right)/\partial\tau}{g_{2,n}\left(\tau\right)}-\frac{g_{1,n}\left(\tau\right)\partial g_{2,n}\left(\tau\right)/\partial\tau}{g_{2,n}\left(\tau\right)}, \end{align} \tag{ F.4 }$$

where, without specifying the oracle $U_{n}\left(\tau\right)$, the required partial derivatives are

$$\begin{align} \frac{\partial g_{1,n}\left(\tau\right)}{\partial\tau} & = \frac{\partial}{\partial\tau}U_{n}\left(\tau\right)\psi\left(0\right) \end{align} \tag{ F.5 }$$

$$\begin{align} \frac{\partial g_{2,n}\left(\tau\right)}{\partial\tau} & = \frac{1}{2}\left(\left\langle U_{n}\left(\tau\right)\psi\left(0\right)|U_{n}\left(\tau\right)\psi\left(0\right)\right\rangle\right)^{-1/2} \left(\left\langle\frac{\partial}{\partial\tau}U_{n}\left(\tau\right)\psi\left(0\right)|U_{n}\left(\tau\right)\psi\left(0\right)\right\rangle\right.\nonumber\\ & \quad +\left.\left\langle U_{n}\left(\tau\right)\psi\left(0\right)|\frac{\partial}{\partial\tau}U_{n}\left(\tau\right)\psi\left(0\right)\right\rangle\right)\nonumber\\ & = \frac{1}{2}g_{2,n}^{-1}\left(\left\langle\frac{\partial g_{1,n}\left(\tau\right)}{\partial\tau}|g_{1,n}\left(\tau\right)\right\rangle+\left\langle g_{1,n}\left(\tau\right)|\frac{\partial g_{1,n}\left(\tau\right)}{\partial\tau}\right\rangle\right)\nonumber\\ & = \frac{1}{2}g_{2,n}^{-1}\left(\left\langle\frac{\partial g_{1,n}\left(\tau\right)}{\partial\tau}|g_{1,n}\left(\tau\right)\right\rangle+\left\langle g_{1,n}\left(\tau\right)|\frac{\partial g_{1,n}\left(\tau\right)}{\partial\tau}\right\rangle\right). \end{align} \tag{ F.6 }$$

In the special case that

$$\begin{align} U_{n} = \textrm{sech}\left(\hat{H}\tau\right) = \frac{2}{\textrm{e}^{\hat{H}\tau}+\textrm{e}^{-\hat{H}\tau}}, \end{align} \tag{ F.7 }$$

the required functions are

$$\begin{align} g_{1,n}\left(\tau\right) & = \frac{2}{\textrm{e}^{\hat{H}\tau}+\textrm{e}^{-\hat{H}\tau}}\left|\psi\left(0\right)\right\rangle \end{align} \tag{ F.8 }$$

$$\begin{align} g_{2,n}\left(\tau\right) & = \sqrt{\left\langle\frac{2}{\textrm{e}^{\hat{H}\tau}+\textrm{e}^{-\hat{H}\tau}}\psi\left(0\right)|\frac{2}{\textrm{e}^{\hat{H}\tau}+\textrm{e}^{-\hat{H}\tau}}\psi\left(0\right)\right\rangle }\nonumber\\ & = 2\sqrt{\left\langle\frac{1}{\textrm{e}^{\hat{H}\tau}+\textrm{e}^{-\hat{H}\tau}}\psi\left(0\right)|\frac{1}{\textrm{e}^{\hat{H}\tau}+\textrm{e}^{-\hat{H}\tau}}\psi\left(0\right)\right\rangle} \end{align} \tag{ F.9 }$$

such that the gradient of the first function is

$$\begin{align} \frac{\partial g_{1,n}\left(\tau\right)}{\partial\tau} & = \frac{\partial}{\partial\tau}\textrm{sech}\left(\hat{H}\tau\right)\psi\left(0\right)\nonumber\\ & = \frac{2\hat{H}\left(\textrm{e}^{\hat{H}\tau}-\textrm{e}^{-\hat{H}\tau}\right)}{\left(\textrm{e}^{-\hat{H}\tau}+\textrm{e}^{\hat{H}\tau}\right)^{2}}\psi\left(0\right)\nonumber\\ & = \frac{\hat{H}\left(\textrm{e}^{\hat{H}\tau}-\textrm{e}^{-\hat{H}\tau}\right)}{\textrm{e}^{-\hat{H}\tau}+\textrm{e}^{\hat{H}\tau}}\textrm{sech}\left(\hat{H}\tau\right)\psi\left(0\right)\nonumber\\ & = \frac{\hat{H}\left(\textrm{e}^{\hat{H}\tau}-\textrm{e}^{-\hat{H}\tau}\right)}{\textrm{e}^{-\hat{H}\tau}+\textrm{e}^{\hat{H}\tau}}g_{1,n}\left(\tau\right)\nonumber\\ & = H\tanh\left(\hat{H}\tau\right)g_{1,n}\left(\tau\right) \end{align} \tag{ F.10 }$$

and the gradient of the second function is

$$\begin{align} \frac{\partial g_{2,n}\left(\tau\right)}{\partial\tau} & = \frac{1}{2}g_{2,n}^{-1}\left(\left\langle\hat{H}\tanh\left(\hat{H}\tau\right)g_{1,n}\left(\tau\right)|g_{1,n}\left(\tau\right)\right\rangle + \left\langle g_{1,n}\left(\tau\right)|\hat{H}\tanh\left(\hat{H}\tau\right)g_{1,n}\left(\tau\right)\right\rangle\right)\nonumber\\ & = g_{2,n}^{-1}\textrm{Re}\left(\left\langle g_{1,n}\left(\tau\right)\left|\hat{H}\tanh\left(\hat{H}\tau\right)\right|g_{1,n}\left(\tau\right)\right\rangle\right). \end{align} \tag{ F.11 }$$

The evolution of the wavefunction is then

$$\begin{align} \frac{\partial\left|\psi\left(\tau\right)\right\rangle }{\partial\tau} & = \frac{\hat{H}\tanh\left(\hat{H}\tau\right)g_{1,n}\left(\tau\right)}{g_{2,n}\left(\tau\right)} -\frac{g_{1,n}\left(\tau\right)g_{2,n}^{-1}\textrm{Re}\left\langle\left\langle g_{1,n}\left(\tau\right)\left|\hat{H}\tanh\left(\hat{H}\tau\right)\right|g_{1,n}\left(\tau\right)\right\rangle\right\rangle }{g_{2,n}\left(\tau\right)}\nonumber\\ & = \hat{H}\tanh\left(\hat{H}\tau\right)\left|\psi\left(\tau\right)\right\rangle-\textrm{Re}\left\langle\left\langle\psi\left(\tau\right)\left|\hat{H}\tanh\left(\hat{H}\tau\right)\right|\psi\left(\tau\right)\right\rangle\right\rangle\left|\psi\left(\tau\right)\right\rangle\nonumber\\ & = \hat{H}\tanh\left(\hat{H}\tau\right)\left|\psi\left(\tau\right)\right\rangle-\textrm{Re}\left[\left\langle\hat{H}\tanh\left(\hat{H}\tau\right)\right\rangle\right]\left|\psi\left(\tau\right)\right\rangle\nonumber\\ & = \hat{H}\tanh\left(\hat{H}\tau\right)\left|\psi\left(\tau\right)\right\rangle-\left\langle\hat{H}\tanh\left(\hat{H}\tau\right)\right\rangle\left|\psi\left(\tau\right)\right\rangle\nonumber\\ & = \left(\hat{H}\tanh\left(\hat{H}\tau\right)-E_{3}\left(\tau\right)\right)\left|\psi\left(\tau\right)\right\rangle,\quad \textrm{and}\nonumber\\ E_{3}\left(\tau\right) & = \left\langle\psi\left(\tau\right)\left|\hat{H}\tanh\left(\hat{H}\tau\right)\right|\psi\left(\tau\right)\right\rangle , \end{align} \tag{ F.12 }$$

in direct analogy to the exponential and double-exponential oracles. The McLachlan variational principle corresponding to the resulting minimization problem is

$$\begin{align} \delta\left\Vert\left(\partial/\partial\tau-\left[\hat{H}\tanh\left(\hat{H}\tau\right)-E_{3}\right]\right)\left|\psi\left(\tau\right)\right\rangle\right\Vert^{2} = 0, \end{align} \tag{ F.13 }$$

which, again in direct analogy to the exponential and double-exponential oracles, is equivalent to the linear equation

$$\begin{align} \sum_{m}A_{k,m}\dot{\theta}_{m} = C_{k}, \end{align} \tag{ F.14 }$$

where

$$\begin{align} A_{k,m} & = \textrm{Re}\left(\frac{\partial\left\langle\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\right|}{\partial\theta_{k}}\frac{\partial\left|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\right\rangle }{\partial\theta_{m}}\right), \end{align} \tag{ F.15 }$$

$$\begin{align} C_{k} & = \textrm{Re}\left(\frac{\partial\left\langle\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\right|}{\partial\theta_{k}}H\tanh\left(H\tau\right)\left|\phi\left(\boldsymbol{\theta}\left(\tau\right)\right)\right\rangle\right). \end{align} \tag{ F.16 }$$

The Hamiltonian for a chemical system in the second quantization picture has the following general form

$$\begin{align} \hat{H} = \sum_{ij}h_{ij}\hat{a}^\dagger_i\hat{a}_j+\sum_{ijkl}V_{ijkl}\hat{a}^\dagger_i\hat{a}^\dagger_k\hat{a}_l\hat{a}_j, \end{align} \tag{ H.1 }$$

where $\hat{a}^\dagger_i$ is a creation operator that creates an electron on the ith orbital, $\hat{a}_i$ is an annihilation operator which removes an electron from the ith orbital, and h_ij and V_ijkl are the one-electron and two-electron interaction coefficients, respectively, which are determined for specific systems. The antisymmetric property of electrons is fulfilled by the anti-commutation relation of the creation and annihilation operators $\{\hat{a}_i,\hat{a}^\dagger_j\} = \delta_{ij},\,\{\hat{a}^\dagger_i,\hat{a}^\dagger_j\} = 0$. The above anti-commutation relation precludes direct encoding of a chemical Hamiltonian on a quantum computer, since the operating units of a quantum computer (i.e. qubits) obey the commutation relation of spins. The remedy to this discrepancy is to perform a fermion-spin mapping, such as the Jordan–Wigner transformation [54] presented here as an example. The transformation maps the fermionic creation and annihilation operators to qubit raising and lowering operators ${\hat{\sigma}}^{\pm} = \hat{X}\pm i\hat{Y}$ with a string of $\hat{Z}$ operators to enforce the fermionic anti-commutation properties: $\hat{a}_j \rightarrow(\prod_{l = 1}^{j-1}-\hat{Z}_l)\hat{\sigma}^{-}_j,$ $\hat{a}^\dagger_j \rightarrow(\prod_{l = 1}^{j-1}-\hat{Z}_l)\hat{\sigma}^{+}_j$ With the above transformation, the fermionic anti-commutation relation is preserved. For other fermion-spin mapping approaches, the reader is referred to the literature [4]. An alternative method to decompose gates for molecular systems is provided in [55]. Here the atomic orbital basis (sto3g basis set) is used to ground-state energy is determined at the Hartree–Fock level of theory via PySCF [56].

As the number of high-quality qubits inside a quantum processing unit (QPU) grows over time, it is expected that eventually no classical supercomputer will be able to simulate, verify, and cross-check the inner working mechanism and data obtained from the QPU. This is commonly taken to be as crossing over the 'quantum advantage' threshold. Once such an event occurs, from a practical point-of-view, it is beneficial to make use of existing quantum hardware that is already well-calibrated to simulate subsets of new QPU designs. QCADs of superconducting qubits [33] and photonic chips [34] have recently been proposed and experimentally realized in a superconducting qubit architecture [35], but not yet with imaginary-time-like quantum simulation. Here, we show that, with the proposed QIPA, we are able to optimize for the ground-state-energy search of a flux-tunable superconducting transmon system.

Given an arbitrary classical electrical circuit diagram composed of inductors, capacitors, and Josephson junctions, one can quantize such circuit into a quantum Hamiltonian [57] via the Legendre transformation. Once we obtain the quantum Hamiltonian, the task is to translate it into a language that a quantum computer can understand, such as Pauli words or strings. Let us consider the case of a flux-tunable transmon system shown in the main text figure 1 prior to conversion as an example. The system has the following Hamiltonian:

$$\begin{align} \hat{H}_{\textrm{transmon}} = 2e^2\hat{N}^2/\mathbf{C}_\mathrm{cap}-2E_{\textrm{J}}\left|\cos\left(2\pi f{\,}\right)\right|\cos\hat{\varphi} . \end{align} \tag{ I.1 }$$

Here, e is the electron charge. The normalized external flux $f = \Phi_{\textrm{ext}}/\Phi_0$ is derived from the external magnetic flux $\Phi_{\textrm{ext}}$ that penetrates the loop formed by the two Josephson junctions of the transmon. The Josephson energy of the two junctions is equal and given by E_J while $\mathbf{C}_\mathrm{cap}$ is the capacitance. The magnetic flux quantum Φ₀ is a fundamental constant that describes the smallest amount of flux that a superconducting loop can sustain. Here, $\hat\varphi$ and $\hat{N}$ are the phase and number operators, respectively and fulfill the commutation relation $[\hat{\varphi},\hat{N}] = i$. And, the following relations follow:

$$\begin{align} \left[{e^{i\hat{\varphi}}},\hat{N}\right] = -{e^{i\hat{\varphi}}},\quad {e^{\pm i\hat{\varphi}}}\vert{n}\rangle = \vert{n\pm 1}\rangle, \end{align} \tag{ I.2 }$$

where $\vert{n_j}\rangle$ are the eigenstates of $\hat{N}$. We notice that the operators $e^{\widehat{{\pm i}{\varphi_j}}}$ are similar to the usual bosonic creation and annihilation operators, without the square root prefactor. They are, in fact, the Susskind–Glogower phase operators [58]. When we write down the transmon Hamiltonian in the charge number basis, we use equation (I.2) and assign the operators as:

$$\begin{align} \hat{N} & = \sum_{n = 0}^{d-1}\left(n-\frac{d}{2}\right)\vert n\rangle\langle n\vert, \end{align} \tag{ I.3 }$$

$$\begin{align} \cos\hat{\varphi} & = \frac{1}{2}\sum_{n = 0}^{d-2}\left(\vert n\rangle\langle n+1\vert+\vert n+1\rangle\langle n\vert\right), \end{align} \tag{ I.4 }$$

$$\begin{align} \sin\hat{\varphi} & = \frac{i}{2}\sum_{n = 0}^{d-2}\left(\vert n\rangle\langle n+1\vert-\vert n+1\rangle\langle n\vert\right). \end{align} \tag{ I.5 }$$

In general, the number of Cooper pairs can take on infinitely many integer values. However, for practical purposes, we are only interested in low-lying energy states. In that case, we can truncate the Hilbert space as described above by introducing a finite maximum number of excitations $d = 2^k$. The number of data qubits used in the quantum simulation of the transmon qubit is $k\in\mathbb{N}$.

Next, we convert the eigenstates of the number operator $\hat{N}$ into the computational basis states of k data qubits by representing the integer charge number in a preferred encoding [59–61]. This implies a truncation of the physical space to the subspace spanned by 2^k Cooper pair numbers. There are combinatorially many ways to map such a state space to a set of qubits. For all the numerical quantum simulation experiments presented in this work, we have employed the Gray code [33, 34, 61] due to its resource-efficient representation of tridiagonal quantum matrix operators (table I1).

Table I1. Qubit encodings (standard binary and Gray code) of elementary operators used in this study, with a truncation of d = 16. In our numerical experiments, we utilize the Gray code [61]. Here, $X,Y,Z$ correspond to standard Pauli matrices and the subscript represents the qubit index number.

d = 16	Std. Binary	Gray
N	−0.5 I	−0.5 I
	−4.0 Z3	−4.0 Z3
	−2.0 Z2	−2.0 $Z_{2}Z_{3}$
	−1.0 Z1	−1.0 $Z_{1}Z_{2}Z_{3}$
	−0.5 Z0	−0.5 $Z_{0}Z_{1}Z_{2}Z_{3}$
$\cos\varphi$	+0.5 X0	+0.5 X0
	+0.25 $X_{0}X_{1}$	+0.25 X1
	+0.25 $Y_{0}Y_{1}$	−0.25 $Z_{0}X_{1}$
	+0.125 $X_{0}X_{1}X_{2}$	+0.125 X2
	+0.125 $X_{0}Y_{1}Y_{2}$	−0.125 $Z_{1}X_{2}$
	+0.125 $Y_{0}X_{1}Y_{2}$	+0.125 $Z_{0}X_{2}$
	−0.125 $Y_{0}Y_{1}X_{2}$	−0.125 $Z_{0}Z_{1}X_{2}$
	+0.0625 $X_{0}X_{1}X_{2}X_{3}$	+0.0625 X3
	+0.0625 $X_{0}X_{1}Y_{2}Y_{3}$	−0.0625 $Z_{2}X_{3}$
	+0.0625 $X_{0}Y_{1}X_{2}Y_{3}$	+0.0625 $Z_{1}X_{3}$
	−0.0625 $X_{0}Y_{1}Y_{2}X_{3}$	−0.0625 $Z_{1}Z_{2}X_{3}$
	+0.0625 $Y_{0}X_{1}X_{2}Y_{3}$	+0.0625 $Z_{0}X_{3}$
	−0.0625 $Y_{0}X_{1}Y_{2}X_{3}$	−0.0625 $Z_{0}Z_{2}X_{3}$
	−0.0625 $Y_{0}Y_{1}X_{2}X_{3}$	+0.0625 $Z_{0}Z_{1}X_{3}$
	−0.0625 $Y_{0}Y_{1}Y_{2}Y_{3}$	−0.0625 $Z_{0}Z_{1}Z_{2}X_{3}$
$\sin\varphi$	−0.5 Y0	−0.5 $Y_{0}Z_{1}Z_{2}Z_{3}$
	−0.25 $X_{0}Y_{1}$	−0.25 $Y_{1}Z_{2}Z_{3}$
	+0.25 $Y_{0}X_{1}$	+0.25 $Z_{0}Y_{1}Z_{2}Z_{3}$
	−0.125 $X_{0}X_{1}Y_{2}$	−0.125 $Y_{2}Z_{3}$
	+0.125 $X_{0}Y_{1}X_{2}$	+0.125 $Z_{1}Y_{2}Z_{3}$
	+0.125 $Y_{0}X_{1}X_{2}$	−0.125 $Z_{0}Y_{2}Z_{3}$
	+0.125 $Y_{0}Y_{1}Y_{2}$	+0.125 $Z_{0}Z_{1}Y_{2}Z_{3}$
	−0.0625 $X_{0}X_{1}X_{2}Y_{3}$	−0.0625 Y3
	+0.0625 $X_{0}X_{1}Y_{2}X_{3}$	+0.0625 $Z_{2}Y_{3}$
	+0.0625 $X_{0}Y_{1}X_{2}X_{3}$	−0.0625 $Z_{1}Y_{3}$
	+0.0625 $X_{0}Y_{1}Y_{2}Y_{3}$	+0.0625 $Z_{1}Z_{2}Y_{3}$
	+0.0625 $Y_{0}X_{1}X_{2}X_{3}$	−0.0625 $Z_{0}Y_{3}$
	+0.0625 $Y_{0}X_{1}Y_{2}Y_{3}$	+0.0625 $Z_{0}Z_{2}Y_{3}$
	+0.0625 $Y_{0}Y_{1}X_{2}Y_{3}$	−0.0625 $Z_{0}Z_{1}Y_{3}$
	−0.0625 $Y_{0}Y_{1}Y_{2}X_{3}$	+0.0625 $Z_{0}Z_{1}Z_{2}Y_{3}$

Prime factorization of biprimes is essential to modern Rivest–Shamir–Adleman (RSA) encryption algorithms [62] and is seen as a classic test of the power of quantum computing to address problems that are computationally intensive on classical computers [63–66]. Prime factorization algorithms for NISQ quantum computers are essential to demonstrate the promise of quantum computers for this task and complement current approaches such as Shor's algorithm [1, 63], variational quantum factoring (VQF) [67, 68], exact search [27, 65, 69], QITE [18–20, 66, 70–74], and quantum annealing [64, 75–80]. Recently QITE has been used to identify prime factors via global optimization [66], and here we show such an imaginary-time-like approach can be further accelerated by using QIPA.

To solve the prime factorization problem for a given biprime (product of two prime numbers) $N = q^*\times p^*$, we consider the Hamiltonian

$$\begin{align} \hat H_N\left(q,p\right) = d\left(N;q,p\right)^2,\quad d\left(N;q,p\right) = N-q\times p, \end{align} \tag{ J.1 }$$

defined on the space of prime numbers $q,p\unicode{x2A7D}\sqrt{N}$. This Hamiltonian is a non-negative function that attains its global minimum $\hat H_N(q,p) = 0$ for the unique pair of solutions $q = q^*$ and $p = p^*$. Binary representation [33, 34, 61, 64, 66, 67] of the prime factors yields

$$\begin{align} d\left(N;q,p\right) = N-\sum_{j,k = 0}^Lq_jp_k2^{j+k}, \end{align} \tag{ J.2 }$$

where $(q_L,\dots,q_0),(p_L,\dots,p_0)\in\{0,1\}^L$ are the binary representations of q and p, respectively, and $L = \lfloor\log_2(N/2)\rfloor+1$. Assuming that $p^*,q^*\gt2$ (otherwise N is an even number), we further restrict the search space by setting $q_0 = p_0 = 1$ such that (J.2) can also be written in terms of the combined parameter $\vec x = (q_1,\dots,q_L,p_1,\dots,p_L)$, as follows:

$$\begin{align} d\left(N;\vec x\right) = N-\left(1+\sum_{j = 1}^Lx_j2^j\right)\times\left(1+\sum_{k = 1}^{L}x_{L+k}2^k\right), \end{align} \tag{ J.3 }$$

where $x_l^2 = x_l$ holds for all $1\unicode{x2A7D} l\unicode{x2A7D} 2L$ since $x_l\in\{0,1\}$. Equivalently, in terms of the scaled spin parameters $s_l = 2x_l-1$, we write

$$\begin{align} d\left(N;\vec s\right) & = N-\left(2^L+\sum_{j = 1}^Ls_j2^{j-1}\right)\times\dots \left(2^L+\sum_{k = 1}^{L}s_{L+k}2^{k-1}\right), \end{align} \tag{ J.4 }$$

where $s_l^2 = 1$ holds since $s_l\in\{-1,1\}$. For example, for N = 15 we obtain L = 2 and $d(15;\vec s) = 15-(4+s_1+$ $2s_2)\times(4+s_3+2s_4) = -1-4s_1-8s_2-4s_3-s_1s_3-2s_2s_3-8s_4-2s_1s_4-4s_2s_4$, which gives the Hamiltonian $H_{15}(\vec s) = d(15;\vec s)^2 = 186+48s_1+96s_2+84s_1s_2+48s_3+34s_1s_3+68s_2s_3+32s_1s_2s_3+96s_4+68s_1s_4+$ $136s_2s_4+64s_1s_2s_4+84s_3s_4+32s_1s_3 s_4+64s_2s_3s_4+16s_1s_2s_3s_4$, with twin global minima. Alternatively, with previous information about the prime factors of 15, the integer can be factorized with fewer qubits via the test Hamiltonian [66] $\hat{H} = 196-52\hat{Z}_2-52\hat{Z}_0-56\hat{Z}_2\hat{Z}_0-96\hat{Z}_1-48\hat{Z}_2\hat{Z}_1+16\hat{Z}_0\hat{Z}_1+128\hat{Z}_0\hat{Z}_1\hat{Z}_2$, which features a unique ground state $\vert{011}\rangle$ that corresponds to the correct factorization of the number 15 to 3 and 5. The resulting circuit implementation is given in figure G2.

Zoom In Zoom Out Reset image size