Time Evolution of Schrödinger’s Equation Using Quantum Phase Estimation and Non-Self-Adjoint Hamiltonian

1. Introduction

This paper demonstrates the potential of the scalable QPE algorithm to capture time evolution of the phase derived by Schrödinger’s equation. Notably, this paper introduces a novel method for phase estimation when we are dealing with a non-self-adjoint Hamiltonian. This innovation holds the dimension of the original Hamiltonian unchanged and brings significant promises for tackling intricated quantum systems. We begin by reviewing essential concepts such as modal decomposition in quantum computations, Quantum Fourier Transform (QFT) required to implement QPE, and the inherent probabilistic nature of quantum measurement described by von Neumann measurement.

Quantum computers have the potential to revolutionize and solve longstanding computational problems in cryptography, search, optimization, communication, and simulation of quantum systems ​[1,2,3]​. Shor's algorithm for prime factorization ​[4,5]​, Grover’s algorithm for unstructured search ​[6]​, and the HHL algorithm ​[3]​, which is highly relevant to finite element method (FEM) in engineering, are a few notable quantum algorithms. However, unlike classical computing machines, quantum computer hardware is still nascent, and suffers for performing accurate computation of elementary arithmetic operations. Gehlot et al. ​[7]​ presented the quantum-classical architecture for solving linear differential equation using a proper encoding of an arbitrary state vector. Solving linear differential equations enables the simulation of the linear state observer discussed in ​[8]​ using Positive Operator Valued Measure (POVM) of a closed quantum system. The authors in ​[1]​ presented a table comparing the runtime of classical algorithms for Linear System Problem (LSP) against their quantum counterparts. They indicated that the number of operations required by the best classical algorithm, the Conjugate Gradient (CG) method, is higher than HHL, the slowest quantum algorithm.

The quantum Fourier transform is an efficient quantum algorithm for computing the Fourier transform and it is pivotal in the quantum phase estimation algorithm to estimate eigenvalues of a unitary operator. The development of QFT has been explained in detail in many studies such as ​[6,9,10]​. These studies note that QFT generally does not accelerate the computation of the Fourier transform over classical methods. QFT is the basis for another quantum subroutine algorithm called quantum phase estimation. Briefly, in QPE, a word describing a binary fraction is encoded as the phase of the unitary operator. It is important to note that the length of the binary fraction adjusts to the desired accuracy of the estimated phase.

The role of phase estimation shows up when discussing the eigenvalues of a linear system. By reformulating a dynamical computational model into an eigenvalue problem, QPE has the potential to solve widespread and multi-domain problems. A quantum computational problem exploits QPE to extract eigenvalues of the observer associated with a physical system. References ​[9,11]​ discuss the detailed derivation of QFT for phase estimation. Moreover, they provided a detailed illustration of the steps to recover the phase of a unitary matrix using a single ancilla qubit circuit.

1.1. Eigenvalue

The mathematical framework of a huge number of physical problems ends up with either Ordinary Differential Equations (ODE) or linear system problems. For instance, the governing equations of mechanical systems typically result in an ODE or a Partial Differential Equations (PDE) model that live in the field of reals. On the other hand, Schrödinger’s equations, that describe the dynamics and quantum probability for a closed system, result in linear ODE equations that reside in complex field. Besides, Finite Element Models (FEM), that form the basis of today’s engineering problems, ultimately result in solving large scale systems of linear equations. Diagonalization is a practical approach to solve these problems. In the quantum paradigm, the state sits in the separable Hibert space $\mathcal{H}$, meaning that any arbitrary vector state can be written in orthonormal bases ​[12,13,14]​.

Let’s consider a self-adjoint matrix $A\in {\mathbb{C}}^{N\times N}$. There exist orthonormal eigenvectors $v_i\in {\mathbb{C}}^N$ associated with eigenvalues ${\lambda }_i\in \mathbb{C}$ such that $A_{N\times N}={\sum }_i^N{\lambda }_i|v_i⟩⟨v_i|$ where $N\in \mathbb{N}$ is reserved for the dimension of matrix $A$. Suppose linear systems can be written as $A\left|x\right.⟩=|b⟩$, after writing matrix $A$ and vector $\left|b\right.⟩\in {\mathbb{C}}^N$ in orthonormal bases, this problem can be rearranged as $\left|x\right.⟩={\sum }_i^N{\lambda }_i^{-1}|v_i⟩⟨v_i|\left|b\right.⟩$[3]. For example, consider the first order time-invariant (LTI) system like $\dot{x}=Ax$. The solution would be $x\left(t\right)=e^{At}{x}_0$, where matrix $A$ is an $N\times N$ self-adjoint matrix and $x_0$is the initial states at time $t_0$. It can be written as $\left|x\left(t\right)\right.⟩={\sum }_i^N{e}^{{\lambda }_i}\left|v_i\right.⟩\left.⟨v_i\right|\left|x_0\right.⟩,$using quantum computing notation. This shows that diagonalization (modal decomposition), a typical method to solve challenging computational problems, is an approach used extensively in quantum computing.

1.2. Measurement and Simulation

According to the postulates of quantum mechanics, the time evolution of a closed quantum system is described by a unitary operator $U\left|\Psi \right.⟩=e^{-j2\pi \theta }|\Psi ⟩$(here,$\theta$ is the phase), which preserves the norm of the state vector and ensures that quantum processes are reversible in principle [9,13,14,15].

Measurement in quantum mechanics, specifically the von Neumann measurement (also known as projective measurement), describes the process by which observing a quantum system causes its state to collapse onto one of the eigenvalues of the measured observable. Quantum measurement presents a unique challenge. There exists an ensemble of the realizations of the system under the identical set of conditions leading to the preparation of a single quantum mechanical system. Therefore, what is measured is the expected value of the observable [15]:

$${⟨A⟩}_{\mathsf{\Psi }}=⟨\mathsf{\Psi }\left|A\right|\mathsf{\Psi }⟩={\displaystyle \sum _i^N}{\lambda }_i|{\left.{\alpha }_i\right|}^2\ni \mathsf{\Psi }=\sum _i^N{\alpha }_i{v}_i$$

(1)

where $\Psi$ is a state and $v_i$ and ${\left|{\alpha }_i\right|}^2$ are the eigenvector of the observable $A$ and the quantum probability at which the system collapses, respectively. ${⟨A⟩}_{\Psi }$ addresses that if system sits at state $\Psi$, measurement leads to one of the eigenvalues of the matrix $A$.

1.3. Dynamics of Quantum Particles

The Schrödinger’s equation for a closed system (i.e., a system that does not exchange energy with the environment) is described by

$$j\hslash {\displaystyle \frac{\partial \left|\mathsf{\Psi }\right.⟩}{\partial t}}=H\left|\mathsf{\Psi }\right.⟩$$

(3)

where $H$is a self-adjoint Hamiltonian of the system indicating the total energy of the system. The self-adjoint operator ensures that the energy eigenvalues are real, $\hslash$ is Plank’s constant, and $|\Psi ⟩$ represents the quantum state vector belonging to Hilbert space $\mathcal{H}$. Schrödinger’s equation describes how a quantum state evolves over time in a closed system. Given an initial state and the Hamiltonian, the equation determines the future states of the system. The solution of this linear Schrödinger’s equation is

$$\left|\mathsf{\Psi }\right.⟩=e^{-j{\displaystyle \frac{Ht}{\hslash }}}\left|{\mathsf{\Psi }}_0\right.⟩$$

(2)

where ${\Psi }_0$ is the initial state of the system at time $t_0$.

If the system is not closed, it will interact with the surrounding thermal bath, and additional terms will be present in the equation describing the evolution of the state over time. One representation is the Lindblad master equation [14], which contains additional perturbation terms with a Hamiltonian representing the environment and plays an important role in diverse fields such as quantum optics and quantum information. However, in this paper, we only deal with closed systems.

Time Discretization

Digital quantum computers operate as discrete, digital computation machines. Consequently, to implement qubit dynamics in a quantum circuit, we must discretize the qubit's Ordinary Differential Equation (ODE). Sample time is defined as $\Delta t$. Under these conditions, we can derive the discrete time model of the qubit dynamics as

$$\left|\mathsf{\Psi }\left(t_{k+1}\right)\right.⟩=e^{-j{\displaystyle \frac{H\mathsf{\Delta }t}{\hslash }}}\left|\mathsf{\Psi }\left(t_k\right)\right.⟩\ni \mathsf{\Delta }t=t_{k+1}-t_k$$

(4)

where $t_k$ indicates time at step $k$. This model allows us to represent the continuous evolution of the qubit state in discrete time steps, which is essential for digital quantum computation. The resulting discrete time model provides a framework for simulating and controlling qubit behavior in quantum circuits, bridging the gap between the continuous nature of quantum mechanics and the discrete operations of quantum computers (We assume $\hslash$ = $1$ for simplicity and this does not impact any results presented in this paper).

1.4. Paper Organization

In section 2, we discuss quantum phase estimation and how to implement the phase kick-back algorithm for QPE in subsection 2.1. Subsection 2.2 addresses the post-processing needed to convert the output of the QPE quantum circuit into a phase in the decimal system, incorporating classical computations. Section 3 is dedicated to an illustrative example and simulation results. Before presenting the examples, subsection 3.1 defines the conventions for quantum computing operations in Qiskit. Subsection 3.2 contains a numerical example of QPE, demonstrating how to obtain the phase of a given matrix, besides demonstrating the scalability of the QPE algorithm, and perform post-processing for all phases of the given matrix. Additionally, we present a histogram related to all estimated phases of the example. Section 4 focuses on validating and verifying the results obtained by QPE, including the derivation of phases implementable through classical computation. The next section covers the time evolution of the Schrödinger equation with a self-adjoint Hamiltonian, explaining the phase at different time steps of the given Hamiltonian in the Schrödinger equation. So far, all discussions have been about self-adjoint Hamiltonians. In section 6, we propose an efficient method to handle QPE for a non-self-adjoint Hamiltonian. In contrast to the HHL algorithm, which doubles the dimensionality, this method operates with the original dimension of the Hamiltonian matrix. The final section is the conclusion and future work.

2. Quantum Phase Estimation Algorithm with a Self-Adjoint Hamiltonian

Phase estimation is widely applied to address a range of problems, including hidden subgroup identification, graph isomorphism, quantum walk analysis, quantum sampling, adiabatic computing, order finding, and large number factorization. The implementation of QPE has been described in numerous studies [11,16,17]. This widespread adoption and continuous refinement underscore QPE's status as a cornerstone technique in the field of quantum computation. Following section is dedicated to extract the phase of the unitary operator:

$$U\left|\mathsf{\Psi }\right.⟩=\lambda \left|\mathsf{\Psi }\right.⟩=e^{j2\pi \theta }\left|\mathsf{\Psi }\right.⟩\ni U^{-1}=U^{*}\&0\le \theta \le 1\&\left|\lambda \right|=1$$

(5)

Here, $U$ is a unitary operator with complex eigenvalues of magnitude $1$. The parameter $\theta$ is known as the phase which is supposed to be extracted by the phase kick-back algorithm below. Since this exponential function is periodic, its magnitude is $1$.

2.1. Phase Kick-Back Algorithm for QPE for Self-Adjoint Hamiltonian

Figure 1 shows quantum circuit behind the QPE algorithm. Here, $m$quantum registers initiate with state $\left|0\right.⟩,$with vector representation as $\left[\begin{array}{c}1\\ 0\end{array}\right]$. The number of these ancilla qubits determines the accuracy of the estimated phase. To represent in binary fraction, the phase would be $\theta =0.{\theta }_1{\theta }_2\dots {\theta }_m$, where ${\theta }_i$ is the binary digit that can be either $0$ or $1$. In other words, $m$ determines the number of the binary digits in estimated phase, $\theta$. Last qubits, represented by $|v⟩$, carry the eigenvector associated with the unitary operator $U$. If $N=2^n$ and the dimension of the unitary operator is $N\times N$, then $n$qubits are required to represent the state $|v⟩$ in the circuit. A Hadamard unitary gate takes each state of ancilla qubit to its superposition state, which can be interpreted as a representation in the frequency basis. Then, there is a controlled unitary gate, shown as $U^{2^0}$, $U^{2^1}$, and so on, on the way of each qubit which generates the associated binary digit. The Inverse Quantum Fourier Transform (IQFT) then processes this kicked-back phase to produce the final estimation. The term 'phase kick-back' is often associated with this QPE algorithm because the controlled-U operations 'kick back' the phase information onto the ancilla qubits. It should be mentioned that $R_i$ represents the controlled rotation matrix

$$R_i=\left[\begin{array}{cc}1& 0\\ 0& e^{j{\displaystyle \frac{2\pi }{2^i}}}\end{array}\right],R_i^{-1}=\left[\begin{array}{cc}1& 0\\ 0& e^{-j{\displaystyle \frac{2\pi }{2^i}}}\end{array}\right]$$

(6)

where $R_i^{-1}$ is inverse of the controlled rotation matrix.

It is important to note that the qubits representing the phase are ordered from least to most significant, with ${\theta }_m$ being the least significant qubit and ${\theta }_1$ being the most significant qubit. In other words, ${\theta }_1$ is the first digit number after the decimal point in binary representation of $\theta$and ${\theta }_m$ is the farthest digit from decimal point.

2.2. Post Processing

The result of executing the QPE circuit, for the $1024$ of runs, produces a histogram plot which illustrates the probability distribution of all possible measured states. The next step is to convert the highest probable state from binary representation to decimal system to interpret the results. The binary representation $0.{\theta }_1{\theta }_2\dots {\theta }_m$ converts to a decimal value as follows:

$$0.{\theta }_1{\theta }_2\dots {\theta }_m={\displaystyle \frac{{\theta }_1}{2}}+{\displaystyle \frac{{\theta }_2}{2^2}}+\dots +{\displaystyle \frac{{\theta }_m}{2^m}}\ni {\theta }_i=0\mathrm{o}\mathrm{r}1$$

(7)

The decimal value representation of the phase falls within the range $0\le \theta \le 1$.

3. Example and Simulation Results

QASM, the quantum simulator provided by Qiskit (an opensource software), paves the way to run a quantum algorithm on the local computer. Simulator at the level of circuit governed by quantum computation rules lets us verify the correctness of quantum code and quantum circuit behavior without needing to access the actual quantum hardware. In simulation each shot, which collapse at the state indicating eigenvalue of the observer, represents a single realization of the quantum system [18,19].

When working with quantum circuit notation, it is important to carefully consider Qiskit’s operational sequence against the standard matrix multiplication order. In quantum computing, operations on qubits are represented by unitary matrices, and the order in which these matrices are multiplied can significantly affect the final state of the qubits. Therefore, understanding the convention Qiskit uses for this multiplication is essential for constructing and simulating quantum circuits correctly. This ensures that our theoretical calculations and the actual implementation in Qiskit are consistent. The following subsection delves into the specifics of this order convention. It explains how Qiskit sequences operations and how this affects the resulting state of the quantum system.

3.1. Convention in Qiskit to Find Significant Qubits

We use two circuits, as shown in Figure 2, to illustrate the order of operations for both sequential application and tensor products of logic gates in Qiskit.

Figure 2a shows a quantum circuit containing two unitary logic gates: a Hadamard gate and an X (NOT) gate. Knowing the matrix representations of both gates, we will implement the circuit in Qiskit to obtain the final state and simultaneously compute the product of the gate matrices manually. By comparing these two approaches, we can gain Qiskit insight into how Qiskit orders gate operation.

$$\mathrm{H}\mathrm{X}\left|0\right.⟩=\left({\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\right)\left[\begin{array}{c}1\\ 0\end{array}\right]={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{c}1\\ -1\end{array}\right]={\displaystyle \frac{1}{\sqrt{2}}}\left(\left|0\right.⟩-\left|1\right.⟩\right)$$

(8.1)

$$\mathrm{X}\mathrm{H}\left|0\right.⟩=\left(\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]{\displaystyle \frac{1}{\sqrt{2}}}\right)\left[\begin{array}{c}1\\ 0\end{array}\right]={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{c}1\\ 1\end{array}\right]={\displaystyle \frac{1}{\sqrt{2}}}\left(\left|0\right.⟩+\left|1\right.⟩\right)$$

(8.2)

This analysis reveals that Qiskit implementation of sequential gate operations aligns with a right-to-left order of matrix multiplication. It means that the rightmost gate in the circuit is applied first.

The same procedure is applied to analyze the tensor product of quantum gates. We then compare Qiskit’s implementation with manual calculations for tensor products:

$$\left(\mathrm{H}\otimes \mathrm{X}\right)\left|00\right.⟩=\left({\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]\otimes \left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\right)\left|00\right.⟩={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cccc}0& 1& 0& 1\\ 1& 0& 1& 0\\ 0& 1& 0& -1\\ 1& 0& -1& 0\end{array}\right]\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right]={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{c}0\\ 1\\ 0\\ 1\end{array}\right]={\displaystyle \frac{1}{\sqrt{2}}}\left(\left|01\right.⟩+\left|11\right.⟩\right)$$

(9.1)

$$\left(\mathrm{X}\otimes \mathrm{H}\right)\left|00\right.⟩=\left(\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\otimes {\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]\right)\left|00\right.⟩={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cccc}0& 0& 1& 1\\ 0& 0& 1& -1\\ 1& 1& 0& 0\\ 1& -1& 0& 0\end{array}\right]\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right]={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{c}0\\ 0\\ 1\\ 1\end{array}\right]={\displaystyle \frac{1}{\sqrt{2}}}\left(\left|10\right.⟩+|11⟩\right)$$

(9.2)

Here, we see that Qiskit’s implementation of the tensor product coincides with a down-to-up order of matrix multiplication. And therefore, when calculating the effect of tensor products of multiple gates, Qiskit applies the bottom-most gate first, followed by the gates above it in sequence.

3.2. Phase Estimation for Unitary Matrix $\left(2^n\times 2^n\right)$

In this section, we present an example of quantum phase estimation for a $4\times 4$ unitary operator. Our implementation demonstrates the scalability of the QPE algorithm.

$$U=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& e^{j{\displaystyle \frac{2\pi }{5}}}& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& e^{-j{\displaystyle \frac{2\pi }{5}}}\end{array}\right]$$

(10)

Since the given matrix is diagonal, we can show that the true values of the phases associated with a given unitary operator are ${\theta }^{\left(1\right)}=0,{\theta }^{\left(2\right)}=0.2,{\theta }^{\left(3\right)}=0.5$ and ${\theta }^{\left(4\right)}=0.8$. Figure 3 depicts the execution of the QPE algorithm in Qiskit for a unitary operator, using $6$ ancilla qubits and $1024$ total shots (or runs). The histogram in Figure 3a shows the probability distribution of the estimated first phase, ${\theta }^{\left(1\right)}$. The most probable state is $|000000⟩$, with a probability of 1 (or $1024$ shots out of $1024$ total shots). Figure 3b illustrates the histogram plot of the second phase. The most probable state is state $|001101⟩$, with $900$ shots out of $1024$ shots. The two histogram plots left, Figure 3c,d are estimating third and fourth phase with most probable state $|100000⟩$ and $|110011⟩$, in order. It should be noted that the first two $"00"$ digits in the states labels on x-axis refer to the two qubits of the unitary operator’s eigenvector.

The next step in the QPE process is to interpret the obtained states and extract the corresponding phases. This requires post-processing using classical computing techniques. The results of the post-processing step are as follows:

$${\theta }^{\left(1\right)}=\left|000000\right.⟩\stackrel{BinaryFraction}{=}0.000000\stackrel{Decimal}{=}0$$

(11.1)

$${\theta }^{\left(2\right)}=\left|001101\right.⟩\stackrel{BinaryFraction}{=}0.001101\stackrel{Decimal}{=}{\displaystyle \frac{0}{2}}+{\displaystyle \frac{0}{2^2}}+{\displaystyle \frac{1}{2^3}}+{\displaystyle \frac{1}{2^4}}+{\displaystyle \frac{0}{2^5}}+{\displaystyle \frac{1}{2^6}}=0.203125$$

(11.2)

$${\theta }^{\left(3\right)}=\left|100000\right.⟩\stackrel{BinaryFraction}{=}0.100000\stackrel{Decimal}{=}{\displaystyle \frac{1}{2}}+{\displaystyle \frac{0}{2^2}}+{\displaystyle \frac{0}{2^3}}+{\displaystyle \frac{0}{2^4}}+{\displaystyle \frac{0}{2^5}}+{\displaystyle \frac{0}{2^6}}=0.5$$

(11.3)

$${\theta }^{\left(4\right)}=\left|110011\right.⟩\stackrel{BinaryFraction}{=}0.110011\stackrel{Decimal}{=}{\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{2^2}}+{\displaystyle \frac{0}{2^3}}+{\displaystyle \frac{0}{2^4}}+{\displaystyle \frac{1}{2^5}}+{\displaystyle \frac{1}{2^6}}=0.796875$$

(11.4)

In this step, the quantum circuit computed phases are compared to the previously provided true values to evaluate the accuracy of the QPE algorithm.

For completeness, the next section provides the expressions for obtaining the true values of phase using classical computing.

4. Classical Approach to Obtain Phase

Let's consider a unitary operator $U$ and its corresponding eigenvector $v$. We can represent them as follows:

$$U=\left[\begin{array}{ccc}{u}_{11}& \dots & u_{1N}\\ ⋮& ⋮& ⋮\\ u_{N1}& \dots & u_{NN}\end{array}\right]=\left[\begin{array}{c}{u}^{\left(1\right)}\\ ⋮\\ u^{\left(N\right)}\end{array}\right],v=\left[\begin{array}{c}{v}_1\\ ⋮\\ v_N\end{array}\right]$$

(12)

where $U$ is an $N\times N$ matrix, and each $u_{ij}$ represents an element of the matrix $U$. $u^{\left(k\right)}$ indicates $k^{th}$row of $U$. The eigenvalue problem can be expressed as

$$U\left|v\right.⟩=e^{j2\pi \theta }\left|v\right.⟩\to \left[\begin{array}{c}{u}^{\left(1\right)}\\ ⋮\\ u^{\left(N\right)}\end{array}\right]\left[\begin{array}{c}{v}_1\\ ⋮\\ v_N\end{array}\right]=e^{j2\pi \theta }\left[\begin{array}{c}{v}_1\\ ⋮\\ v_N\end{array}\right]$$

(13)

Then, pick out any arbitrary row $u^{\left(k\right)}$. Now, we can express the eigenvalue equation as:

$$e^{j2\pi \theta }=\sum _{i=1}^N{u}^{\left(k\right)}\left({\displaystyle \frac{v_i}{v_k}}\right)$$

(14)

This expression allows us to calculate the phase θ using any row of the unitary operator $U$ and the corresponding components of its eigenvector $v$. Further, it provides a method to verify the results obtained from the quantum phase estimation algorithm using classical computation.

5. Time Evolution of Schrödinger’s Equation with Self-Adjoint Hamiltonian

For a closed system, we can derive a discrete-time version of the Schrödinger’s equation, which is particularly useful for numerical simulations. For a closed system, discrete version of Schrödinger’s equation (considering ℏ=1) can be expressed as

$$\left|\mathsf{\Psi }\left(t_{k+1}\right)\right.⟩=e^{-jH\mathsf{\Delta }t}\left|\mathsf{\Psi }\left(t_k\right)\right.⟩=e^{-jH{n}_t\mathsf{\Delta }t}\left|\mathsf{\Psi }\left(t_0\right)\right.⟩$$

(15)

where an integer number $n_t$ is reserved for time step counter, representing the number of time steps that have passed since the initial state. As long as, Hamiltonian is Self-Adjoint, the time evolution operator $e^{-jH{n}_t\mathsf{\Delta }t}$ is guaranteed to remain unitary $\left(e^{-jH{n}_t\mathsf{\Delta }t}{\left(e^{-jH{n}_t\mathsf{\Delta }t}\right)}^{\mathrm{*}}=e^{-jH{n}_t\mathsf{\Delta }t+jH{n}_t\mathsf{\Delta }t}=I\right)$. In the following, we provide a numerical example of time evolution of Schrödinger’s equation.

Consider a self-adjoint Hamiltonian $H=\left[\begin{array}{cc}5& 2\\ 2& -6\end{array}\right]$, and a time step $\mathsf{\Delta }t=0.1$. By generating the unitary operator as $e^{-jH{n}_t\mathsf{\Delta }t}$ at different time steps, first phases of the operator would be figure below.

Figure 4a illustrates the estimated phase of the Schrödinger operator over a series of time steps. It compares the phase evolution as calculated using both quantum and classical approaches. The x-axis represents the time steps, while the y-axis shows the phase values. This visualization helps in understanding how the phase evolves over time under different computational methods. Figure 4b depicts the difference between the phase estimates obtained through the quantum approach and those obtained through the classical approach. The x-axis represents the time steps, and the y-axis shows the difference in phase values. This plot is helpful in assessing the accuracy and reliability of the quantum method in comparison to the classical method and in highlighting any significant deviations between the two approaches.

6. Schrödinger’s Equation with Non-Self-Adjoint Hamiltonian

What if Hamiltonian is not self-adjoint? Can we still obtain its eigenvalue using quantum computation? In this section, we want to propose an efficient method to overcome this challenge. The decomposition of a non-self-adjoint matrix into the sum of two self-adjoint matrices can be expressed as

$$H=\left({\displaystyle \frac{H+H^{*}}{2}}\right)+j\left({\displaystyle \frac{H-H^{*}}{2j}}\right)=P+jQ\ni P=P^{*},Q=Q^{*}$$

(16)

Consequently, the exponential of the matrix can be written as:

$$e^{-jH\mathsf{\Delta }t}=e^{-j(P+jQ)\mathsf{\Delta }t}=e^{-jP\mathsf{\Delta }t+Q\mathsf{\Delta }t}$$

(17)

Horn et al. [20] demonstrated that for small time steps $\left(\mathsf{\Delta }t\to 0\right)$, we can approximate:

$$e^{-jP\mathsf{\Delta }t+Q\mathsf{\Delta }t}\approx e^{-jP\mathsf{\Delta }t}{e}^{Q\mathsf{\Delta }t}\mathrm{a}\mathrm{s}\mathsf{\Delta }t\to 0$$

(18)

Since $P$ is self-adjoint, $e^{-jP\mathsf{\Delta }t}$is a unitary operator and can be implemented into Quantum circuit for phase estimation. The problematic part of this decomposition lies in $e^{Q\mathsf{\Delta }t}$. This term is not unitary and therefore complicates its direct implementation in quantum circuit. Small tricks are key to overcome this challenge. Since diagonalization of both $e^{Q\mathsf{\Delta }t}$ and $e^{jQ\mathsf{\Delta }t}$ yields similar structure, exponential of eigenvalue of $Q$ multiply $\mathsf{\Delta }t,e^{{\lambda }_q\mathsf{\Delta }t}$versus $e^{j{\lambda }_q\mathsf{\Delta }t}$. Given this similarity, we can implement $e^{jQ\mathsf{\Delta }t}$ in another quantum circuit separately to eventually obtain the phase of matrix $Q$.

The last step is to correlate eigenvalues of the original non-self-adjoint matrix $A$ with decomposed matrices when the original matrix is normal ($A^{\mathrm{*}}A=A{A}^{\mathrm{*}}$). If $A$ is normal, it is diagonalizable with a set of orthonormal eigenvectors [20]. We can substitute matrix $A$ with matrices $P$and $Q$ (Eq. (19.1)). It implies that $U^{\mathrm{*}}PU$ and $U^{\mathrm{*}}QU$are both diagonal as expressed in Eq$.\left(19.2\right)$. Since $P$ and $Q$ are self-adjoint, their eigenvalues would be real and we name them as ${\lambda }_p$ and ${\lambda }_q$, eigenvalues of $P$ and $Q$ respectively. Comparing Eq.$\left(19.1\right)$ and Eq$.\left(19.2\right)$ results in Eq$.\left(19.3\right)$.

$$U^{*}AU=\mathsf{\Lambda }\to U^{*}\left(P+jQ\right)U=\mathsf{\Lambda }$$

(19.1)

$$U^{*}PU+j{U}^{*}QU=diag\left[{\lambda }_p+j{\lambda }_q\right]$$

(19.2)

$$\mathsf{\Lambda }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[\lambda \right]=diag\left[{\lambda }_p+j{\lambda }_q\right]$$

(19.3)

$U$ is reserved to a set of orthonormal eigenvectors of matrix $A$. $\mathsf{\Lambda }$ indicates a diagonal matrix filled with eigenvalues of matrix $A$.

In contrast to the HHL algorithm, which requires converting a non-self-adjoint matrix into a self-adjoint matrix through stacking matrix and its conjugate transpose leading to doubling the dimension, this method doesn’t need any changes in dimension.

7. Conclusions and Future Work

In this study, we demonstrated the scalability of the QPE algorithm through the phase estimation of a 4×4 unitary operator. We also illustrated the time evolution of Schrödinger’s equation and proposed an efficient method for phase estimation of non-self-adjoint Hamiltonians without altering the dimension of the original matrix. Our future work will focus on exploring open systems and their environmental interactions as described by the Lindblad equation [14].