Feedback-Based Quantum Algorithm for Excited States Calculation

Salahuddin Abdul Rahman$^1$, Özkan Karabacak$^2$, Rafal Wisniewski $^1$

1 Automation and Control section, Department of electronic systems, Aalborg University, Aalborg, Denmark

2 Department of Mechatronics Engineering, Kadir Has University,  Istanbul, Turkey

§ ABSTRACT

Recently, feedback-based quantum algorithms have been introduced to calculate the ground states of Hamiltonians, inspired by quantum Lyapunov control theory. This paper aims to generalize these algorithms to the problem of calculating an eigenstate of a given Hamiltonian, assuming that the lower energy eigenstates are known. To this aim, we propose a new design methodology that combines the layer-wise construction of the quantum circuit in feedback-based quantum algorithms with a new feedback law based on a new Lyapunov function to assign the quantum circuit parameters. We present two approaches for evaluating the circuit parameters: one based on the expectation and overlap estimation of the terms in the feedback law and another based on the gradient of the Lyapunov function. We demonstrate the algorithm through an illustrative example and through an application in quantum chemistry. To assess its performance, we conduct numerical simulations and execution on IBM's superconducting quantum computer.

§ INTRODUCTION

Calculating ground states and excited states of Hamiltonians is paramount in various domains, including quantum chemistry, condensed matter physics and combinatorial optimization. It is one of the leading applications of quantum computers where a quantum advantage is expected to be beyond classical algorithms in the noisy intermediate-scale quantum (NISQ) era <cit.>. Variational quantum algorithms (VQAs) are the leading algorithms for NISQ devices that can approximately calculate ground states and excited states <cit.>. These algorithms are tailored to meet the demands of NISQ devices. However, the efficacy of VQAs faces challenges, including the intricate design of ansatz structures and the need to tackle computationally demanding classical optimization problems for parameter updates within parameterized quantum circuits <cit.>. Some of the proposals to tackle the challenge of ansatz design are based on a layer-wise ansatz construction, named as adaptive derivative assembled pseudo-Trotter variational quantum eigensolver (ADAPT-VQE) algorithm <cit.> and adaptive derivative assembled problem tailored quantum approximate optimization algorithm (ADAPT-QAOA)<cit.>. To tackle the challenge of the classical optimization problem,  suitable classical optimizers were proposed <cit.>.

In <cit.>,  Magann et al. introduced feedback-based algorithm for quantum optimization (FALQON) as a novel approach for solving quadratic unconstrained binary optimization (QUBO) problems. FALQON constructs the quantum circuit layer by layer and determines the parameters through measurements of qubits from the preceding layer. This methodology circumvents classical optimization and exhibits a monotonic enhancement of approximate solutions with increasing circuit depth. In addition, the feedback-based quantum algorithm (FQA) was proposed in <cit.> to generalize FALQON to prepare the ground states of Hamiltonians, specifically for molecular Hamiltonians expressed in the second quantization and for Fermi-Hubbard model. In <cit.>, a randomization technique is incorporated to enhance the performance of feedback-based quantum algorithms. FALQON can also be utilized as an initialization technique for the quantum approximate optimization algorithm (QAOA) parameters, potentially increasing its performance <cit.>.

Investigating the excited states of Hamiltonians holds significance across diverse applications, such as quantum chemistry and combinatorial optimization. In quantum chemistry, determining the energy spectrum of a quantum system is crucial, especially when examining drug discovery and catalysis reaction rates. The dynamics of molecules are characterized by their energy spectra; therefore, understanding these spectra is crucial in molecular design <cit.>. Similarly, in combinatorial optimization, calculating excited states is of paramount importance. For example, in problems such as the shortest vector problem (SVP), the solution is encoded as the first excited state of the emerging Hamiltonian, obtained by transforming SVP translated as a QUBO problem into an Ising Hamiltonian <cit.>. In literature, several quantum algorithms have been proposed to find the excited states, such as VQAs <cit.>, imaginary-time evolution algorithms <cit.>, quantum annealing algorithms <cit.> and others <cit.>.

In <cit.>, Higgott et al. introduced the variational quantum deflation (VQD) algorithm, which extends the applicability of the Variational Quantum Eigensolver (VQE) to compute excited states. They utilized Hotelling's deflation method <cit.> to construct an effective Hamiltonian, which possesses its lowest eigenstate as the first excited state of the original Hamiltonian, and its lowest eigenvalue corresponds to the energy of that particular state. Accordingly, they define a new cost function to minimize the energy of the effective Hamiltonian instead of the original Hamiltonian. In this paper, using a similar methodology and utilizing tools from quantum Lyapunov control (QLC) theory, we propose the feedback-based quantum algorithm for excited states preparation (FQAE). We combine the layer-wise construction of the quantum circuit in feedback-based quantum algorithms with a new feedback law based on a new Lyapunov function to assign the quantum circuit parameters. We propose a suitable Lyapunov function using a new operator that encodes the $m$'th excited state of the original Hamiltonian as its ground state. We give two different approaches to evaluate the controllers using both quantum and classical computers. In the first approach, we evaluate the controller by reducing it to different terms, including expected values and overlap terms, and then estimate them using a quantum computer. In the second approach, we utilize the fact that the controller is related to the gradient of the Lyapunov function and use gradient evaluation techniques to evaluate the controller, such as the finite difference approximation method <cit.> and the parameter-shift-rule (PSR) <cit.>. We demonstrate the algorithm through an illustrative example and through an application in quantum chemistry. We assess its performance, convergence properties, and robustness against sampling noise through extensive numerical simulations. Additionally, we demonstrate proof of principle execution on IBM's superconducting quantum computer.

The remainder of this paper is structured as follows. Section 2 reviews QLC and FALQON. Section 3 introduces the proposed FQAE and gives a detailed analysis of how to compute the controller using both classical and quantum computers in a hybrid manner. Section 4 investigates the efficiency of FQAE through an illustrative example and an application in quantum chemistry of computing the spectrum of the hydrogen molecule Hamiltonian. Finally, the outlook and conclusion are given in Section 5.

§ PRELIMINARIES

This section presents a summary of FALQON for preparing ground states and tackling QUBO problems <cit.>. We first examine QLC and then elucidate its relationship with FALQON.

§.§ Quantum Lyapunov Control

Consider the Hilbert space $\mathcal{H}=\mathbb{C}^N$ along with its corresponding orthonormal basis
$\left.\mathcal{A}=\{|q\rangle\right\}_{q \in\left\{0, \ldots, N-1\}\right.}$  and the set of quantum states denoted by $\mathcal{Z}=\{\ket{\psi} \in \mathbb{C}^N : \braket{\psi}=\| \ket{\psi} \|^2=1 \}$. Subsequently, all operators will be represented in the basis $\mathcal{A}$.

Consider a quantum system whose dynamics are governed by the controlled time-dependent Schrödinger equation

\begin{equation}
i|\dot{\psi}(t)\rangle=H(t)|\psi(t)\rangle, \quad H(t)=H_0+H_c(t), \quad H_c(t)=\sum_{l=1}^r H_l u^{(l)}(t),
\label{Model}
\end{equation}

where $\hbar$ is normalized to $1$, $\{u^{(l)}(t)\}_{l=1}^r$
represents the control inputs, and $H_0$ and $H_c$ are the drift (cost) and control (mixer) Hamiltonian, respectively. The Hamiltonians $H_0$ and $H_l$ are assumed to be non-commuting and time-independent: $[H_0, H_l]\neq0$, $l=1,\dots,r$. Assume that for the simplicity of the analysis, we only have one control input in the control Hamiltonian, i.e.,

\begin{equation}
\end{equation}

Hence, the model becomes the following:

\begin{equation}
\label{model1}
\end{equation}

FALQON strives to solve the following control problem.

Problem 1: Given a Hamiltonian $H_0$, the main objective is to find a control law in a feedback form, $u(\ket{\psi(t)})$, that guarantees the convergence of the quantum system (<ref>), from any initial state to the ground state of the Hamiltonian $H_0$, i.e., the state $\ket{\psi_f}=\text{argmin}_{\ket{\psi}\in \mathcal{H}} \bra{\psi}H_0\ket{\psi}$.

We choose the matrix $H_0$ to be diagonal in the $\mathcal{A}$ basis with eigenvalues $E_0, \dots, E_{N-1}$ and corresponding eigenvectors $\ket{E_0}, \dots, \ket{E_{N-1}}$, and introduce our assumptions as follows:

Assumption 1.  The drift Hamiltonian $H_0$ has distinct eigenvalue gaps, i.e., $E_i-E_j \neq E_k-E_l$ for all $(i,j) \neq (k,l).$

Assumption 2.  The control Hamiltonian $H_1$ is fully connected, which means that all off-diagonal elements of $H_1$ are non-zero, i.e., $\bra{E_i} H_1 \ket{E_j} \neq 0$ for all $i \neq j$.

Consider a Lyapunov function of the form

\begin{equation}
V(\ket{\psi}) = \bra{\psi(t)}H_0\ket{\psi(t)}.
\label{V1}
\end{equation}

The derivative of $V(\ket{\psi(t)})$ along the trajectories of system (<ref>) is given by

\begin{align}
\dot{V}(\ket{\psi(t)}) & =  \bra{\psi(t)}  \mathrm{i}[H_1,H_0] \ket{\psi(t)} u(t) .
\label{vdot1}
\end{align}

Designing $u(t)$ as

\begin{equation}
u(t)= - K   g(\bra{\psi(t)}  \mathrm{i}[H_1,H_0] \ket{\psi(t)}),
\label{controller1}
\end{equation}

where $K>0$ and $g$ represents a continuous function that meets the conditions $g(0) = 0$ and $xg(x) > 0$ for all $x \neq 0$, ensures the condition $\dot{V}\leq0$. It is known in the literature (see, e.g., <cit.>) that applying the controller given by (<ref>) to the system (<ref>) guarantees asymptotic convergence of almost all initial states to the ground state of the Hamiltonian $H_0$ given Assumptions 1-2 are satisfied.

§.§ Feedback-Based Quantum Optimization Algorithm

We now demonstrate the connection between FALQON and continuous QLC discussed in the preceding subsection, in line with the content presented in <cit.>.
We start by considering the quantum evolution propagator $U(t)$, which is defined as the solution to Eq. (<ref>) as follows:

\begin{equation}
U(t)=\tau e^{-i \int_{0}^{t} H(s) d s},
\label{propagator}
\end{equation}

where $\tau$ is the time-ordering operator. By decomposing it into $p$ piece-wise constant time intervals of length $\Delta t$, such that the total evolution time $T=p\Delta t$ we obtain

\begin{equation}
U(T, 0) \approx \prod_{k=1}^{p} e^{-i H(k \Delta t) \Delta t},
\label{U(T,0)}
\end{equation}

where $\Delta t$ is chosen sufficiently small such that $H(t)$ remains approximately constant within each interval $\Delta t$. By Trotterization, we can simplify this further to:

\begin{equation}
U(T, 0)  \approx \prod_{k=1}^{p} \bigl(e^{-i u(k \Delta t)  H_{1}\Delta t}  e^{-iH_{0}\Delta t} \bigr).
\end{equation}

Hence, we get a digitized formulation of the evolution in the form

\begin{align}
\ket{\psi_p}   &=\prod_{k=1}^{p} \bigl( e^{-i u(k \Delta t)  H_{1}\Delta t}e^{-iH_{0}\Delta t} \bigr) \ket{\psi_0} =  \prod_{k=1}^{p} \bigl(U_1(u_k)U_0 \bigr) \ket{\psi_0},
\label{evolution}
\end{align}

where we use the notation $u_k=u(k \Delta t)$, $\ket{\psi_k}=\ket{\psi(k \Delta t)}$, $U_0=e^{-iH_{0}\Delta t}$, $U_1(u_k)=e^{-i u(k \Delta t)  H_{1}\Delta t}$. Note that in each discrete time step, the controller's values correspond to the parameters of the circuit, and we will use these terms interchangeably.

The unitary $U_0$ can be efficiently implemented as a quantum circuit given the drift Hamiltonian $H_0$ is specified as a sum of Pauli strings $H_0=\sum_{q=1}^{N_0} c_qO_q$, where $N_0$ is given as a polynomial function of the number of qubits $n$, $\{c_q\}_{q=1}^{N_0}$ are real scalar coefficients,  and a Pauli string is a Hermitian operator in the form $O_q=O_{q,1} \otimes O_{q,2} \otimes ... \otimes O_{q,n}$ with $O_{q,d} \in \{ I,X,Y,Z \}$. Furthermore, to be able to implement unitary $U_1(u_k)$ as a quantum circuit efficiently, the Hamiltonian $H_1$ is designed as a linear combination of Pauli strings as $ H_1 = \sum_{q=1}^{N_1} \Bar{c}_k \Bar{O}_k $, where $N_1$ is a polynomial function of the number of qubits. For such a control Hamiltonian design, the Hamiltonians $H_0$ and $H_1$ are decomposed in terms of a polynomial number of Pauli strings. Therefore, the quantum circuit that implements the evolution defined in (<ref>) can be efficiently implemented on a quantum computer. For a detailed analysis of the implementation of the quantum circuit of unitaries $U_0$ and $U_1(u_k)$, we refer the reader to <cit.>. Figure <ref> shows the circuit implementation of a unitary in the general form of $e^{-iO_{q} \Delta t}$. This circuit can simply be generalized to the case of a linear combination of Pauli strings (see <cit.>).

Quantum circuit implementation for the unitary in the general form $e^{-iO_{q} \Delta t}$, where $O_{q}=O_{q,1} \otimes O_{q,2} \otimes \cdots \otimes O_{q,n}$ is a Pauli string, $U=U_1 \otimes U_2 \otimes \cdots \otimes U_n$ and $U_i=\left\{\begin{array}{ll}
R_y(-\pi/2) , & \text { if } O_{q,i}=X, \\
R_x(\pi/2), & \text { if } O_{q,i}=Y, \\
I, & \text { if } O_{q,i}=Z.
\end{array} \right.$

We have so far approximated the continuous-time evolution by a sequential product of discrete time steps, and the resulting quantum circuit that implements $U(T,0)$ simulates the propagator represented by Eq. (<ref>) where the condition $\dot{V}\leq0$ can be guaranteed to hold by choosing $\Delta t$ to be sufficiently small. For a detailed analysis of how to choose $\Delta t$, see Appendix A in <cit.>. For the control law, a discrete version of the controller given by (<ref>) is adopted as follows:

\begin{equation}
u_{k+1}= -K    \bra{\psi_k} \mathrm{i} [H_1,H_0]  \ket{\psi_k},
\label{udis1}
\end{equation}

where the function $g(\cdot)$ is chosen as the identity function. To evaluate the controller, we need to estimate the term $\bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k}$. Expanding this term by means of Pauli strings, we get

\begin{equation}
\label{expectation}
\bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k} = \sum_{q=1}^{N_2} \hat{c}_q \bra{\psi_k} \hat{O}_q \ket{\psi_k},
\end{equation}

where we employed the product rule of Pauli strings, which asserts that $O^qO^j=\delta_{q,j}I + \mathrm{i} \epsilon_{qjk} O^k$, where $O^1=X$, $O^2=Y$ and $O^3=Z$. Here, $\delta_{qj}$ represents the Kronecker delta, and $\epsilon_{qjk}$ denotes the Levi Civita symbol. Therefore, the value of $\bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k}$ can be determined by estimating the expectations of $N_2$ Pauli strings, where the specific value of $N_2$ is determined by the characteristics of the Hamiltonians $H_0$ and $H_1$. Note that since $N_0$ and $N_1$ are polynomial functions of the number of qubits, $N_2$ will also be a polynomial function of the number of qubits.

The procedure for implementing FALQON is detailed in Algorithm 1 and depicted in Figure (<ref>). The time step is initially set to $\Delta t$, and the procedure is seeded with an initial value of $u_1 = u_{\text{init}}$, which is set to $ u_{\text{init}}=0$. Following this, a group of qubits is initialized to an easy-to-prepare initial state $\ket{\psi_0}$, and a single circuit layer is applied to prepare the state $\ket{\psi_1}$. Subsequently, the controller $u_2$ for the next layer of the quantum circuit is estimated using the quantum computer. Next, an additional layer is added to the circuit, and this process is repeated iteratively for a total of p layers. The resulting iteratively designed quantum circuit to approximate the ground state $\prod_{k=1}^{p} \bigl(U_1(u_k)U_0 \bigr)$, together with its parameters $\{u_k\}_{k=1}^p$, forms the output of the algorithm.

FALQON <cit.>
[1]
$H_0$, $H_1$, $\Delta t$, $p$, $\ket{\psi_0}$
Circuit parameters $\{u_k\}_{k=1}^p$
Set $u_1=0$

for $k = 1$ to $p-1$ do
Prepare the initial state $\ket{\psi_0}$
Prepare the state $\ket{\psi_k}=\prod_{l=1}^{k} \bigl(U_1(u_k)U_0 \bigr)\ket{\psi_0}$
Estimate the term $\bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k}$ by estimating each of the expectations in (<ref>) on the quantum computer
Calculate $u_{k+1}$ using (<ref>)
end for

This figure is adapted from <cit.>. The figure demonstrates the procedural steps involved in FALQON. In the first step, the algorithm is fed with an initial guess of the value of the controller, the state $\ket{\psi_1}$ is prepared in the quantum computer and the controller for the next layer $u_2$ is calculated. Over the following $k$ iterations, the quantum circuit expands by adding a layer composed of $U_1(u_k)U_0$, with $u_k$ calculated from the preceding layer using Equation (<ref>). The process concludes upon reaching the maximum layer count $p-1$.

§ FEEDBACK-BASED QUANTUM ALGORITHM FOR EXCITED STATES CALCULATION

In this section, we introduce the problem of calculating excited states of a given Hamiltonian based on QLC and propose an approach to solving this problem by defining a new operator that encodes the m'th excited state of the original Hamiltonian as its ground state. We then design the control law such that the trajectories converge to this operator's ground state. We demonstrate that employing this method enables an efficient implementation of the quantum circuit and facilitates efficient controller computation.

We start by considering the same model as (<ref>). In this model, we assumed that we have one control input. This assumption can be extended to the case of multiple control inputs, as given in Appendix A. We strive to solve the following control problem.

Problem 2: Given a Hamiltonian $H_0$ with its lower energy eigenstates $\{\ket{\phi^0}, ..., \ket{\phi^{m-1}}\}$ are known, our main objective is to find a control law in a feedback form, $u(\bra{\psi(t)} A \ket{\psi(t)})$ where $A$ is an observable, that assures the convergence of the quantum system (<ref>), from any initial state to the m'th excited state of the Hamiltonian $H_0$.

Note that in Problem 2 we assume the lower energy eigenstates are known.

We assume that Assumption 1-2 hold. Note that the controller's design is constrained to the form of an expected value of an observable $u(\bra{\psi(t)} A \ket{\psi(t)})$ to facilitate efficient evaluation using quantum computers.

We consider a Lyapunov function of the form

\begin{equation}
V(\ket{\psi}) = \bra{\psi(t)}H_0\ket{\psi(t)} + \sum_{j=0}^{m-1} \alpha_j \abs{\braket{\psi(t)}{\phi^j}}^2,
\label{V}
\end{equation}

where $\{\ket{\phi^0}, ..., \ket{\phi^{m-1}}\}$ are the lower energy eigenstates of the Hamiltonian $H_0$, and $\{ \alpha_0, ..., \alpha_{m-1} \}$ are constants that should be chosen to be sufficiently large. More details on choosing these values will be provided subsequently.

(<ref>) can be expressed as

\begin{align}
&= \bra{\psi(t)}H_0\ket{\psi(t)} +  \bra{\psi(t)} \sum_{j=0}^{m-1} \alpha_j \ketbra{\phi^j} \ket{\psi(t)}, \nonumber \\
& = \bra{\psi(t)} \bigl(H_0 + \sum_{j=0}^{m-1} \alpha_j \ketbra{\phi^j} \bigr)\ket{\psi(t)}, \nonumber \\
& = \bra{\psi(t)} P \ket{\psi(t)}.
\label{LyapunovFunc}
\end{align}

Here, we introduced a new Hermitian operator $P$ defined as

\begin{equation}
P := H_0+\sum_{j=0}^{m-1} \alpha_j M_j,
\label{P}
\end{equation}

where $M_j=\ketbra{\phi^j}$ is the projector on the j'th excited state. Note that the operator $P$ commutes with $H_0$, i.e., $[P,H_0]=0$, hence, it is diagonal in the basis $𝒜$, with eigenvalues $p_0,p_1,...,p_N-1$. In addition, given that the parameters $α_j$ are chosen sufficiently large ($α_j > E_m - E_j$, $j=0,1,...,m-1$), the ground state of the operator $P$, i.e., $|ψ_g⟩=argmin_|ψ⟩∈ℋ ⟨ψ|P|ψ⟩$ is the same as the m'th excited state of the Hamiltonian $H_0$. Therefore, we propose a new control problem equivalent to Problem~2 as follows. \\ \\
\textbf{Problem 3:} \textit{Given the Hamiltonian $H_0$ with its $m-1$ lowest energy eigenstates $\{\ket{\phi^0}, ..., \ket{\phi^{m-1}}\}$, find a control law in a feedback form $u(\bra{\psi(t)} A \ket{\psi(t)})$ where $A$ is an observable, that guarantees the convergence of the quantum system \eqref{model1}, from any initial state to the ground state of the operator $H_0 + \sum_{j=0}^{m-1} \alpha_j \ketbra{\phi^j}$.} \\ \\
The derivative of $V(\ket{\psi(t)})$ along the trajectories of system \eqref{Model} is given by
\begin{align}
\dot{V}(\ket{\psi(t)}) & =  \bra{\psi(t)} \mathrm{i} [H_1,P] \ket{\psi(t)} u(t).
\label{vdot}
\end{align}
We design $u(t)$ such that $\dot{V}\leq0$
\begin{equation}
u(t)= - K   g(\bra{\psi(t)}  \mathrm{i}[H_1,P] \ket{\psi(t)}),
\label{controller}
\end{equation}
where $K>0$ and $g$ represents a continuous function that meets the conditions $g(0) = 0$ and $xg(x) > 0$ for all $x \neq 0$.

Applying the controller given by \eqref{controller} to the system \eqref{model1} ensures asymptotic convergence of almost all initial states to the ground state of the operator $P$, given in addition to Assumptions~1-2, the following assumption is also satisfied \cite{cong2013survey,grivopoulos2003lyapunov}:
\\
\textbf{\textit{Assumption 3.}} The operator $P$ has a non-degenerate spectrum, i.e., $p_i \neq p_j$ for all $i \neq j$.

The choice of the shifting parameters $\{\alpha_j\}_{j=0}^{m-1}$ in the design of the operator $P$ is crucial to guarantee that Problem 3 is equivalent to Problem 2. In the following, we give a detailed analysis of the choice of the parameters $\{\alpha_j\}_{j=0}^{m-1}$.

For the operator $P$ defined as in \eqref{P}, and for a general state expressed in the basis $\mathcal{A}$ as $\ket{\psi}=\sum_{q=0}^{N-1} \beta_q\ket{q}$,
we have
\begin{equation}
\label{Pe}
\expval{P}{\psi} = \sum_{q=0}^{m-1}\left|\beta_q\right|^2\left(E_q+\alpha_q\right)+\sum_{q=m}^{N-1}\left|\beta_q\right|^2 E_q,
\end{equation}
where $E_0, E_1, ..., E_{m-1} $ are the eigenvalues of the lower energy eigenstates. To guarantee that the ground state of the operator $P$ is the targeted m'th excited state of the Hamiltonian $H_0$, we need to choose $\alpha_q > E_m-E_q$ for all $q \in \{0, 1, \ldots, m-1\}$. Ideal performance is achieved by choosing these parameters large enough to make Problem~2 and Problem~3 equivalent and avoid excessively high values of $\alpha_q$ as they may adversely affect the algorithm's efficiency.  In our analysis, we assume that the lower energy eigenstates are known, and accordingly, we can find their corresponding eigenvalues. However, the energy of the targeted eigenstate $E_m$ is unknown. Therefore, we cannot choose $\alpha_q$ directly. In the following, we adopt three alternative approaches for choosing the shifting parameters $\alpha_q$ similar to those proposed for designing the cost function for VQD algorithm \cite{higgott2019variational}.

The first approach is to find an upper bound on the energy gap $G$ defined as $G := E_{\max}-E_{\min} \geq E_m-E_q$, where $E_{\max}$ and $E_{\min}$ are the maximum and minimum eigenvalues of $H_0$, respectively. For a drift Hamiltonian $H_0$ specified as a sum of Pauli strings $H_0=\sum_{q=1}^{N_0} c_qO_q$, an upper bound can be found as
\begin{equation}
G \leq 2||H|| \leq 2 \sum|c_q|.
\label{G}
\end{equation}
This upper limit provides insights into selecting the appropriate value for $\alpha_q$.

In the second approach, we possess an accurate estimate of the gap $G$ by evaluating $E_{\min}$ and $E_{\max}$ using a quantum algorithm, where $-H_0$ is used to evaluate $E_{\max}$.

The third approach involves an iterative adjustment of the parameters $\alpha_j$. We initialize $\alpha_q$ with a small guess, and if convergence to the lower energy eigenstate is observed, indicating that $\alpha_q < E_m - E_q$, we double its value and repeat the process. This iterative procedure guarantees $\alpha_q > E_m - E_q$ after $O(\log_2(E_m - E_q))$ repetitions of the algorithm. Although this approach requires multiple evaluations, it can offer improved performance by choosing smaller parameter values.

Following the same Trotterization procedure presented in Subsection \ref{SFQA}, we get the same quantum circuit as in \eqref{evolution}. For the controller, we adopt the following feedback law:
\begin{equation}
u_{k+1}= -K    \bra{\psi_k} \mathrm{i} [H_1,P]  \ket{\psi_k},
\label{udis}
\end{equation}
which represents a discrete version of the controller given by Eq. \eqref{controller}, and the function $g(\cdot)$ is chosen as the identity function.

Since FQAE has the same form of the quantum circuit as FALQON, the construction of the quantum circuit of FQAE follows the same procedures as for FALQON presented in Subsection \ref{SFQA}. However, the controller evaluation is different from that of FALQON. Since the operator $P$ cannot be efficiently expressed in terms of Pauli strings, we will not expand the term $\bra{\psi_k} \mathrm{i} [H_1,P]  \ket{\psi_k}$ in a Pauli strings' basis and therefore we cannot directly evaluate the controller given by \eqref{udis} using the same procedure as in FALQON. In the following subsection, we give different approaches on how to evaluate the controller in a hybrid way using both quantum and classical computers. The details of the algorithmic steps to implement FQAE are given in Algorithm 2.

% %--------------------------------------------------------------
%-------------------- Algorithm ------------------------------
\begin{algorithm} [H]
\caption{FQAE}\label{FQAE}
\begin{algorithmic}[1]
\Input{$H_0$, $H_1$, $\Delta t$, $p$, $\ket{\psi_0}$, $\{\alpha_j\}_{j=1}^{m-1}$}
\Output{Quantum circuit parameters $\{u_k\}_{k=1}^p$}
\State{Design the operator $P$ using Eq. \eqref{P}}
\State{Set $u_1=u_\text{init}$}
\State{\textbf{for} $k = 1$ to $p-1$ \textbf{do}}
\State Prepare the initial state $\ket{\psi_0}$
\State Prepare the state $\ket{\psi_k}=\prod_{l=1}^{k} \bigl(U_1(u_k)U_0 \bigr)\ket{\psi_0}$
\State {Calculate the controller for the next layer $u_{k+1}$ using one of the approaches detailed in Subsection~4.A. }
\State{\textbf{end for}}
\end{algorithmic}
\end{algorithm}
\noindent
\textit{Remark 1:}
Our analysis assumes that the lower energy eigenstates, their eigenvalues, and the circuits that prepare them are known. If the lower energy eigenstates are not known, we can use FQAE in an iterative manner by calculating the lowest energy eigenstate $\ket{E_0}$. Subsequently, utilizing $\ket{E_0}$, we can calculate $\ket{E_1}$. Repeating this process for all lower energy eigenstates, we can evaluate the m'th excited state.\\ \\
\textit{Remark 2:}
In addition to being an easy-to-prepare state, various choices for the initial state could be employed to enhance the algorithm's performance. For instance, in \cite{larsen2023feedback}, the initial state $\ket{\psi_0}$ was selected as the ground state of the control Hamiltonian $H_1$, similar to quantum annealing. Another choice involves utilizing a warm-starting technique such as \cite{egger2021warm}, where the initial state corresponds to the solution of a relaxed combinatorial optimization problem. Alternatively, one could adopt a VQA with a short-depth circuit, using its output as the initial state. We defer a comprehensive analysis of different initial state choices to future work. \\ \\
\textit{Remark 3.} When seeking to compute the ground state of the Hamiltonian $H_0$, since there are no eigenstates with lower energy, according to equation \eqref{P}, we obtain $P=H_0$. This design of the operator $P$ leads to the recovery of FQA. Consequently, FQA can be regarded as a special case of FQAE, wherein the operator $P$ is defined as $H_0$, and the targeted eigenstate is the ground state of the Hamiltonian $H_0$.

\subsection{Controller Computation} \label{cc}
This subsection provides a detailed analysis of the controller's evaluation using a hybrid approach with quantum and classical computers. We propose two alternative methods for evaluating the controller.

\subsubsection{Expectation and Overlap Estimation-Based Approach for Evaluating the Controller}
We introduce our first approach for evaluating the controller. In this approach, we reduce the evaluation of the controller given by \eqref{udis} into several terms, including expected values of the form $\bra{\psi} B \ket{\psi}$ and overlap terms of the form $\bra{\psi} C \ket{\phi}$ where $B$ and $C$ are Pauli strings, and $\ket{\psi}$ and $\ket{\phi}$ are arbitrary quantum states. We estimate these terms on a quantum computer and then use these estimates to compute the controller by post-processing on a classical computer. To evaluate the overlap terms, we use the Hadamard test \cite{aharonov2006polynomial}. The Hadamard test serves as a subroutine in quantum computing that enables estimation of the real and imaginary parts of the inner product $\bra{\psi} U \ket{\psi}$ for an arbitrary quantum state $\ket{\psi}$ and unitary $U$. The quantum circuit for the Hadamard test is depicted in Figure~\ref{HT}.
\begin{figure}[H]
\centering
\includegraphics[ width=0.7 \linewidth]{Figs/HTT.pdf}
\caption{Basic Hadamard test, where $a \in \{0, 1\}$, $H$ is the Hadamard gate, $S=e^{-i\pi Z/4}$ is the phase gate, and $U$ is the unitary we are interested in its expected value with respect to the state $\ket{\psi}$. For $a = 0$, $\expval{Z}$ corresponds to the real part of the inner product $\text{Re}\bra{\psi} U \ket{\psi}$, and for $a = 1$, $\expval{Z}$ corresponds to the imaginary part of the same inner product $\text{Im}\bra{\psi} U \ket{\psi}$.}
\label{HT}
\end{figure}

For the first method, using the controller defined in \eqref{udis} and the operator $P$ defined in \eqref{P} and using the bilinearity property of the commutator, we obtain:
\begin{align}
u_{k+1} & = -K   \bra{\psi_k}  \mathrm{i} [H_1,P]  \ket{\psi_k}, \nonumber \\
& = -K    \bra{\psi_k} \mathrm{i}  [H_1,H_0+\sum_{j=0}^{m-1} \alpha_j M_j]  \ket{\psi_k}, \nonumber \\
& = -K    \Bigl( \bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k} +  \bra{\psi_k}  \mathrm{i} [H_1,\sum_{j=0}^{m-1} \alpha_j M_j]  \ket{\psi_k} \Bigr), \nonumber \\
&= -K    \Bigl( \bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k} + \bra{\psi_k}  \mathrm{i} (H_1 \sum_{j=0}^{m-1} \alpha_j M_j)\ket{\psi_k} - \bra{\psi_k}  \mathrm{i} ((\sum_{j=0}^{m-1} \alpha_j M_j) H_1)  \ket{\psi_k} \Bigr), \nonumber \\
&= -K    \Bigl( \bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k} +\bra{\psi_k}  \mathrm{i} (H_1 \sum_{j=0}^{m-1} \alpha_j M_j)\ket{\psi_k}+(\bra{\psi_k}  \mathrm{i} (H_1 \sum_{j=0}^{m-1} \alpha_j M_j)\ket{\psi_k})^* \Bigr), \nonumber\\
&= -K    \Bigl( \bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k} +2\text{Re}\{\bra{\psi_k}  \mathrm{i} (H_1 \sum_{j=0}^{m-1} \alpha_j M_j)\ket{\psi_k}\}\Bigr),
\nonumber \\
&= -K    \Bigl( \bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k} +2\cdot\text{Re}\{ \mathrm{i}\cdot\sum_{j=0}^{m-1} \alpha_j \bra{\psi_k} H_1 \ket{\phi^j} \braket{\phi^j}{\psi_k}\}\Bigr),
\label{uest}
\end{align}
As outlined in the following steps, we estimate the terms $\bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k}$ , $\bra{\psi_k} H_1 \ket{\phi^j}$ and $\braket{\phi^j}{\psi_k}$ using hybrid quantum and classical processors.
\begin{enumerate}
\item Estimating $\bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k}$: \\
Expanding this term in terms of Pauli strings, we get
\begin{equation}
\label{vk1}
\bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k} = \sum_{q=1}^{N_2} \hat{c}_q \bra{\psi_k} \hat{O}_q \ket{\psi_k}.
\end{equation}
Therefore, by estimating expectations of $N_2$ Pauli strings, we can estimate the term $\bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k}$. The quantum circuit used to estimate this term is shown in Figure~\eqref{sub1}.
\item Estimating $\bra{\psi_k} H_1 \ket{\phi^j}$ for all $j \in \{0,1, ..., m-1\}$: \\
To estimate $ \bra{\psi_k} H_1 \ket{\phi^j}$, we decompose this term as
\begin{align}
\label{cB}
\bra{\psi_k} H_1 \ket{\phi^j} &=  \bra{0} U_k^\dagger (\sum_{q=1}^{N_1} \Bar{c}_q \Bar{O}_q) U^{(j)} \ket{0}= \sum_{q=1}^{N_1} \Bar c_q \bra{0} U_k^\dagger \Bar{O}_q U^{(j)} \ket{0},
\end{align}
where $U_k$ is the circuit to prepare the state $\ket{\psi_k}=U_k \ket{0}$, $U^{(j)}$ is the circuit to prepare the state $\ket{\phi^j}=U^{(j)} \ket{0}$ (See Remark~1).
Each of the terms in Eq. \eqref{cB} is then estimated on the quantum computer using the Hadamard test shown in Figure~\eqref{sub2}.

\item Estimating $\braket{\phi^j}{\psi_k}$ for all $j \in \{0,1, ..., m-1\}$:

The term $\braket{\phi^j}{\psi_k}$ can be estimated using the Hadamard test shown in Figure~\eqref{sub3}, as follows:
\begin{align}
\braket{\phi^j}{\psi_k} &= \bra{0} U^{(j) \dagger} U_k \ket{0}.
\end{align}
\begin{figure}[H]
\centering
\begin{subfigure}{0.9\linewidth} % Adjusted width
\centering
\includegraphics[ width=0.5 \linewidth]{Figs/E1.pdf}

\caption{The quantum circuit used to estimate the term $\bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k}$.}
\label{sub1}
\end{subfigure}

\vspace{1em} % Add some vertical space between the upper and lower figures

\begin{subfigure}{0.45\linewidth} % Adjusted width
\centering
\includegraphics[ width=1 \linewidth]{Figs/O1.pdf}
\caption{The Hadamard test circuit to estimate the term $\bra{0} U_k^\dagger \Bar{O}_q U^{(j)} \ket{0}$. }
\label{sub2}
\end{subfigure}
\hspace{0.05\linewidth} % Adjust horizontal spacing
\begin{subfigure}{0.45\linewidth} % Adjusted width
\centering
\includegraphics[ width=1\linewidth]{Figs/O2.pdf}
\caption{The Hadamard test circuit to estimate the term $\bra{0} U^{(j) \dagger} U_k \ket{0}$. }
\label{sub3}
\end{subfigure}

\caption{The quantum circuits used to estimate the different terms of the controller $u_{k+1}$.}
\label{EO}
\end{figure}

\item After estimating the term $\bra{\psi_k}  \mathrm{i}[H_1,H_0]  \ket{\psi_k}$, and the terms  $\bra{\psi_k} H_1 \ket{\phi^j}$ and $\braket{\phi^j}{\psi_k}$ for all $j \in \{0,1, ..., m-1\}$ using the quantum computer, these values are fed to a classical computer to compute the controller \eqref{uest} as a classical post-processing step. A schematic of the first approach is given in Figure~\ref{FA}.

\begin{figure}[H]
\centering
\includegraphics[ width=1 \linewidth]{Figs/C1.pdf}
\caption{The expectation and overlap estimation-based approach for calculating the controller. }
\label{FA}
\end{figure}

\end{enumerate}
Utilizing several Hadamard tests in this approach poses challenges for near-term quantum devices. In the subsequent subsection, we propose an alternative approach that circumvents the use of Hadamard tests and offers greater feasibility for these devices.

\subsubsection{Gradient-Based Approach for Evaluating the Controller}
As shown in the previous subsection, evaluating the controller involves several Hadamard tests, which require controlled versions of the unitaries that result in deep circuits, making the approach challenging for the current NISQ devices. In this subsection, we introduce an alternative method for evaluating the controller by establishing a connection between the controller and the gradient of the Lyapunov function. We illustrate the relationship between the controller and the Lyapunov function's gradient and propose two evaluation methods: a finite difference approximation and an analytical approach utilizing PSR.

At the k'th step of Algorithm \ref{FQAE}, the algorithm starts from the initial state $\ket{\psi_0}$ and $k$ layers are applied to get the state
\begin{equation}
\ket{\psi_k} = \prod_{l=1}^{k} (U_1(u_l)U_0)\ket{\psi_0} = U_1(u_k)U_0 \prod_{l=1}^{k-1} (U_1(u_l)U_0)\ket{\psi_0}
\label{l-1}.
\end{equation}
From \eqref{LyapunovFunc} and using \eqref{l-1}, we let $\ket{\psi_0}$ and $u_k$ vary and fix $u_1,u_2,...,u_{k-1}$, hence we get
\begin{align}
V_k(\ket{\psi_0},u_k) &\hat = \bra{\psi_k} P \ket{\psi_k}, \nonumber \\
&=        \bra{\psi_0} \Bigl(\prod_{l=1}^{k-1}U_0^\dagger U_1(u_l)^\dagger \Bigr) U_0^\dagger U_1(u_k)^\dagger  P U_1(u_k)U_0 \Bigl(\prod_{l=1}^{k-1} U_1(u_l)U_0 \Bigr)\ket{\psi_0}, \nonumber \\
&= \bra{\psi_{k-1}} U_0^\dagger U_1(u_k)^\dagger P U_1(u_k)U_0 \ket{\psi_{k-1}}, \nonumber \\
&= \bra{\varphi} U_1(u_k)^\dagger P U_1(u_k) \ket{\varphi},
\end{align}
where $\ket{\varphi}=U_0 \ket{\psi_{k-1}}$. \\
For the controlled unitary $U_1(u_k) = e^{-i\Delta t{u_k} H_1}$, the gradient of this unitary with respect to the controller $u_k$ is given by:
\begin{align}
%\nabla_{u_k} U_1(u_k) &= -\mathrm{i} \Delta t H_1 U_1(u_k).
\frac{\partial}{\partial u_k} U_1(u_k) &= -\mathrm{i} \Delta t H_1 U_1(u_k).
\label{dUdu}
\end{align}
Hence, the gradient of $V$ with respect to $u_k$, and using \eqref{dUdu}, is given as:
\begin{align}
\frac{\partial}{\partial u_k} V_k(u_k) &=  \frac{d}{du_k} \langle \phi| U_1^\dagger(u_k) P U_1(u_k) |\phi \rangle, \nonumber \\
&= \mathrm{i} \Delta t \langle \phi| \left[ U_1^\dagger(u_k) H_1 P U_1(u_k) - U_1^\dagger(u_k) P H_1 U_1(u_k) \right] |\phi \rangle, \nonumber \\
&= \mathrm{i} \Delta t \langle \phi| U_1^\dagger(u_k) (H_1 P - P H_1) U_1(u_k) | \phi \rangle, \nonumber \\
&= \mathrm{i} \Delta t \langle \phi| U_1^\dagger(u_k) [H_1,P] U_1(u_k) | \phi \rangle, \nonumber \\
&= \mathrm{i} \Delta t \bra{\psi_{k-1}} U_0^\dagger U_1^\dagger(u_k) [H_1,P] U_1(u_k) U_0  \ket{\psi_{k-1}}, \nonumber \\
&=  \Delta t \bra{\psi_{k}}  \mathrm{i}[H_1,P]  \ket{\psi_{k}}.
\label{dVdu}
\end{align}
From Equations \eqref{udis} and \eqref{dVdu}, we get
\begin{align}
\label{u_dv}
u_{k+1} &= - \frac{K}{\Delta t} \frac{\partial}{\partial u_k} V_k(u_k), \nonumber \\
& = - \frac{K}{\Delta t}( \frac{\partial}{\partial u_k} \bra{\psi_k}H_0\ket{\psi_k} + \frac{\partial}{\partial u_k} \sum_{j=0}^{m-1} \alpha_j \abs{\braket{\psi_k}{\phi^j}}^2 ).
\end{align}

Therefore, to evaluate the controller for the next layer, we need to evaluate $\frac{\partial}{\partial u_k} V_k(u_k)$, which can be evaluated numerically on the quantum computer using the finite difference approximation of the gradient as follows \cite{baydin2018automatic}:
\begin{equation}
\frac{\partial}{\partial u_k} V_k(u_k) \approx \frac{V_k\left(\boldsymbol{u_k}+\epsilon \right)-V_k\left(\boldsymbol{u_k}-\epsilon\right)}{2 \epsilon}.
\label{FDa}
\end{equation}
This means we need to evaluate the Lyapunov function twice to evaluate the controller. A schematic diagram of the gradient-based approach utilizing the finite-difference approximation technique is shown in Figure~\eqref{FD}. The steps for evaluating the controller are given as follows.

\begin{enumerate}
\item Compute $V_k(u_k+\epsilon)$ using the expression:
\begin{align}
V_k(u_k+\epsilon) &= \bra{\psi_k^{+\epsilon}}H_0\ket{\psi_k^{+\epsilon}} + \sum_{j=0}^{m-1} \alpha_j \abs{\braket{\phi^j}{\psi_k^{+\epsilon}}}^2,
\label{V_dis}
\end{align}
where $\ket{\psi_k^{+\epsilon}}=U_1(u_k+\epsilon)U_0 \prod_{l=1}^{k-1} (U_1(u_l)U_0)\ket{\psi_0}=U_k^{+\epsilon}\ket{0}$. The computation of the term $\bra{\psi_k^{+\epsilon}}H_0\ket{\psi_k^{+\epsilon}} = \sum_{q=1}^{N_0} c_q \bra{\psi_k^{+\epsilon}} O_q \ket{\psi_k^{+\epsilon}}$ involves the estimation of $N_0$ expectations. For the term $\sum_{j=0}^{m-1} \alpha_j \abs{\braket{\psi_k^{+\epsilon}}{\phi^j}}^2$, we need to evaluate the overlap between the evolved state and the lower energy eigenstates. While the Hadamard test \cite{aharonov2006polynomial} can be used for this purpose, a more efficient method, as suggested in \cite{havlivcek2019supervised}, involves expressing each overlap term as $ \abs{\braket{\phi^j}{\psi_k^{+\epsilon}}}^2 = \abs{\bra{0} U^{(j) \dagger} U_k^{+\epsilon} \ket{0}}^2$. Consequently, by preparing the state $U^{(j) \dagger} U_k^{+\epsilon} \ket{0}$, the overlap can be estimated by the fraction of all-zero bitstrings obtained from measuring this state in the computational basis. Figure~\ref{Ge} shows the quantum circuits used to estimate both terms.

\begin{figure}[h]
\centering

\begin{subfigure}{0.45\linewidth}
\centering
\includegraphics[width=\linewidth]{Figs/Ge1.pdf}
%         \begin{quantikz}
% \lstick{$\ket{\psi} \: $} &\gate{ U_1(u_k+\epsilon)U_0 \prod_{l=1}^{k-1} \big(U_1(u_l)U_0\big)}\qwbundle[alternate]{} &\meter{}\qwbundle[alternate]{} &
% \end{quantikz}
\caption{The quantum circuit used to estimate the term $\bra{\psi_k^{+\epsilon}}H_0\ket{\psi_k^{+\epsilon}}$. }
\label{subfig:O1}
\end{subfigure}
\hspace{0.001\linewidth} % Adjust horizontal spacing
\begin{subfigure}{0.45\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{Figs/Ge2.pdf}
%         \begin{quantikz}
% \lstick{$\ket{\psi} \: $} &\gate{ \prod_{l=1}^{k} \bigl(U_1(u_k+\epsilon)U_0 \bigr)}\qwbundle[alternate]{} &  \gate{ U^{(j) \dagger}} \qwbundle[alternate]{}& \meter{}\qwbundle[alternate]{} &
% \end{quantikz}
\caption{The quantum circuit used to estimate the term $\abs{\braket{\phi^j}{\psi_k^{+\epsilon}}}^2$.}
\label{subfig:O2}
\end{subfigure}

\caption{The quantum circuits used to estimate the different terms of the controller $u_{k+1}$.}
\label{Ge}
\end{figure}

\item Compute $V_k(u_k-\epsilon)$ using similar procedures to the above by replacing the state $\ket{\psi_k^{+\epsilon}}$ with the state  $\ket{\psi_k^{-\epsilon}}=U_1(u_k-\epsilon)U_0 \prod_{l=1}^{k-1} U_1(u_l)U_0\ket{\psi_0}=U_k^{-\epsilon}\ket{0}$.
\item Compute the controller using Equations \eqref{u_dv} and \eqref{FDa} .
\end{enumerate}

\begin{figure}[H]
\centering
\includegraphics[ width=0.7 \linewidth]{Figs/C2.pdf}
\caption{The gradient-based approach for calculating the controller using finite difference approximation to evaluate the gradient.}
\label{FD}
\end{figure}
The gradient can also be computed using PSR, which gives an analytical formula for evaluating the gradient on a quantum computer \cite{mari2021estimating}. For the case of the cost Hamiltonian $H_1$ having only two distinct eigenvalues, PSR is given as
\begin{equation}
\frac{\partial}{\partial u_k} V_k(u_k) =\frac{1}{2} (V\left({u_k}+\pi \right)-V\left({u_k}-\pi\right)).
\end{equation}

In \cite{wierichs2022general}, PSR is generalized to a more general type of Hamiltonian. However, it will need more computations. A schematic diagram of the gradient-based approach utilizing PSR is shown in Figure~\eqref{PSR}.

\begin{figure}[H]
\centering
\includegraphics[ width=0.7 \linewidth]{Figs/C3.pdf}
\caption{The gradient-based approach for calculating the controller using PSR to evaluate the gradient.}
\label{PSR}
\end{figure}

Estimating the gradient using the finite difference approximation can encounter numerical instabilities and convergence challenges \cite{baydin2018automatic}. Furthermore, in \cite{harrow2021low}, it has been demonstrated that using the analytical gradient outperforms any finite difference method. Therefore, a better solution can be using PSR. However, as shown in \cite{wierichs2022general}, PSR depends on the control Hamiltonian $H_1$, and for a more general Hamiltonian, more computation cost is needed to evaluate the gradient. In contrast, the computational cost of the finite difference method is independent of the specific Hamiltonian $H_1$.

\section{Applications}
To demonstrate the performance of FQAE, we give an illustrative example and apply the algorithm to a quantum chemistry problem where we determine the spectrum of the Hamiltonian of the hydrogen molecule.
\subsection{An illustrative example}
Consider an Ising Hamiltonian in the form
\begin{equation}
H_0=\sum_{q<j}^n J_{q, j} Z_q Z_j+\sum_{q=1}^n J_q Z_q,
\end{equation}
where $n$ is the number of qubits, $J_{q,j}$ and $J_q$ are real parameters with $J_{q,j}=J_{j,q}$, $Z_q$ is the Pauli Z operator applied to the $q^{th}$ qubit.

For our analysis, we set $n=2$, $J_{1,2}=0.5$, $J_1=1$, $J_2=2$. The resulting Hamiltonian is $H_0= Z_1+2Z_2+0.5Z_1Z_2 = \text{diag}(3.5,-1.5,0.5,-2.5)$. We apply FQAE to prepare the first excited state of this Hamiltonian, where the ground state of this Hamiltonian is the state $\ket{\phi^0}=\ket{11}$ and its energy is $-2.5$.

The cost Hamiltonian is designed as $H_c= \sum_{q=1}^2 u^{(q)}H_q =  u^{(1)}H_1+u^{(2)}H_2 = u^{(1)}Y_1+u^{(2)}Y_2$. The initial state is chosen to be the equal superposition state $\ket{\psi_0} = H^{\otimes 2} \ket{0}= \ket{++}$. The quantum circuits representing the operators $U_0$, $U_1(u^{(1)}_k)$, $U_2(u_k^{(2)})$, $\Bar{U}$ and $U^{(0)}$ are shown in Figure (12). %\eqref{QC}.
\newcommand*{\figuretitle}[1]{%
{\centering%   <--------  will only affect the title because of the grouping (by the
\textbf{#1}%              braces before \centering and behind \medskip). If you remove
\par\medskip}%            these braces the whole body of a {figure} env will be centered.

\begin{figure}[H]
\centering
\begin{subfigure}{0.45\linewidth} % Adjusted width
$U^{(0)} \equiv$
\centering
\begin{quantikz}[thin lines]
\qw & \gate{X} & \qw \\
\qw & \gate{X} & \qw
\end{quantikz}
\caption{The quantum circuit to prepare the ground state $\ket{11}=U^{(0)} \ket{00}$.}
\label{fig:subfig1}
\end{subfigure}
\hspace{0.05\linewidth} % Adjust horizontal spacing
\begin{subfigure}{0.45\linewidth} % Adjusted width
\centering
$\Bar{U} \equiv$
\begin{quantikz}[thin lines]
\qw & \gate{H} & \qw \\
\qw & \gate{H} & \qw
\end{quantikz}
\caption{The quantum circuit to prepare the initial state $\ket{\psi_0}=\Bar{U} \ket{0}$.}
\label{fig:subfig2}
\end{subfigure}

\vspace{1em} % Add some vertical space between the upper and lower figures

\begin{subfigure}{0.45\linewidth} % Adjusted width
\centering
$U_0 \equiv$
\begin{quantikz}[thin lines]
\qw & \gate{R_z(2J_{1} \Delta t)} & \ctrl{1} & \qw & \ctrl{1} & \qw  \\
\qw & \gate{R_z(2J_{2} \Delta t)} & \targ{} & \gate{R_z(2J_{1,2} \Delta t)} & \targ{} & \qw
\end{quantikz}
\caption{The quantum circuit to implement the unitary $U_0$.}
\label{fig:subfig3}
\end{subfigure}
\hspace{0.05\linewidth} % Adjust horizontal spacing
\begin{subfigure}{0.45\linewidth} % Adjusted width
\centering
$U_1(u_k^1)U_2(u_k^2) \equiv$
\begin{quantikz}[thin lines]
\qw & \gate{R_y(2u_k^1 \Delta t)} & \qw \\
\qw & \gate{R_y(2u_k^2 \Delta t)} & \qw
\end{quantikz}
\caption{The quantum circuit to implement the cost operator $U_c=U_1(u_k^1)U_2(u_k^2)$.}
\label{fig:subfig4}
\end{subfigure}

\caption{This figure shows the quantum circuits of the different unitary operators needed to implement the different steps of FQAE.}
\label{fig:subfigures4}
\end{figure}

To determine the value of $\alpha_0$, we employ the first approach, calculating the upper bound on the energy gap via equation \eqref{G}, yielding $G \leq 7$, and subsequently assigning $\alpha_0$ as $7$.
The operator $P$ is constructed as $P=H_0 + \alpha_0 \ketbra{11}$ and Hence, we get $P=\text{diag}(3.5,-1.5,0.5,4.5)$. Note that the ground state of the operator $P$ is the state $\ket{\phi^1}=\ket{01}$ which is the first excited state of the Hamiltonian $H_0$. The time step $\Delta t$ is chosen to be $0.08$, and the controllers' gains are chosen to be $K_1 = K_2 =1.5$. To execute the quantum circuit, the Qiskit Aer quantum simulator \cite{aleksandrowicz2019qiskit} is used. The simulation results for preparing the first excited state for the Hamiltonian $H_0$ using the exact statevector simulator are shown in Figure~\ref{SR1}.

\begin{figure}[H]
\centering
\includegraphics[ width=1 \linewidth]{Figs/SS1.pdf}
\caption{Simulation results of FQAE applied to prepare the first excited state $\ket{\phi^1}=\ket{01}$ of the Hamiltonian $H_0$. The layer $k$ is plotted versus: (a) The control inputs, (b) The Lyapunov function of the evolved state $V_k= \bra{\psi_k} P \ket{\psi_k}$, and (c) The fidelity of the evolved state with the eigenstates of the Hamiltonian $H_0$. }
\label{SR1}
\end{figure}
Figure~\ref{SR1} shows that FQAE can calculate the first excited state. The Lyapunov function exhibits a monotonic behaviour, which aligns with our controller design in \eqref{controller}. It is also seen that the fidelity with respect to the first excited state $\abs{\bra{\psi_k} \ket{\phi^1}}^2$ increases with the increase in the circuit layers approaching 1. Note that in practice, evaluating the Lyapunov function at each step is unnecessary. Only an estimation of the controllers is needed to construct the next layer. The convergence of the controllers toward zero can indicate the algorithm's convergence toward the targeted eigenstate, which can be used to set the maximum depth of the quantum circuit.

In practice, the controllers will be evaluated through a finite number of samples. To assess the impact of sampling noise on FQAE's performance, we execute the algorithm for the two controller evaluation approaches proposed in Section~4 and for a different number of sample sizes. The results are compared to the exact values using statevector simulator(i.e. $m \rightarrow \infty$).

\begin{figure}[H]
\centering
\includegraphics[ width=1 \linewidth]{Figs/comparison.pdf}
\caption{Simulation results of FQAE applied to prepare the first excited state of the Hamiltonian $H_0$ for a different number of samples and using the two proposed approaches to evaluate the controller: the expectation and overlap estimation-based approach and the gradient-based approach utilizing PSR. The layer $k$ is plotted versus the mean trajectory (solid line) and the corresponding standard deviation (shaded area) of the controllers $u_k^{(1)}$ and $u_k^{(2)}$ and the Lyapunov function $V_k=\bra{\psi_k}P\ket{\psi_k}$. In addition, the exact noiseless simulation ($m \rightarrow \infty$) is plotted for comparison (black solid line).}
\label{u1u2}
\end{figure}

Figure~\ref{u1u2} shows that even with a small number of samples $m=100$, the estimated controller follows the exact controller, indicating that FQAE is robust against sampling noise. It also shows that the gradient-based approach for evaluating the controller using PSR has better robustness against the sampling noise than the expectation and overlap estimation-based approach.

We conclude this subsection by presenting, as a proof of principle, an experimental demonstration of FQAE on superconducting quantum computers. We run FQAE on the publicly available IBM quantum computer \textit{ibm\_osaka} through the cloud-based service \cite{IBM}. The gradient-based approach using PSR is adopted to evaluate the controller, and the number of shots is set to $4000$. The results are presented in Figure~\ref{experiemental}.

\begin{figure}[H]
\centering
\includegraphics[ width=1 \linewidth]{Figs/experimental.pdf}
\caption{ Experimental results of running FQAE on the IBM superconducting quantum computer \textit{ibm\_osaka}. The gradient-based approach for evaluating the controller is used for this experimental run, where PSR is used to evaluate the controller. The number of shots is set to $4000$.  The layer k is plotted versus: (a) The first control input $u_k^{(1)}$, (b) The second control input $u_k^{(2)}$, and (c) The Lyapunov function $V_k$.}
\label{experiemental}
\end{figure}

From Figure~\ref{experiemental}, we see that the controller's design is robust against noise and can recover the exact controllers calculated through the exact statevector simulator. We also see that the Lyapunov function monotonically decreases toward its global minimum value, which is $\bra{01} P \ket{01} = -1.5$ till the 10th layer. These experimental results underscore the promise of FQAE's practical applicability in the near future.

\subsection{Application to Molecular Hamiltonians}
Using transformations such as the Bravyi-Kitaev or the Jordan-Wigner (JW) transformations, the quantum simulation of fermionic systems can be translated into qubit operations where the Hamiltonian is generally expressed as a sum of Pauli strings as follows (see \cite{helgaker2013molecular} for details).
\begin{equation}
H_M = \sum_q ^{N_m} h_q \hat{O}_q,
\end{equation}
where $h_j$ are real scalar coefficients and $N_M$ is the number of Pauli strings.

As an example, we chose the hydrogen molecule. The fermionic Hamiltonian of hydrogen, modelled using the minimal STO-3G basis and varying inter-nuclear separations, can be converted into qubit representation through the Bravyi-Kitaev transformation. This yields a two-qubit Hamiltonian dependent on the bond length parameter, as follows (see \cite{colless2018computation} for details)
\begin{equation}
H_0(R) = h_0(R) I + h_1(R)Z_1 + h_2(R)Z_2 + h_3(R)Z_1Z_2 +h_4(R)Y_1Y_2+h_5(R)X_1X_2,
\end{equation}
where $R$ is the interatomic distance and the coefficients $h_i(R)$ are real valued functions of the interatomic distance.

For numerical simulations, we use the values of the coefficients given in Table 1 in the supplementary information of \cite{colless2018computation}. The numerical values of the coefficients for the distance value $R=1.05$ are given as $h_0=-0.5626$, $h_1=-0.248783$,  $h_2=-0.248783$, $h_3=0.00850998$, $h_4=0$ and $h_5=0.199984$.

We find the exact eigenstates and their corresponding energies of the hydrogen molecule Hamiltonian using exact diagonalization. We run two different simulations of FQAE to prepare the first and second excited eigenstates of $H_0$ using the statevector simulator. We choose the shifting parameters $\alpha_1 = 1.8$, $\alpha_2 =0.9$, the time step $\Delta t=0.55$, $K_1=K_2=1$ and the initial state to be $\ket{\psi_0}=\ket{01}$. The cost Hamiltonian is designed as $H_c = \sum_{q=1}^2 u^{(q)}H_q = u^{(1)}H_1+u^{(2)}H_2 = u^{(1)} Z_1+u^{(2)}Z_2$. The initial guess for the controllers is chosen to be $[u^{(1)}_0,u^{(2)}_0]=[0,0]$. The simulation results are shown in Figure~\ref{H2}.

\begin{figure}[H]
\centering
\includegraphics[ width=0.95 \linewidth]{Figs/SS2j2.pdf}
\caption{Illustration of FQAE to prepare the first excited state $\ket{E_1}$ (to the left) and the second excited state $\ket{E_2}$ (to the right). The layer $k$ is plotted versus: (a) The control inputs, (b) The Lyapunov function of the evolved state and (c) The fidelity with the targeted eigenstate.}
\label{H2}
\end{figure}

It is seen in Figure~\ref{H2} that FQAE can efficiently calculate the first and second excited states of the Hamiltonian $H_0$. %In a depth of 15 layers, a fidelity of higher than 0.8 is achieved.
The Lyapunov function monotonically decreases over circuit layers and converges to the operator $ P$'s ground state, which is the targeted eigenstate of $H_0$.

In addition, we run FQAE to calculate the energy of the first and second excited states for a range of R values $R \in [0.1,3.95]$ and for a depth of 25 layers. In this case, we design the control Hamiltonian to be $H_c = \sum_{q=1}^4 u^{(q)}H_q = u^{(1)}H_1+u^{(2)}H_2+ u^{(3)}H_3 + u^{(4)}H_4 = u^{(1)} Y_1+u^{(2)}Y_2 + u^{(3)}Z_1 + u^{(4)}Z_2$. For the values of $R \in [0.1,3.95]$ we set $\Delta t=0.55$, while for $R \in [0.1,0.45]$ we set $\Delta t= 0.15$ since these values of $\Delta t$ ensure controllers' convergence as suggested by simulation results. Since the Hamiltonians $H_0$ and $H_c$ have non-commuting terms, we use the first-order Lie-Trotter Suzuki decomposition to implement the unitaries $U_0$ and $U_c$ as quantum circuits in our simulations. The simulation results are shown in Figure~\ref{H2_all}. Note that the ground and highest excited states can be calculated using FQA.
\begin{figure}[H]
\centering
\includegraphics[ width=0.8 \linewidth]{Figs/H2_all2.pdf}
\caption{Demonstration of FQAE to calculate the energy of the first excited state $E_1$ and the second excited state $E_2$ of different values of $R$. The Internuclear distance $R$ is plotted versus the energy calculated using FQAE, with solid lines indicating the exact energy levels of the excited states.}
\label{H2_all}
\end{figure}

\section{Conclusion and Future Work}
In this paper, we have proposed FQAE, which extends FQA for calculating excited states. A comprehensive investigation of the algorithm's efficient implementation and controller evaluation is provided, presenting two alternative approaches for the controller evaluation. Among these approaches, the gradient-based method utilizing PSR yielded the most promising results in terms of performance and computation cost, emerging as the preferred option for near-term applications. The effectiveness of FQAE is demonstrated through an illustrative example and an application in quantum chemistry. Performance evaluations, including numerical simulations and execution on IBM's superconducting quantum computer, confirm the algorithm's potential for practical quantum computations.

As highlighted in Section~2, a significant challenge of the feedback-based quantum algorithms lies in the stringent assumptions imposed on the Hamiltonians $H_0$ and $H_1$, making their practical fulfilment challenging. However, as shown in \cite{magann2022feedback, magann2022lyapunov,larsen2023feedback,magann2023randomized}, and through numerical simulations, asymptotic convergence towards the desired eigenstate and good performance remain achievable even without strict adherence to these assumptions. Additionally, to enhance algorithmic performance, a range of heuristics has been proposed in \cite{magann2022feedback}, offering pathways to achieve convergence towards the target eigenstate in scenarios where adherence to the assumptions is unattainable. Furthermore, the randomization method proposed in \cite{magann2023randomized} helps in convergence to the targeted eigenstate and avoids convergence to sub-optimal solutions. Future endeavours may involve relaxing some of these assumptions, perhaps by incorporating mid-circuit measurements, as explored in \cite{clausen2023measurement}.

Another main challenge of FQAE, which is the bottleneck of the algorithm for being applicable for near-term devices, is the deep circuit required to run the algorithm. One approach to tackle this problem could be, as highlighted in Remark 2, by adopting a warm starting technique for the initial state of the algorithm. Another approach to overcome this challenge could be adopting a fixed-time control strategy as presented in \cite{li2022lyapunov}. The fixed-time Lyapunov control technique could guarantee convergence to the excited state in a fixed time and independently of the initial state. This could help in designing fixed-depth quantum circuits. Therefore, valuable future work would be adapting FQAE to such control designs.

We highlight that feedback-based quantum algorithms remain applicable in near-term devices as a warm-starting technique for VQAs. In \cite{magann2022lyapunov}, it was demonstrated that employing FQA as an initialization technique for QAOA can potentially improve its performance. Our work demonstrates the potential of quantum control theory to inspire the development of efficient quantum algorithms.

\section{Appendix A: Multiple control inputs}
In this appendix, we extend the analysis for the case of multiple control inputs, i.e., $H_c(t)=\sum_{l=1}^r H_l u^{(l)}(t)$ for $r \geq 2$. For this case, the model is given as:
\begin{equation}
i|\dot{\psi}(t)\rangle=(H_0+\sum_{l=1}^r H_l u^{(l)}(t))|\psi(t)\rangle,
\label{model_mult}
\end{equation}
where $\{u^{(l)}(t)\}_{l=1}^r$ are the control inputs each associated with $H_l$.

The derivative of $V(\ket{\psi(t)})$ along the trajectories of system \eqref{model_mult} then becomes:
\begin{equation}
\dot{V}(\ket{\psi(t)}) = \sum_{l=1}^r  \bra{\psi(t)} \mathrm{i} [H_1,P] \ket{\psi(t)} u^{(l)}(t) .
\end{equation}
Designing the control laws $\{u^{(l)}(t)\}_{l=1}^r$ such that $\dot{V}\leq0$ yields
\begin{equation}
u^{(l)}(t)= - K_l  g( \bra{\psi(t)}  \mathrm{i}[H_1,P] \ket{\psi(t)} ),
\label{controllers}
\end{equation}
where $K_l>0$ and $g$ is a continuous function satisfying $g(0) = 0$ and $xg(x) >
0$, for all $x \neq 0$.

As shown in \cite{cong2013survey}, the application of the control laws \eqref{controllers} to the system \eqref{model_mult} guarantees the convergence, starting from any initial state to the ground state of the operator $P$, given the same assumptions in Section~1 hold and Assumption~2 being relaxed as follows. \\
\textbf{\textit{Assumption 2 (for multiple control inputs).}} There exists at least an $l \in \{1,2, ..., r\}$ such that $\bra{E_i} H_l \ket{E_j} \neq 0$.

Following a Trotterization procedure as presented in Subsection~2, we get the following modification of \eqref{evolution}
\begin{align}
\ket{\psi_p}    &=  \prod_{k=1}^{p} U_c(u^{(1)}_k, u^{(2)}_k,...,u^{(r)}_k )U_0\ket{\psi_0},
\label{evolution_mult}
\end{align}
where $U_c(u^{(1)}_k, u^{(2)}_k,...,u^{(r)}_k)=e^{-i \Delta t \sum_{l=1}^r u^{(l)}(k \Delta t)  H_{l}}$.

Accordingly, the feedback law \eqref{udis} is modified to
\begin{equation}
u^{(l)}_{k+1}= -K_l    \bra{\psi_k} \mathrm{i} [H_l,P]  \ket{\psi_k},
\label{udis_mult}
\end{equation}
where $K_l > 0$, $l \in \{1,2, ..., r\}$.

\begingroup
\footnotesize
\raggedright
\bibliography{litterature.bib}
\endgroup

\end{document}