<figure><img src="image_1.png" /></figure>

HAL
open science

# A Practical Guide to Quantum Linear Algebra

Alexander Weiss, Paul Bruzzi, Fabian M Faulstich

► To cite this version:

Alexander Weiss, Paul Bruzzi, Fabian M Faulstich. A Practical Guide to Quantum Linear Algebra. 2025. [hal-05090131](https://hal.science/hal-05090131v1)

HAL Id: hal-05090131

https://hal.science/hal-05090131v1

Preprint submitted on 29 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Practical Guide to Quantum Linear Algebra

Alexander Weiss*† and Paul Bruzzi*‡

Project advisor: Fabian M. Faulstich§

**Abstract.** We present a comprehensive and pedagogical study of two core quantum algorithms in quantum linear algebra: the Hadamard test (HT) and quantum phase estimation (QPE). Building on Prof. Lin's lecture notes on Quantum Algorithms for Scientific Computing [1], we detail their mathematical formulations, circuit implementations, and convergence behavior. Both algorithms are benchmarked on simulators and IBM's 127-qubit Eagle processor (ibm_rensselaer) using a simple eigenvalue problem. While QPE performs well in simulation, its hardware performance degrades beyond $t = 3$ ancilla qubits due to circuit depth and two-qubit gate noise. In contrast, HT shows stable, sampling-based convergence and greater resilience to hardware noise. We further apply QPE to the transverse-field Ising model (TFIM), demonstrating the importance of ansatz quality for accurate ground state energy estimation. These results provide a realistic baseline for implementing phase estimation on noisy intermediate-scale quantum (NISQ) era devices.

**1. Introduction.** Quantum linear algebra lies at the heart of many quantum algorithms that promise computational advantages over classical methods [2–4]. As linear algebra underpins a wide range of scientific computing tasks – from solving systems of equations to eigenvalue problems and signal processing – quantum linear algebra is expected to play a foundational role in quantum-enhanced numerical simulation, machine learning, and optimization (see, e.g., Refs. [5–7]).

While the theoretical development of quantum algorithms has seen rapid progress, their practical realization is limited by the current state of quantum hardware. Today's devices are constrained by short coherence times, gate infidelities, and limited qubit connectivity – collectively referred to as noise – which significantly restrict the depth and complexity of executable circuits [8]. As a result, many of the most promising quantum algorithms, while asymptotically efficient, are not yet viable for deployment on existing quantum processors. Bridging this gap between theoretical algorithm design and practical implementation is essential to realizing quantum utility in the near term.

This manuscript pursues two primary objectives: First, we aim to provide a pedagogical introduction to quantum linear algebra, accessible at the advanced undergraduate level. We focus on two fundamental algorithms, namely, the Hadamard test (HT) and quantum phase estimation (QPE). These serve as the conceptual and technical foundations for many advanced quantum algorithms, including those for Hamiltonian simulation, quantum singular value transformation, and solving linear systems of equations. Second, we provide a realistic comparison between the theoretical predictions and the performance of these algorithms on actual quantum hardware. Specifically, we benchmark both HT and QPE using (i) fault-tolerant simulations and (ii) experiments on an IBM's 127-qubit Eagle-class quantum

<footnote>*These authors contributed equally to this work.</footnote>

<footnote>†Department of Physics, Rensselaer Polytechnic Institute, Troy, NY (alexxweiss34@gmail.com)</footnote>

<footnote>‡Department of Computer Science, Rensselaer Polytechnic Institute, Troy, NY (paulbruzzi@gmail.com)</footnote>

<footnote>§Department of Mathematics, Rensselaer Polytechnic Institute, Troy, NY</footnote>
processor (ibm_rensselaer) [9]. This comparison exposes the strengths and limitations of both methods in the current noisy intermediate-scale quantum (NISQ) era and highlights the critical impact of circuit depth, gate noise, and state preparation on algorithmic performance.

This manuscript was inspired by and closely follows the structure of Ref. [1], Prof. Lin Lin's lecture notes on Quantum Algorithms for Scientific Computing. We supplement his exposition with additional clarifications and visual aids to support undergraduate students from STEM backgrounds who are new to quantum computing. In addition, we provide detailed Jupyter notebook implementations of all major concepts and circuits, which are available on our project GitHub repository, see Ref. [10].

The manuscript is organized as follows. Section 2 introduces the foundational elements of quantum computation, including qubits, single-qubit gates, and quantum phase. Section 3 presents the Hadamard test as a simple method for estimating expectation values. Section 4 discusses the discrete Fourier transform and its quantum implementation via the quantum Fourier transform. Section 5 builds on this to construct the quantum phase estimation algorithm. In section 6, we compare the performance of HT and QPE in terms of accuracy and runtime on both simulated and real hardware. Section 7 applies QPE to the transverse-field Ising model (TFIM), illustrating the role of ansatz quality. Conclusions are summarized in section 8.

We gratefully acknowledge the use of IBM Quantum services for this work. The views expressed here are those of the authors and do not represent the official views of IBM or the IBM Quantum team.

## 2. Computation on a Quantum Computer.

### 2.1. Qubits.
Quantum bits, or qubits, are the fundamental units of quantum information. Analogous to classical bits, which can be in one of two states (0 or 1), qubits can exist in two computational basis states denoted as $|0\rangle$ and $|1\rangle$ in Dirac notation, often referred to as spin up and spin down states respectively as in Refs. [3,4].

**Definition 2.1.** The computational basis states of a single qubit are defined as

$$ (2.1) \qquad |0\rangle := \begin{bmatrix} 1 \\ 0 \end{bmatrix} \quad \text{and} \quad |1\rangle := \begin{bmatrix} 0 \\ 1 \end{bmatrix}. $$

A general qubit state is a normalized linear combination (superposition) of these basis states:

$$ (2.2) \qquad |\psi\rangle = c_1|0\rangle + c_2|1\rangle, \quad \text{where } c_1, c_2 \in \mathbb{C}, |c_1|^2 + |c_2|^2 = 1. $$

**Remark 2.2.** Unlike classical bits, qubits may exist in superpositions. For instance,

$$ (2.3) \qquad |+\rangle := \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle) = \frac{1}{\sqrt{2}}\begin{bmatrix} 1 \\ 1 \end{bmatrix} \quad \text{and} \quad |-\rangle := \frac{1}{\sqrt{2}}(|0\rangle - |1\rangle) = \frac{1}{\sqrt{2}}\begin{bmatrix} 1 \\ -1 \end{bmatrix}, $$

are both valid qubit states. See Appendix A.1 for visualizations.

For systems involving multiple qubits, the full state is represented by a tensor product of individual qubit states. This is commonly abbreviated using the following notation:

$$ (2.4) \qquad |\psi_1\rangle \otimes |\psi_2\rangle = |\psi_1\rangle|\psi_2\rangle = |\psi_1\psi_2\rangle. $$
**Example 2.3.** A system with two qubits in states $|0\rangle$ and $|1\rangle$ is represented as

$$ (2.5) \qquad |0\rangle \otimes |1\rangle = \begin{bmatrix} 1 \\ 0 \end{bmatrix} \otimes \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \end{bmatrix} = |01\rangle. $$

**2.2. Measurement and Phase.** Measuring a qubit collapses its superposition into one of the computational basis states, with probabilities determined by the squared modulus of its inner product (overlap) with the measured state.

**Definition 2.4.** The probability of measuring a state $|\phi\rangle$ and finding it in state $|\psi\rangle$ is

$$ (2.6) \qquad \mathbb{P}(X \in \{|\psi\rangle\}) = |\langle\psi|\phi\rangle|^2, $$

where X is a random variable representing the measurement of $|\phi\rangle$.

**Example 2.5.** The probability of measuring $|+\rangle$ and obtaining $|0\rangle$ is

$$ (2.7) \qquad |\langle 0 | + \rangle|^2 = \left| \langle 0 | \left( \frac{1}{\sqrt{2}} (|0\rangle + |1\rangle) \right) \rangle \right|^2 = \left| \frac{1}{\sqrt{2}} (\langle 0 | 0 \rangle + \langle 0 | 1 \rangle) \right|^2 = \left| \frac{1}{\sqrt{2}} \right|^2 = \frac{1}{2}, $$

where we used the orthonormality of the computational basis, i.e., $\langle 0|0\rangle = 1$ and $\langle 0|1\rangle = 0$. In vector notation this reads

$$ (2.8) \qquad |\langle 0 | + \rangle|^2 = \left| \begin{bmatrix} 1 \\ 0 \end{bmatrix}^\dagger \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 1 \end{bmatrix} \right|^2 = \left| \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 1 \end{bmatrix} \right|^2 = \left| \frac{1}{\sqrt{2}} \right|^2 = \frac{1}{2}, $$

where “†” denotes the adjoint. This means that a qubit initially in the $|+\rangle$ state will collapse to the $|0\rangle$ state upon measurement with a probability of 1/2.

**Remark 2.6.** Since the modulus square is invariant under rotation, there exist multiple states that collapse to the same state with the same probability, i.e, when multiplying the expansion coefficients of a state by $e^{i2\pi\varphi}$ with $\varphi \in [0, 1)$. The quantity $e^{i2\pi\varphi}$ is called *phase* and $2\pi\varphi$ is called *phase angle* of the state. Here, $\varphi$ acts as the fraction of a full rotation of $2\pi$ around the complex unit circle. There are two types of phase, *global phase* which is the phase between different state vectors and *relative phase* which is the phase between bit strings in the same state vector. The idea of relative and global phase is expanded upon in **Example 2.9**.

**2.3. Single-Qubit Gates.** Quantum gates manipulate the state of qubits. The most common single-qubit gates include the Pauli-X (bit flip), Pauli-Z (phase flip), $R_z$ (phase rotation), and the Hadamard gate (superposition).

**Definition 2.7.** The matrix representations of the Pauli-X, Pauli-Z, $R_z$, and Hadamard gate in the computational basis are:

$$ (2.9) \quad X := \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}, \quad Z := \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}, \quad R_z(\theta) = \begin{bmatrix} e^{-i\theta/2} & 0 \\ 0 & e^{i\theta/2} \end{bmatrix} \quad \text{and} \quad H := \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}. $$
The following table summarizes their action on the computational basis states:

<table><caption>**Table 1**
The result of the X-gate, Z-gate, R<sub>z</sub>(θ)-gate, and H-gate acting on the |0⟩ and |1⟩ states.</caption><thead><tr><th>Input</th><th>X</th><th>Z</th><th>R<sub>z</sub>(θ)</th><th>H</th></tr></thead><tbody><tr><td>|0⟩</td><td>|1⟩</td><td>|0⟩</td><td>e<sup>-iθ/2</sup>|0⟩</td><td><sup>1</sup>⁄<sub>√2</sub>(|0⟩ + |1⟩)</td></tr><tr><td>|1⟩</td><td>|0⟩</td><td>−|1⟩</td><td>e<sup>iθ/2</sup>|1⟩</td><td><sup>1</sup>⁄<sub>√2</sub>(|0⟩ − |1⟩)</td></tr></tbody></table>

**Remark 2.8.** The effect of applying a phase gate (e.g., Z-gate, R<sub>z</sub>(θ)-gate) does not change the measurement outcome, but changes the quantum state by a complex phase. These phase factors can be interpreted as eigenvalues:

$$ (2.10) \qquad U|\psi\rangle = \lambda|\psi\rangle, \quad \text{with} \quad \lambda = \langle\psi|U|\psi\rangle. $$

Since quantum states are normalized, |λ| = 1, implying λ = e<sup>i2πϕ</sup> (c.f. **Remark 2.6**).

**Example 2.9.** Expanding on further on **Remark 2.6**, the action of R<sub>z</sub>(θ) on |+⟩ yields:

$$ (2.11) \qquad R_z(\theta)|+\rangle = \frac{1}{\sqrt{2}} \left( e^{-i\theta/2}|0\rangle + e^{i\theta/2}|1\rangle \right) = \frac{e^{-i\theta/2}}{\sqrt{2}} \left( |0\rangle + e^{i\theta}|1\rangle \right). $$

revealing a global phase of e<sup>-iθ/2</sup> and a relative phase of e<sup>iθ</sup> between the basis components.
See **Appendix B** for more information about single-qubit gates.

**2.4. Circuit Diagrams.** A circuit diagram (e.g. Figure 1) is a graphical representation that describes the workflow of a quantum algorithm. Quantum circuits are read from left to right as operations or gates are applied over time on lines representing the qubits. In general, one line in a circuit diagram corresponds to a single qubit; however, if a line has a slash through it, e.g. the third line in Figure 1, it represents multiple qubits. Applying a gate to a qubit is represented as a box on the qubit’s line. Some gate operations depend on one or more control qubit states. These controls are denoted by a point drawn over the controlling qubit(s) with a line connected to the gate it is controlling; see |π<sub>2</sub>⟩ in Figure 1. To denote the swapping of two qubit states, a (×) is drawn over the swapping qubits with a line connecting them; see |π<sub>3</sub>⟩ in Figure 1. Applying a gate to multiple qubits is shown with a box over all the qubits such as at |π<sub>4</sub>⟩ in Figure 1. Measurement is depicted with a “meter-like” symbol in a box on the qubit being measured, as shown in the first line at the end of the circuit. Finally, a double line denotes classical information, such as after a measurement has been performed.

**3. The Hadamard Test.** Although quantum measurements discard global and relative phase information (see **Remark 2.6**), many quantum algorithms critically rely on access to this phase. The Hadamard Test (HT) is a simple yet powerful circuit primitive for estimating the expectation value ⟨ψ|U|ψ⟩ of a unitary operator U with respect to a quantum state |ψ⟩. Though elementary in structure, the Hadamard Test forms the conceptual backbone of numerous state-of-the-art quantum algorithms for linear algebra, including those for matrix functions, Hamiltonian simulation, and singular value transformation. Because the expectation value is generally complex, its real and imaginary parts must be estimated separately
<figure><img src="image_3.png" /><figcaption>**Figure 1.** An example quantum circuit diagram depicting a multi-qubit slash on the third line in the state $|\psi\rangle$, a Hadamard ($H$) gate on the first qubit (line), a controlled $X$-gate controlled on the first qubit targeting the second, a swap operation between the first and second qubits, an arbitrary multi-qubit $U$ gate, and a measurement operator on the first qubit.</figcaption></figure>

using two circuit variants: the real Hadamard test and the imaginary Hadamard test as in Refs. [11, 12].

### 3.1. The HT Circuit.

**Definition 3.1.** Let $U \in \mathbb{C}^{n \times n}$ be unitary and define the phase gate $S = \sqrt{Z}$. The real and imaginary Hadamard tests (HT) are defined as the circuits depicted below in Figure 2.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** (Left) The real Hadamard test circuit and (Right) the imaginary Hadamard test circuit. $H$ is the Hadamard gate, $S$ is the phase gate, $U$ is a unitary gate. The topmost qubit starting in the state $|0\rangle$ is the control qubit, and the lower qubit(s) starting in the state $|\psi\rangle$ are the computation bits representing the eigenstate.</figcaption></figure>

The measurement of the control qubit in the real HT is represented by a random variable $X_{\text{Re}}$ which is permitted to be either $|0\rangle$ or $|1\rangle$. Instead of referencing this random variable directly, we use the shorthand notation:

$$ (3.1) \qquad \mathcal{P}_{\text{Re}|0\rangle} := \mathbb{P}(X_{\text{Re}} \in \{|0\rangle\}) \quad \text{and} \quad \mathcal{P}_{\text{Re}|1\rangle} := \mathbb{P}(X_{\text{Re}} \in \{|1\rangle\}). $$

Similarly, a random variable $X_{\text{Im}}$ being $|0\rangle$ or $|1\rangle$ describes the measurement of the control qubit in the imaginary HT, where we use the shorthand notation:

$$ (3.2) \qquad \mathcal{P}_{\text{Im}|0\rangle} := \mathbb{P}(X_{\text{Im}} \in \{|0\rangle\}) \quad \text{and} \quad \mathcal{P}_{\text{Im}|1\rangle} := \mathbb{P}(X_{\text{Im}} \in \{|1\rangle\}). $$

**Proposition 3.2.** The measurement probabilities for the real Hadamard test are

$$ (3.3) \qquad \mathcal{P}_{\text{Re}|0\rangle} = \frac{1}{2}(1 + \text{Re}\langle\psi|U|\psi\rangle) \quad \text{and} \quad \mathcal{P}_{\text{Re}|1\rangle} = \frac{1}{2}(1 - \text{Re}\langle\psi|U|\psi\rangle). $$

*Proof.* Consider the real HT (left panel in Figure 2 in Definition 3.1) with initial state

$$ (3.4) \qquad |\pi_0\rangle = |0\rangle \otimes |\psi\rangle = |0\rangle|\psi\rangle. $$
A Hadamard gate is then applied to the control qubit, which yields the $|+\rangle$ state, i.e.,

$$ (3.5) \qquad |\pi_1\rangle = H|0\rangle I|\psi\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)|\psi\rangle = |+\rangle|\psi\rangle. $$

Next, a controlled unitary gate applies the U operation to $|\psi\rangle$ if the control qubit is $|1\rangle$, i.e.,

$$ (3.6) \qquad |\pi_2\rangle = (|0\rangle\langle 0| \otimes I + |1\rangle\langle 1| \otimes U) |\pi_1\rangle = \frac{1}{\sqrt{2}}|0\rangle I|\psi\rangle + \frac{1}{\sqrt{2}}|1\rangle U|\psi\rangle. $$

The last operation is a Hadamard gate, yielding the final state

$$ (3.7) \qquad |\pi_3\rangle = |0\rangle \frac{I+U}{2} |\psi\rangle + |1\rangle \frac{I-U}{2} |\psi\rangle, $$

where the probability of measuring the control qubit in $|0\rangle$ is

$$ (3.8) \qquad \mathcal{P}_{\text{Re}|0\rangle} = |\langle 0|\pi_3\rangle|^2 = \left| \frac{I+U}{2} |\psi\rangle \right|^2 = \frac{2 + \langle \psi|(U+U^\dagger)|\psi\rangle}{4} = \frac{1}{2}(1 + \text{Re}\langle \psi|U|\psi\rangle). $$

Similarly, the probability of measuring the control qubit in $|1\rangle$ is

$$ (3.9) \qquad \mathcal{P}_{\text{Re}|1\rangle} = |\langle 1|\pi_3\rangle|^2 = \left| \frac{I-U}{2} |\psi\rangle \right|^2 = \frac{2 - \langle \psi|(U+U^\dagger)|\psi\rangle}{4} = \frac{1}{2}(1 - \text{Re}\langle \psi|U|\psi\rangle). \quad \blacksquare $$

Having established the measurement probabilities for the real Hadamard test, we can characterize the real part of the expectation value.

**Corollary 3.3.** *The real part of the expectation value is given by*

$$ (3.10) \qquad \text{Re}\langle\psi|U|\psi\rangle = \mathcal{P}_{\text{Re}|0\rangle} - \mathcal{P}_{\text{Re}|1\rangle} = 2\mathcal{P}_{\text{Re}|0\rangle} - 1 = 1 - 2\mathcal{P}_{\text{Re}|1\rangle}. $$

*Proof.* Note that Proposition 3.2 yields

$$ (3.11) \qquad 2\mathcal{P}_{\text{Re}|0\rangle} - 1 = (1 + \text{Re}\langle\psi|U|\psi\rangle) - 1 = \text{Re}\langle\psi|U|\psi\rangle $$

and

$$ (3.12) \qquad \text{Re}\langle\psi|U|\psi\rangle = 1 - (1 - \text{Re}\langle\psi|U|\psi\rangle) = 1 - 2\mathcal{P}_{\text{Re}|1\rangle}. $$

Combining equations (3.11) and (3.12) yields,

$$ (3.13) \qquad \text{Re}\langle\psi|U|\psi\rangle = \frac{1}{2}(1 + \text{Re}\langle\psi|U|\psi\rangle) - \frac{1}{2}(1 - \text{Re}\langle\psi|U|\psi\rangle) = \mathcal{P}_{\text{Re}|0\rangle} - \mathcal{P}_{\text{Re}|1\rangle}. $$

For a more detailed step-by-step proof, we refer the reader to Appendix C.2. ■

The analysis for the imaginary Hadamard test follows analogously. We here report the measurement probabilities together with the resulting expectation value, and refer the reader to Appendix C.3 for the detailed proofs.

**Proposition 3.4.** *The measurement probabilities for the imaginary Hadamard test are*

$$ (3.14) \qquad \mathcal{P}_{\text{Im}|0\rangle} = \frac{1}{2}(1 + \text{Im}\langle\psi|U|\psi\rangle) \quad \text{and} \quad \mathcal{P}_{\text{Im}|1\rangle} = \frac{1}{2}(1 - \text{Im}\langle\psi|U|\psi\rangle). $$

**Corollary 3.5.** *The imaginary part of the expectation value is given by*

$$ (3.15) \qquad \text{Im}\langle\psi|U|\psi\rangle = \mathcal{P}_{\text{Im}|0\rangle} - \mathcal{P}_{\text{Im}|1\rangle} = 2\mathcal{P}_{\text{Im}|0\rangle} - 1 = 1 - 2\mathcal{P}_{\text{Im}|1\rangle}. $$
**3.2. Estimating Eigenphases.** Returning to the task of phase estimation, recall that the eigenvalues of $U$ are of the form

$$ (3.16) \qquad \lambda = e^{i\theta} = \cos\theta + i\sin\theta, $$

and thus $\mathrm{Re}\langle\psi|U|\psi\rangle = \cos\theta$. So the real probabilities can be written as

$$ (3.17) \qquad \mathcal{P}_{\mathrm{Re}|0\rangle} = \frac{1}{2}(1 + \cos\theta) \quad \text{and} \quad \mathcal{P}_{\mathrm{Re}|1\rangle} = \frac{1}{2}(1 - \cos\theta), $$

hence,

$$ (3.18) \qquad \theta = \arccos(2\mathcal{P}_{\mathrm{Re}|0\rangle} - 1) = \arccos(1 - 2\mathcal{P}_{\mathrm{Re}|1\rangle}). $$

Similarly, $\mathrm{Im}\langle\psi|U|\psi\rangle = \sin\theta$. Thus, $\theta$ can also be estimated by

$$ (3.19) \qquad \theta = \arcsin(2\mathcal{P}_{\mathrm{Im}|0\rangle} - 1) = \arcsin(1 - 2\mathcal{P}_{\mathrm{Im}|1\rangle}). $$

Finally, $\theta$ can be estimated using both the real and imaginary parts of the expectation value. Each part can be thought of as a vector on the complex unit circle with the real part on the x-axis and the imaginary part on the y-axis. Taking the arccos of the imaginary part or the arccos of the real part is effectively calculating the vector of the entire complex expectation value. This means that the arctan of the imaginary part over the real part also gives $\theta$. Therefore, the angle $\theta$ can be estimated with the following three methods:

$$ (3.20) \qquad \theta = \arccos(\mathrm{Re}\langle\psi|U|\psi\rangle) = \arcsin(\mathrm{Im}\langle\psi|U|\psi\rangle) = \arctan\left(\frac{\mathrm{Im}\langle\psi|U|\psi\rangle}{\mathrm{Re}\langle\psi|U|\psi\rangle}\right). $$

*Remark 3.6.* Note that using both the real and imaginary parts of the expectation value to estimate $\theta$ with arctan increases the error of the estimation, since it uses data from two random variables rather than just one.

**3.3. HT Phase Estimation Example.** To demonstrate the Hadamard test in practice, consider the unitary matrix

$$ (3.21) \qquad U = \begin{bmatrix} 1 & 0 \\ 0 & e^{i\theta} \end{bmatrix}, \quad \text{with} \quad \theta = \frac{1}{2} + \frac{1}{2^4} = 0.5625. $$

We will use this simple matrix as a running example throughout this article. Our numerical tests are performed using Qiskit 1.3.1 alongside AerSimulator [13] for the simulated quantum data. We estimate $\theta$ by applying both the real and imaginary Hadamard tests and measuring 100 shots, see Figure 3.

Theoretical expectation values (c.f. (3.17)) follow from:

$$ (3.22) \qquad \mathcal{P}_{\mathrm{Re}|0\rangle} = \frac{1}{2}(1 + \cos(0.5625)) \approx 0.9230 \quad \text{and} \quad \mathcal{P}_{\mathrm{Re}|1\rangle} = \frac{1}{2}(1 - \cos(0.5625)) \approx 0.0770 $$

whereas the measured values from a simulation with 100 shots yield

$$ (3.23) \qquad \mathcal{P}_{\mathrm{Re}|0\rangle} = 0.90 \quad \text{and} \quad \mathcal{P}_{\mathrm{Re}|1\rangle} = 0.10. $$
<figure><img src="image_5.png" /><figcaption>Figure 3. Counts for the real (left) and imaginary (right) Hadamard tests for 100 shots with a phase angle $\theta = 0.5625$.</figcaption></figure>

Similarly, for the imaginary circuit, the computed probabilities are

$$ (3.24) \quad \mathcal{P}_{\text{Im}|0\rangle} = \frac{1}{2}(1 + \sin(0.5625)) \approx 0.7667 \quad \text{and} \quad \mathcal{P}_{\text{Im}|1\rangle} = \frac{1}{2}(1 - \sin(0.5625)) \approx 0.2333, $$

whereas the measured values from a simulation with 100 shots yield

$$ (3.25) \quad \mathcal{P}_{\text{Im}|0\rangle} = 0.75 \quad \text{and} \quad \mathcal{P}_{\text{Im}|1\rangle} = 0.25. $$

This illustrates the discrepancy between the theoretical and experimental values using 100 shots for both HTs. As expected, both probabilities always add to 1 (c.f. Appendix C.4). Provided these probabilities, the real and imaginary expectation values can be found using (3.10) and (3.15), respectively. This yields

$$ (3.26) \quad \text{Re}\langle\psi|U|\psi\rangle = 0.80 \quad \text{and} \quad \text{Im}\langle\psi|U|\psi\rangle = 0.50, $$

whereas the theoretical expectation values obtained from (3.22) and (3.24) are

$$ (3.27) \quad \text{Re}\langle\psi|U|\psi\rangle = 0.8460 \quad \text{and} \quad \text{Im}\langle\psi|U|\psi\rangle = 0.5334. $$

Thus, the experimental and theoretical expectation values are approximately what we expect. Finally, either or both of these values can be used to estimate the phase angle $\theta$ using any equation in (3.20).

*Remark 3.7.* Since trigonometric functions are cyclic, inverse trigonometric functions have multiple possible values of $\theta$ that they can return. The standard for these types of function is to return the value closest to zero, which is not necessarily the true phase. Therefore, the outputs of both arccos and arcsin are used to determine which quadrant the phase is in. The function arccos ranges from 0 to $\pi$, and arcsin from $-\pi/2$ to $\pi/2$. arccos is used to
determine whether the angle is in the positive x or negative x part of the graph, and arcsin is used to determine if it is in the positive y or negative y part of the graph. Using the angle returned from arccos, if arcsin gives a negative value, then the angle must be found in the negative y part of the graph. To do this, the estimated angle can be set to $2\pi - \theta$. Using the angle returned from arcsin, if arccos gives a value greater than $\pi/2$, then the angle must be determined in the negative x part of the graph. To do this, the estimated angle is set to $\pi - \theta$. For a graphical explanation we refer the reader to Ref. [14].

We emphasize that the above experiment simulated 100 shots, leading to an estimated phase angle $\theta_{est} = 0.6435$ (c.f. $\theta = 0.5625$) when using cos to determine the angle and sin to determine the sign. Note that the accuracy of the estimation depends directly on the number of shots. In the following section, we will connect the convergence behavior to Monte-Carlo type estimators, providing a rate of convergence for the HTs.

**3.4. Error Analysis.** In quantum computing, repeated measurements—known as shots—are used to estimate the probability of a quantum system collapsing into a given state. This is inherently a statistical process: each measurement samples from an underlying probability distribution that is not directly accessible. To infer this distribution, we employ tools from statistical estimation, most notably the empirical distribution function (EDF), which approximates the true distribution based on a finite number of observations. This forms the foundation for quantifying the accuracy and convergence of algorithms like the HTs.

Each single-qubit measurement of the control qubit yields a Bernoulli random variable $X \in \{0, 1\}$, with success probability $p = P(X = 0)$. For $N_s$ independent shots, the EDF is defined as

$$ (3.28) \qquad \hat{F}(x) = \frac{1}{N_s} \sum_{i=1}^{N_s} \mathbb{1}_{\{X_i \le x\}}, $$

where $X_i \sim X$ are i.i.d. and are corresponding to each measurement result. Here, $\mathbb{1}_{\{X_i \le x\}}$ is the indicator function which returns 1 when $X_i \le x$ and 0 otherwise.

It is often necessary to estimate the probability of a specific outcome. This is a special case of the EDF where it simplifies to the empirical probability estimator which focuses on the proportion of trials yielding a specific outcome. For example, the probability of measuring 0 on a quantum computer with $N_s$ shots is estimated by

$$ (3.29) \qquad \hat{p}_0 = \frac{1}{N_s} \sum_{i=1}^{N_s} \mathbb{1}_{\{0\}}(X_i). $$

This is also called a Monte Carlo estimator which uses repeated sampling to predict the probability of a variety of outcomes. Recall that for the HTs, only a single qubit is measured, and thus the $X_i$'s are Bernoulli distributed. Tying single-qubit measurement to common probability theory, the weak law of large numbers (see, e.g., Ref. [15, Section 10.2]) yields that $\hat{p}_0$ in (3.29) converges to $p_0 = E(X)$ in probability.
**Proposition 3.8.** The average of $N_s$ single-qubit measurements converges in probability. Moreover, the convergence rate of the empirical estimator $\hat{p}_0$ can be calculated from the variance

$$ (3.30) \qquad \mathbb{V}(\hat{p}_0) = \frac{p_0(1-p_0)}{N_s}, $$

which yields the standard deviation

$$ (3.31) \qquad \sigma = \sqrt{\frac{p_0(1-p_0)}{N_s}}. $$

This translates to taking $N_s$ single-qubit measurements as follows.

**Proposition 3.9.** Let $\epsilon > 0$. The number of shots required to recover $p_0$ with precision $\epsilon$ via $\hat{p}_0$ scales as

$$ (3.32) \qquad \epsilon \ge \sigma = \sqrt{\frac{p_0(1-p_0)}{N_s}} \Rightarrow \epsilon^2 \ge \frac{p_0(1-p_0)}{N_s} \Rightarrow N_s \ge \frac{p_0(1-p_0)}{\epsilon^2} \in \Omega\left(\frac{1}{\epsilon^2}\right). $$

The variance of the control qubit's output is dependent on $p_0$ and thus $\sigma$ in (3.31). Plotting $\sigma$ against $N_s$ in Figure 4 for $10^3-10^5$ shots reveals that the error approaches zero as the number of shots increases when experimentally tested. Moreover, the estimated angle approaches $\theta = 0.5625$ as the number of shots increases, see the table in Figure 4.

<figure><img src="image_6.png" /><figcaption>Error w.r.t. Number of Shots</figcaption></figure>

<table><thead><tr><td>Shot Count</td><td>Estim. Angle θ</td><td>Error</td></tr></thead><tbody><tr><td>100</td><td>0.4027</td><td>0.1598</td></tr><tr><td>500</td><td>0.4864</td><td>0.0761</td></tr><tr><td>1000</td><td>0.5394</td><td>0.0231</td></tr><tr><td>5000</td><td>0.5830</td><td>0.0205</td></tr><tr><td>10000</td><td>0.5698</td><td>0.0073</td></tr><tr><td>50000</td><td>0.5569</td><td>0.0056</td></tr><tr><td>100000</td><td>0.5620</td><td>0.0005</td></tr></tbody></table>

**Figure 4.** (Left) A log-log scaled plot of the standard deviation of the Hadamard test $\sigma$ against the number of shots $N_s$ converging like $N_s^{-1/2}$. (Right) A table showing the estimated angle $\theta$ converging to the expected value $\theta = 0.5625$ and its error approaching 0 as the number of shots increases.

**4. Quantum Fourier Transform.** Fourier analysis plays a central role in scientific computing, offering a way to transform problems from the time or spatial domain to the frequency domain, where many tasks – such as filtering, compression, or solving differential equations –
become more tractable [16]. In quantum computing, a related transformation known as the quantum Fourier transform (QFT) serves as a key subroutine in algorithms including Shor's algorithm, quantum phase estimation, and quantum signal processing (see, e.g., Refs. [17–19]). This section introduces both the classical Discrete Fourier Transform (DFT) and its quantum counterpart, the QFT, and outlines their mathematical foundations and circuit implementations.

**4.1. Discrete Fourier Transform.** We begin by recalling the classical Discrete Fourier Transform, which maps a complex-valued vector in the time domain into the frequency domain.

**Definition 4.1.** The (forward) discrete Fourier transform (DFT) is defined as the map,

$$ (4.1) \qquad \mathcal{F}: \mathbb{C}^N \to \mathbb{C}^N; \quad x = \begin{bmatrix} x_0 \\ \vdots \\ x_{N-1} \end{bmatrix} \mapsto \tilde{x} = \begin{bmatrix} \tilde{x}_0 \\ \vdots \\ \tilde{x}_{N-1} \end{bmatrix}, $$

where

$$ (4.2) \qquad \tilde{x}_k = \frac{1}{\sqrt{N}} \sum_{j=0}^{N-1} x_j e^{2\pi i \frac{jk}{N}}. $$

**Example 4.2.** Consider the continuous signal

$$ (4.3) \qquad x(t) = \sin(2\pi t) + \sin(2\pi(11.2)t) + \sin(2\pi(18.4)t) $$

sampled at a rate of 100Hz, i.e., 100 samples per second. Taking 1 second of equally spaced samples yields the plot in the left panel of Figure 5. Applying the DFT in (4.2) where $x_j$ are the samples, reveals the frequency amplitudes $\tilde{x}_k$ in the right panel of Figure 5.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** (Left) A sample of $x(t)$ in (4.3), and (Right) its Fourier transform showing pronounced peaks at 1 Hz, 11.2 Hz and 18.4 Hz.</figcaption></figure>
Conversely, given a set of frequencies, the inverse discrete Fourier transform reconstructs the corresponding time-domain signal from its frequency-domain representation.

**Definition 4.3.** The inverse discrete Fourier transform (IDFT) is defined as the map,

$$ (4.4) \qquad \mathcal{F}^{-1} : \mathbb{C}^N \to \mathbb{C}^N ; \quad \tilde{x} = \begin{bmatrix} \tilde{x}_0 \\ \vdots \\ \tilde{x}_{N-1} \end{bmatrix} \mapsto x = \begin{bmatrix} x_0 \\ \vdots \\ x_{N-1} \end{bmatrix}, $$

where

$$ (4.5) \qquad x_j = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} \tilde{x}_k e^{-2\pi i \frac{jk}{N}}. $$

**4.2. QFT Algorithm.** The quantum Fourier transform performs an operation similar to the DFT, but acts on the amplitudes of a quantum state. It is defined over computational basis states, which represent integers in binary.

**Definition 4.4.** Let $j \in \mathbb{N}$ and $n \in \mathbb{N}$ be the number of qubits. The (forward) quantum Fourier transform (QFT) is defined as the map

$$ (4.6) \qquad |j\rangle \mapsto \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} e^{2\pi i \frac{jk}{2^n}} |k\rangle. $$

Here we express $j$ and $k$ as binary numbers to represent the computational basis state vectors. See Appendix A.2 for more details.

**Remark 4.5.** The QFT can be expressed using the unitary matrix

$$ (4.7) \qquad \text{QFT}_n = \frac{1}{\sqrt{2^n}} \sum_{j=0}^{2^n-1} \sum_{k=0}^{2^n-1} e^{2\pi i \frac{jk}{2^n}} |k\rangle\langle j| = \frac{1}{\sqrt{2^n}} \begin{bmatrix} 1 & 1 & 1 & \dots & 1 \\ 1 & \omega & \omega^2 & \dots & \omega^{2^n-1} \\ 1 & \omega^2 & \omega^4 & \dots & \omega^{2(2^n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{2^n-1} & \omega^{2(2^n-1)} & \dots & \omega^{(2^n-1)^2} \end{bmatrix}, $$

where $\omega = e^{2\pi i \frac{1}{2^n}}$.

**Example 4.6.** Taking $n = 1$ in (4.7), the matrix form of $\text{QFT}_1$ is

$$ (4.8) \qquad \text{QFT}_1 = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & e^{\frac{2\pi i}{2}} \end{bmatrix} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & \cos(\pi) + i \sin(\pi) \end{bmatrix} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} = H. $$

Therefore, the QFT for a single qubit is simply the Hadamard gate, i.e., this gate transforms the computational basis states $|0\rangle$ and $|1\rangle$ into the Fourier basis states $|+\rangle$ and $|-\rangle$, respectively.

**Proposition 4.7.** QFT$_n$ acting on a quantum state $|j\rangle$ yields

$$ (4.9) \qquad \text{QFT}_n |j\rangle = \frac{1}{\sqrt{2^n}} \bigotimes_{l=1}^{n} (|0\rangle + e^{2\pi i \frac{j}{2l}} |1\rangle). $$
**Proof.** The map in Definition 4.4 can be rewritten using fractional binary notation where $k$ and $j$ are binary numbers, i.e.,

$$ (4.10) \qquad k = k_1 k_2 \dots k_n 0, \quad j = j_1 j_2 \dots j_n 0. $$

Therefore, QFT$_n$ in (4.6) acting on a state $|j\rangle$ can be written as

$$ (4.11) \qquad \text{QFT}_n |j\rangle = \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} e^{2\pi i \frac{jk}{2^n}} |k\rangle = \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} \prod_{l=1}^{n} e^{2\pi i k_l \frac{j}{2^l}} |k\rangle, $$

where

$$ (4.12) \qquad \frac{jk}{2^n} = k_1 \frac{j}{2^1} + k_2 \frac{j}{2^2} + \dots + k_n \frac{j}{2^n} = \sum_{l=1}^{n} k_l \frac{j}{2^l}. $$

Since $k$ is a binary number, each $k_l$ can only be 0 or 1. Thus the sum over $k$ becomes the direct product of the sum of each individual $k_l$ from 0 to 1, i.e.,

$$ (4.13) \qquad \text{QFT}_n |j\rangle = \frac{1}{\sqrt{2^n}} \sum_{k_1=0}^{1} e^{2\pi i k_1 \frac{j}{2^1}} |k_1\rangle \otimes \sum_{k_2=0}^{1} e^{2\pi i k_2 \frac{j}{2^2}} |k_2\rangle \otimes \dots \otimes \sum_{k_n=0}^{1} e^{2\pi i k_n \frac{j}{2^n}} |k_n\rangle. $$

This can be rewritten as multiple tensor products

$$ (4.14) \qquad \text{QFT}_n |j\rangle = \frac{1}{\sqrt{2^n}} \bigotimes_{l=1}^{n} \left( |0\rangle + e^{2\pi i \frac{j}{2^l}} |1\rangle \right). \qquad \blacksquare $$

**4.3. QFT Circuit Implementation.** The QFT algorithm has an efficient implementation as a quantum circuit consisting of Hadamard gates and controlled phase rotations.

**Proposition 4.8.** Let $R_n = R_z(\frac{\pi}{2^{n-1}})$. The QFT circuit is depicted in Figure 6.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** The QFT circuit where $R_n = R_z(\frac{\pi}{2^{n-1}})$ (Definition 2.7). Note that $R_1 = H$.</figcaption></figure>

**Proof.** If the circuit in Figure 6 performs a Fourier transform, then the final state of the circuit is the expression in (4.14), therefore a circuit without the swap gates – which we will
call the swap-less QFT circuit (QFT) – will swap the qubit order of the output. This implies that the final state of the swap-less QFT circuit would be

$$ (4.15) \qquad \overleftarrow{\text{QFT}}_n |j\rangle = \frac{1}{\sqrt{2^n}} \bigotimes_{l=n}^{1} \left( |0\rangle + e^{2\pi i \frac{j}{2^l}} |1\rangle \right). $$

Using the swap-less QFT circuit and (4.15), we construct an inductive proof verifying the circuit in Figure 6 does indeed perform a Fourier transform. Consider the swap-less QFT circuit with a single qubit. Here, the circuit becomes a single $R_1$ gate which is equivalent to a Hadamard gate. According to (4.15), the initial state of the qubit transforms to,

$$ (4.16) \qquad \overleftarrow{\text{QFT}}_1 |j\rangle = \frac{1}{\sqrt{2^1}} \bigotimes_{l=1}^{1} \left( |0\rangle + e^{2\pi i \frac{j}{2^l}} |1\rangle \right) = \frac{1}{\sqrt{2}} (|0\rangle + e^{2\pi i 0 \cdot j_1} |1\rangle) = \begin{cases} |+\rangle, & \text{if } j=0 \\ |-\rangle, & \text{if } j=1. \end{cases} $$

Since this is exactly the effect of a Hadamard gate, when $n = 1$, the swap-less QFT circuit correctly performs the map in (4.15), showing the induction base.

Since $\overleftarrow{\text{QFT}}_n$ represents the swap-less QFT circuit for n qubits, it then follows from (4.15) that $\overleftarrow{\text{QFT}}_{n+1}$ can be written as

$$ (4.17) \qquad \overleftarrow{\text{QFT}}_{n+1} |j\rangle = \frac{1}{\sqrt{2}} (|0\rangle + e^{2\pi i 0 \cdot j_1 \dots j_{n+1}} |1\rangle) \otimes \overleftarrow{\text{QFT}}_n |j_2 \dots j_{n+1}\rangle. $$

From the inductive hypothesis, the $\overleftarrow{\text{QFT}}_n |j_2 \dots j_{n+1}\rangle$ part of $\overleftarrow{\text{QFT}}_{n+1}$ is correctly performing the map in (4.15). Thus, if the first bit is being transformed correctly, the entire mapping is correct.

The phase on the $|1\rangle$ state in (4.17) can be rewritten as a combination of rotations specific to each binary digit of j, i.e.,

$$ (4.18) \qquad e^{2\pi i 0 \cdot j_1 \dots j_{n+1}} = e^{2\pi i \frac{j_1}{2^1}} e^{2\pi i \frac{j_2}{2^2}} \dots e^{2\pi i \frac{j_{n+1}}{2^{n+1}}} = e^{j_1(\pi i)} e^{j_2(\frac{\pi i}{2})} \dots e^{j_{n+1}(\frac{\pi i}{2^n})}. $$

Since the phase rotation can be rewritten as a product of rotations that are dependent only on a single qubit, each rotation can be directly converted to a single controlled rotation gate. The rotation for each of these gates is defined as

$$ (4.19) \qquad r_n = \frac{\pi}{2^{n-1}}, $$

where $n$ represents the qubit's number/label. Upon inspecting $\overleftarrow{\text{QFT}}_{n+1}$, the gates defined here and in Proposition 4.8 are equivalent. This means that if the swap-less QFT circuit with $n$ qubits performs the map $\overleftarrow{\text{QFT}}_n$, the swap-less QFT circuit with $n+1$ qubits performs the map $\overleftarrow{\text{QFT}}_{n+1}$. Therefore, by induction the swap-less QFT circuit performs the map $\overleftarrow{\text{QFT}}_n$ for all $n \ge 1$. Since the swap-less QFT performs the mapping $\overleftarrow{\text{QFT}}_n$, the circuit in Figure 6 must perform the QFT mapping in Definition 4.4. ■
We refer the interested reader to Appendix D for a more detailed explanation on this circuit’s construction.

**Definition 4.9.** Let $j \in \mathbb{N}$. The inverse quantum Fourier transform (IQFT) is defined as the map

$$ (4.20) \qquad |k\rangle \mapsto \frac{1}{\sqrt{2^n}} \sum_{j=0}^{2^n-1} e^{-2\pi i \frac{kj}{2^n}} |j\rangle. $$

**Proposition 4.10.** Let $R_n^\dagger = R_z(\frac{-\pi}{2^{n-1}})$. The IQFT circuit is given in Figure 7.

<figure><img src="image_9.png" /><figcaption>Figure 7. The IQFT circuit where $R_n^\dagger = R_z(\frac{-\pi}{2^{n-1}})$. Note that $R_1^\dagger = H$.</figcaption></figure>

*Proof.* Since QFT can be represented as a unitary matrix, its inverse is the conjugate transpose of QFT. Therefore, the IQFT circuit is constructed by applying the dagger operation to the QFT circuit in Figure 7.

**5. Quantum Phase Estimation.** The Quantum Phase Estimation (QPE) algorithm is one of the most fundamental procedures in quantum computing. Like the Hadamard test, it estimates the phase $\theta$ in an eigenvalue of a unitary operator $U$, i.e., when

$$ (5.1) \qquad U|\psi\rangle = e^{2\pi i\theta}|\psi\rangle $$

where $|\psi\rangle$ is a known eigenstate. QPE generalizes the Hadamard test by using a register of ancilla qubits to extract multiple bits of $\theta$, thereby increasing the resolution of the estimate. It forms the foundation of many quantum algorithms, including those for factoring, linear systems, and Hamiltonian simulation (see, e.g. Refs. [1, 3, 18, 20]).

**5.1. The QPE Algorithm.**

**Definition 5.1.** Let $U \in \mathbb{C}^{n \times n}$ be unitary. The quantum phase estimation algorithm is defined via the circuit in Figure 8.

**Proposition 5.2.** The final state of the circuit in Figure 8 is

$$ (5.2) \qquad |\psi_3\rangle = \frac{1}{2^n} \sum_{x=0}^{2^n-1} \sum_{k=0}^{2^n-1} e^{2\pi i k \left(\varphi - \frac{x}{2^n}\right)} |x\rangle|\psi\rangle. $$
<figure><img src="image_10.png" /><figcaption>Figure 8. Quantum phase estimation circuit where H is the Hadamard gate, $U \in C^{n \times n}$ is unitary, and QFT$^\dagger$ is the inverse QFT. The top qubits are the ancilla qubits and the bottom qubits represent the state $|\psi\rangle$.</figcaption></figure>

*Proof.* The initial state of the circuit in Figure 8 can be represented as

$$ (5.3) \qquad |\psi_0\rangle = |0\rangle^{\otimes n} |\psi\rangle. $$

A Hadamard gate is then applied to all $n$ ancillas generating a uniform superposition, i.e.,

$$ (5.4) \qquad |\psi_1\rangle = \frac{1}{\sqrt{2^n}} (|0\rangle + |1\rangle)^{\otimes n} |\psi\rangle = |+\rangle^{\otimes n} |\psi\rangle. $$

Next, the controlled $U$ gates are applied. Since $U|\psi\rangle = e^{2\pi i\varphi}|\psi\rangle$, this yields

$$ (5.5) \qquad U^{2^j} |\psi\rangle = (e^{2\pi i\varphi})^{2^j} |\psi\rangle = e^{2\pi i 2^j \varphi} |\psi\rangle. $$

Therefore, by controlling this gate and applying it to all $n$ qubits, it only affects the state for the case where it is $|1\rangle$. Applying the controlled $U^{2^j}$, with $0 \le j \le n-1$, yields

$$ (5.6) \qquad |\psi_2\rangle = \frac{1}{\sqrt{2^n}} (|0\rangle + e^{2\pi i \varphi 2^{n-1}} |1\rangle) \otimes \dots \otimes (|0\rangle + e^{2\pi i \varphi 2^0} |1\rangle) \otimes |\psi\rangle = \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} e^{2\pi i \varphi k} |k\rangle |\psi\rangle, $$

where $k$ is the integer representation of an $n$-bit binary number. The IQFT is then applied resulting in the state

$$ (5.7) \qquad \frac{1}{\sqrt{2^n}} \sum_{k=0}^{2^n-1} e^{2\pi i \varphi k} |k\rangle |\psi\rangle \xrightarrow{\text{QFT}_N^\dagger} \frac{1}{2^n} \sum_{x=0}^{2^n-1} \sum_{k=0}^{2^n-1} e^{2\pi i k (\varphi - \frac{x}{2^n})} |x\rangle |\psi\rangle = |\psi_3\rangle. \qquad \blacksquare $$

*Remark 5.3.* QPE with one ancilla qubit yields the same circuit as the real HT.

*Example 5.4.* To test QPE, 4 ancilla qubits with $10^3$ shots were used on a fault-tolerant Qiskit AerSimulator [13], estimating the phase angles $\varphi = 1/4$ and $\varphi = 2/3$ for the unitary matrix

$$ (5.8) \qquad U = \begin{bmatrix} 1 & 0 \\ 0 & e^{i2\pi\varphi} \end{bmatrix}. $$
<figure><img src="image_11.png" /><figcaption>Figure 9. QPE applied to a state with a phase angle of $\varphi = 1/4$ (left) and $\varphi = 2/3$ (right).</figcaption></figure>

After a sampling, the probability of measuring each specific state was calculated by dividing the number of times each state was counted by the total number of shots.

Since the phase angle $\varphi = \frac{1}{4}$ can be exactly represented with 2 binary digits, using 4 ancilla qubits in QPE returns the phase angle deterministically, see left panel in Figure 9. However, since the phase angle $\varphi = \frac{2}{3}$ cannot be exactly represented with any number of binary digits, the outcome is probabilistic, and thus multiple measurements are required, see right panel in Figure 9.

**5.2. Error Analysis.** QPE yields a distribution over estimates centered around the true phase. However, when the targeted phase has no exact binary bit expansion, QPE returns the nearest binary approximation with high probability. Quantifying this error is the subject of this section. Throughout this section, we make the assumption that

$$ (5.9) \qquad U|\psi_j\rangle = e^{i2\pi\varphi_j}|\psi_j\rangle, $$

where $0 \le \varphi_0 \le \varphi_1 \cdots \le \varphi_{N-1} \le 1$, and that $\varphi_0$ does not have an exact $d$-bit representation. QPE is moreover applied to the initial state $|0^t\rangle|\psi_0\rangle$ where $t > d$.

**Proposition 5.5.** In order to obtain the phase $\varphi_0$ with an accuracy of $\epsilon = 2^{-d}$ and a success probability of at least $1 - \delta$, a total of $t = d + \log_2(\delta^{-1})$ ancilla qubits are necessary to store the phase value. Furthermore, because $T\epsilon = \delta^{-1}$, the simulation run time must be $T = (\epsilon\delta)^{-1}$.
*Proof.* Performing QPE on the initial state $|0\rangle^{\otimes t}|\psi_0\rangle$ gives

$$ (5.10) \qquad \begin{aligned} |0\rangle^{\otimes t}|\psi_0\rangle \xrightarrow{H^{\otimes t} \otimes I} \frac{1}{\sqrt{2^t}} \sum_{j \in [2^t]} |j\rangle|\psi_0\rangle \\ \xrightarrow{U} \frac{1}{\sqrt{2^t}} \sum_{j \in [2^t]} |j\rangle e^{i2\pi j \varphi_0} |\psi_0\rangle \\ \xrightarrow{\text{QFT}^\dagger} \sum_{k' \in [2^t]} \left( \frac{1}{2^t} \sum_{j \in [2^t]} e^{i2\pi j \left(\varphi_0 - \frac{k'}{2^t}\right)} \right) |k'\rangle|\psi_0\rangle. \end{aligned} $$

For simplicity, it can be said that

$$ (5.11) \qquad |0\rangle^{\otimes t}|\psi_0\rangle \xrightarrow{\text{QPE}} \sum_{k'} \gamma_{0,k'} |k'\rangle|\psi_0\rangle, \quad \text{with} \quad \gamma_{0,k'} = \frac{1}{2^t} \sum_{j \in [2^t]} e^{i2\pi j \left(\varphi_0 - \frac{k'}{2^t}\right)}. $$

Recognizing that $\gamma_{0,k'}$ is a geometric series with the common ratio $r = e^{i2\pi(\varphi_0 - \frac{k'}{2^t})}$ and the coefficient $a = \frac{1}{2^t}$, we find that

$$ (5.12) \qquad \gamma_{0,k'} = \frac{1}{2^t} \frac{1-r^{2^t}}{1-r} = \frac{1}{2^t} \frac{1-e^{i2\pi(\varphi_0 - \tilde{\varphi}_{k'})}}{1-e^{i2\pi(\varphi_0 - \tilde{\varphi}_{k'})}}, \quad \text{with} \quad \tilde{\varphi}_{k'} = \frac{k'}{2^t}, $$

Since $e^{i2\pi x}$ is a 1-periodic function, $x \bmod 1$ is the only measurable value. Therefore, when estimating the phase $k'_0$, we seek the precision $\epsilon = \frac{1}{2^d} = \frac{2^{t-d}}{2^t}$ such that

$$ (5.13) \qquad \left|\varphi_0 - \tilde{\varphi}_{k'_0}\right| < \epsilon. $$

Since for any $\theta \in [-\pi, \pi]$

$$ (5.14) \qquad \left|1 - e^{i\theta}\right| = \sqrt{(1 - \cos\theta)^2 + \sin^2\theta} = \sqrt{2(1 - \cos\theta)} = 2 \left|\sin\frac{\theta}{2}\right| \ge \frac{2}{\pi} |\theta|, $$

and

$$ (5.15) \qquad \left|1 - e^{i\theta}\right| \le \left|1\right| + \left|e^{i\theta}\right| = \sqrt{1^2} + \sqrt{\cos^2\theta - (i \sin\theta)^2} = 1 + \sqrt{\cos^2\theta + \sin^2\theta} = 1 + 1 = 2, $$

we find that

$$ (5.16) \qquad \left|\gamma_{0,k'}\right| \le \frac{2}{2^t 2\pi \frac{2}{\pi} |\varphi_0 - \tilde{\varphi}_{k'}|_1} = \frac{1}{2^{t+1} |\varphi_0 - \tilde{\varphi}_{k'}|_1}. $$

The binary measurement outcome $k'_0$ can be viewed as a random variable where the probability of obtaining $\tilde{\varphi}_{k'_0}$ at least a distance $\epsilon$ away from the expected value $\varphi_0$ is

$$ (5.17) \qquad \mathbb{P}\left(\left|\varphi_0 - \tilde{\varphi}_{k'_0}\right| \ge \epsilon\right) = \sum_{|\varphi_0 - \tilde{\varphi}_{k'}| \ge \epsilon} |\gamma_{0,k'}|^2 \le \left(\frac{1}{2^{t+1} |\varphi_0 - \tilde{\varphi}_{k'}|_1}\right)^2 = \frac{1}{4T^2 |\varphi_0 - \tilde{\varphi}_{k'}|_1^2}, $$
where $T = 2^t$. In other words, this expression represents the probability that the estimated phase is not within the desired error. Setting $x = |\varphi_0 - \tilde{\varphi}_{k'}|_1$, this yields

$$ (5.18) \qquad \mathbb{P}(|x| \ge \epsilon) = \sum_{|x| \ge \epsilon} \frac{1}{4T^2 x^2} =: \sum_{|x| \ge \epsilon} f(x). $$

Since $f$ is a symmetric function,

$$ (5.19) \qquad \sum_{|x| \ge \epsilon} \frac{1}{4T^2 x^2} = \sum_{x \ge \epsilon} \frac{2}{4T^2 x^2} = \sum_{x \ge \epsilon} \frac{1}{2T^2 x^2}. $$

As exemplified in Example 5.7, we see the sum in (5.19) is upper bounded by the integral from $\epsilon$ onwards. Noting that

$$ (5.20) \qquad \sum_{x \ge \epsilon} \frac{1}{x^2} \le \int_{\epsilon}^{\infty} \frac{1}{x^2} dx = \frac{1}{\epsilon}, $$

the expression in (5.19) can be bounded, i.e., $T \ge 1$ yields

$$ (5.21) \qquad \sum_{|x| \ge \epsilon} \frac{1}{4T^2 x^2} \le \frac{1}{2T^2 \epsilon} \le \frac{1}{2T \epsilon}. $$

Since $T\epsilon = 2^{t-d}$, we set $t-d = \log_2(\delta^{-1})$ so that $T\epsilon = \delta^{-1}$. For $0 < \delta < 1$, we then find

$$ (5.22) \qquad \mathbb{P}\left(\left|\varphi_0 - \tilde{\varphi}_{k'_0}\right| \ge \epsilon\right) \le \frac{1}{2T\epsilon} = \frac{\delta}{2} \le \delta. $$

Therefore, in order to obtain the phase $\varphi_0$ to an accuracy of $\epsilon = 2^{-d}$ with a success probability of at least $1 - \delta$, a total of $t = d + \log_2(\delta^{-1})$ ancilla qubits are necessary to store the phase value. Furthermore, since $T\epsilon = \delta^{-1}$, the simulation run time must be $T = (\epsilon\delta)^{-1}$. ■

**Example 5.6.** Set $\varphi_0 = 0.35$ (not exactly representable in binary). For $t=6$ and $t=10$, the QPE distribution $|\gamma_{0,k'}|$ in (5.16) becomes increasingly peaked around the closest binary approximation to $\varphi_0$, see Figure 10. As $t \to \infty$, this distribution converges to a Dirac delta. In fact, for the case where $\varphi_0$ has an exact $d$-bit representation (meaning that $\varphi_0 = \tilde{\varphi}_{k'_0}$ for some $k'_0$), then $\gamma_{0,k'} = \delta_{k',k'_0}$ (the Dirac delta function). Inspecting the function in (5.7), we see that it peaks at $x = 2^n\varphi$, which is the input state.

**Example 5.7.** Let $\epsilon = 2^{-5}$, in Figure 11 we visualize $f(x)$ defined in (5.18).

**6. Hadamard Test vs. Quantum Phase Estimation.** The Hadamard Test (HT) and Quantum Phase Estimation (QPE) both serve to estimate the phase $\varphi$ in an eigenvalue of the form $e^{i2\pi\varphi}$ of a unitary operator $U$. In this section, we compare the two algorithms with respect to accuracy and runtime, using both classical simulations and real quantum hardware (IBM Eagle One system at RPI, known as ibm_rensselaer).

To ensure a fair comparison, both algorithms are used to estimate the same eigenphase $\varphi = 1/3$ of the unitary matrix

$$ (6.1) \qquad U = \begin{bmatrix} 1 & 0 \\ 0 & e^{i2\pi\varphi} \end{bmatrix}, $$

to precision $\epsilon$.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** $|\gamma_{0,k'}|$ plotted against $\tilde{\phi}_{k'}$ for $\phi_0 = 0.35$. (Left) With $t = 6$, and (Right) with $t = 10$ ancilla qubits.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** $\sum_{|x| \ge \epsilon} f(x)$ plotted against $x$ for $\epsilon = 2^{-5}$.</figcaption></figure>

**6.1. Accuracy Comparison.** A direct comparison of the error between the Hadamard Test (HT) and Quantum Phase Estimation (QPE) is nontrivial due to the fundamentally different nature of their errors: HT exhibits sampling error, while QPE exhibits approximation error based on binary resolution. To visualize the convergence behavior of each algorithm as a function of its quantum resources, we plot their respective errors on a shared graph. For this comparison, we run each algorithm over 20 independent trials. The HT is tested with shot counts ranging from $10^3$ to $10^5$, while QPE is tested using 1 to 20 ancilla qubits, each run with 1000 shots. All experiments were conducted both on Qiskit's AerSimulator [13] and on ibm_rensselaer, an IBM Quantum System One Eagle processor [13] (see Figure 12). These parameters were selected to span a realistic range of quantum resources given current simulation and hardware constraints.
Up to this point the reported data has been gathered by classical simulation corresponding to a fault-tolerant quantum computer, but does this still hold on a real quantum machine?

As shown in the left panel of Figure 12, the Hadamard test (HT) exhibits excellent agreement between the observed and theoretically predicted error across all tested configurations, even when executed on real quantum hardware. In contrast, Quantum Phase Estimation (QPE) deviates significantly from theoretical predictions when the number of ancilla qubits exceeds $t = 3$. For $t \le 3$, QPE performs as expected, with error decreasing according to the predicted $O(2^{-t})$ scaling. However, for $t > 3$, the error no longer decreases and may even worsen with additional ancilla qubits.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** (Left) A comparison of the approximation error, and (Right) a comparison of the runtime. We present data for the HT and QPE algorithms obtained from a quantum emulator and an IBM Eagle 127-qubit quantum processor.</figcaption></figure>

This degradation is attributed to the growing circuit depth and gate count: Each added ancilla qubit introduces a new layer of controlled $U^{2j}$ operations, many of which are decomposed into multiple two-qubit gates. On current noisy intermediate-scale quantum (NISQ) hardware, this increased gate complexity significantly amplifies noise and decoherence, undermining QPE's theoretical advantages. To reduce the impact of noise, we employed Qiskit's transceiver with various `optimization_level` settings (0-3), which aim to simplify the circuit and reduce the total gate count. Although higher optimization levels (e.g., `optimization_level = 3`) are effective at minimizing noise in many applications, in this experiment they were insufficient to recover QPE performance beyond $t = 3$. This suggests that current hardware remains too noisy to realize the full potential of QPE in regimes that require deeper circuits.

To better understand the performance degradation of QPE on hardware, it is useful to analyze the circuit depth, which serves as a proxy for accumulated gate noise and decoherence effects. While circuit depth captures the overall number of sequential operations, the primary
contributors to hardware error are two-qubit gates, which are more error-prone than single-qubit gates and dominate the noise profile in most NISQ devices [21]. As illustrated in Figure 13, the number of two-qubit gates in QPE circuits grows rapidly with each additional ancilla qubit. Although Qiskit's higher transpiler optimization levels (2 and 3) significantly reduce this gate count compared to lower levels, the improvement is not sufficient beyond $t = 3$. In particular, circuits for $t \ge 4$ remain too deep – despite aggressive optimization – to yield reliable output on current quantum hardware, underscoring the practical limitations of QPE in its standard form on today's machines.

<figure><img src="image_15.png" /><figcaption>**Figure 13.** 2-qubit gate circuit depth plotted against the number of ancilla qubits for QPE. The “noise threshold” represents the shortest 2-qubit circuit depth where the qubit operations create too much noise to correctly perform their operations.</figcaption></figure>

**6.2. Runtime Comparison.** Runtime is measured as the total execution time for a fixed accuracy target. Although HT and QPE have different resource requirements (shots vs. ancilla qubits), they share a dependence on the physical capabilities of the backend.

The right panel in Figure 12 shows a runtime comparison of HT and QPE on a quantum emulator and an IBM Eagle 127-qubit quantum processor. For HT, the simulated experiment requires less runtime per shot compared to its quantum counterpart. Since the total runtime of HT scales directly with the runtime per shot, we observe that the runtime of HT on the quantum computer exhibits an offset, resulting in a runtime consistently higher than that of the quantum emulator. For QPE, the classical simulation runtime is initially lower than the runtime on the quantum computer but increases exponentially. On the quantum computer, the QPE runtime remains close to linear.

Note that for $t = 1$ ancilla qubits, the QPE circuit is structurally equivalent to the real HT circuit. However, in practice, the total runtime of HT is observed to be approximately twice that of QPE at this setting. This discrepancy arises because HT requires two separate circuits – one to estimate the real part and another to estimate the imaginary part of the expectation value – in order to reconstruct the full complex phase. In contrast, QPE encodes the phase directly into a binary register, requiring only one circuit to extract the phase information.
We emphasize that to increase the precision of the estimated phase on the quantum emulator, the runtime of both the HT and QPE will increase. However, on the quantum computer, regardless of the value of t, QPE's runtime remains nearly linear, while HT's runtime exponentially increases with the number of shots.

**Remark 6.1.** Even with QPE's inaccuracy, its runtime is always very low on a quantum machine. Therefore, if the noise can somehow be mitigated better as the ancillary qubits increase, the QPE algorithm will become extremely effective in estimating phases with high binary bit precision. Since QPE is unreliable for binary representations of $d > 3$ bits, unless the phase can be exactly represented with $d \le 3$ bits, the HT will give higher precision. Since HT is always slower than QPE, if the phase can be exactly represented with $d \le 3$ bits, QPE will give higher precision and lower runtime.

**7. Transverse-Field Ising Model.** The transverse-field Ising model (TFIM) is a widely studied quantum spin system that captures fundamental features of quantum phase transitions and many-body quantum dynamics. It consists of a one-dimensional chain of spin-$\frac{1}{2}$ particles interacting with their nearest neighbors along the z-axis and subjected to an external magnetic field in the transverse x-direction.

### 7.1. Hamiltonian and Dynamics.

**Definition 7.1.** Let $n \in \mathbb{N}$ denote the number of particles and $g \in \mathbb{R}_{\ge 0}$ the magnetic field strength. The governing matrix $H \in \mathbb{C}^{2^n \times 2^n}$, i.e., the Hamiltonian, of the one-dimensional transverse-field Ising model with open boundary conditions is defined as

$$ (7.1) \qquad H = - \sum_{i=0}^{n-1} Z_i Z_{i+1} - g \sum_{i=0}^{n} X_i, $$

where $Z_i$ represents the Pauli-Z operator applied to the ith particle and $X_i$ represents the Pauli-X operator applied to the ith particle.

**Remark 7.2.** The physical interpretation of the TFIM in Definition 7.1 is that $n$ describes the number of particles exposed to the transverse external magnetic field. The first term in (7.1) represents the nearest-neighbor coupling, while the second term describes the interaction with the magnetic field, where $g$ quantifies its strength.

The time evolution of a quantum system governed by this Hamiltonian satisfies the time-dependent Schrödinger equation:

$$ (7.2) \qquad i \frac{\partial}{\partial t} |\psi(t)\rangle = H |\psi(t)\rangle, $$

with initial condition $|\psi(t = 0)\rangle = |\psi_0\rangle$. The solution for $|\psi(t)\rangle$ is given by

$$ (7.3) \qquad |\psi(t)\rangle = e^{-itH} |\psi_0\rangle. $$

Hence, simulating time evolution or extracting eigenvalues of $H$ reduces to computing the eigenphases of the unitary operator $U = e^{-itH}$, which can be performed using QPE or related quantum algorithms.
**7.2. Hamiltonian Decomposition.** To implement time evolution on a quantum computer, we must approximate the operator $U = e^{-itH}$. When $H$ is expressed as a sum of non-commuting components, care must be taken in the exponentiation. Suppose

$$ (7.4) \qquad H = A + B, $$

where $A$ and $B$ are unitary operators. Since $e^H \neq e^A e^B$ unless $A$ and $B$ commute (see Ref. [22, Remark 2.4]), one must use the Lie product formula from Ref. [23] to approximate:

$$ (7.5) \qquad U = e^{-itH} = \lim_{L \to \infty} \left( e^{-itA/L} e^{-itB/L} \right)^L. $$

This is known as the Trotter-Suzuki decomposition. For discrete $L$, this is an approximation. However, in the case of the TFIM, it can be exactly simulated with a quantum circuit as in Ref. [24]. Thus, TFIM serves as a valuable testbed for quantum algorithm development, allowing exploration of algorithmic behavior under both commuting and non-commuting dynamics. See Appendix E for a more detailed explanation of trotter error.

**7.3. State Preparation and Initialization.** While QPE can estimate eigenvalues of any unitary operator, its success depends on preparing an initial state with sufficient overlap with the target eigenstate (typically the ground state). This initial guess is known as the ansatz.

To explore this sensitivity, consider the TFIM with $n = 4$ spins and magnetic field strength $g = 10$. Different initial states were prepared by applying a single-qubit $R_y(\theta)$ gate to each qubit, with $\theta$ ranging from $0$ to $2\pi$. This rotation spans a full set of superpositions about the $y$ axis between $|0\rangle$ and $|1\rangle$, thereby modifying the overlap with the true ground state. Figure 14 illustrates that the accuracy of the QPE estimate varies with the ansatz, demonstrating the importance of physically informed initialization strategies.

<figure><img src="image_16.png" /><figcaption>**Figure 14.** The calculated ground state phase of the TFIM with $g = 10$, $n = 4$, and $t = 4$ ancilla qubits, using a varying rotation about the $y$ axis applied to the initial guess.</figcaption></figure>
For QPE to yield reliable estimates, the initial state must have significant overlap with the target eigenstate. In the case of the TFIM, physically motivated initializations can be derived from the structure of the Hamiltonian in (7.1). When $g \ll 1$, the magnetic field term becomes negligible and the coupling term dominates, aligning spins along the z-axis. In this regime, initializing all qubits to $|0\rangle$ results in a large overlap with the ground state.

Conversely, when $g \gg 1$, the magnetic field term dominates, favoring alignment with its eigenvectors—namely, the *x*-basis states $|+\rangle$ and $|-\rangle$. This can be achieved by applying Hadamard gates to each qubit.

For intermediate values of $g$, where neither term is negligible, the ground state structure becomes more elaborate. However, the Hadamard-initialized state remains a reasonable heuristic, as it captures partial alignment along both axes. The quality of any ansatz ultimately depends on how much each term in (7.1) contributes to the system's energy.

**7.4. Simulation of TFIM with QPE.** We exemplify the simulation of the TFIM with QPE using a Qiskit state vector emulator for the system parameters $n = 4$ and $g = 0$, $g = 1$, and $g = 10^8$. The system Hamiltonian was normalized and scaled by a factor of 1/8 and increased by a factor of 1/2, making the Hamiltonian positive definite, to scale the eigenvalues to be in the range 3/8 to 5/8. Following the state preparation outlined in subsection 7.3, we initialize all qubits with $|0\rangle$ when $g = 0$, all qubits with $|+\rangle$ when $g = 1$ and $g = 10^8$. The ground state energy of the following examples is 0.375 which is closest to the 4-bit binary number 0.0110.

<figure><img src="image_17.png" /><figcaption>**Figure 15.** The QPE State vectors of TFIM with (Left) $g = 0$, (Center) $g = 1$, and (Right) $g = 10^8$.</figcaption></figure>

In all subfigures of Figure 15, the computational basis state with the highest measurement probability corresponds to the ground state energy of the TFIM, indicating that QPE successfully identifies the correct eigenvalue. However, in the center panel, corresponding to the intermediate case where $g = 1$, the probability of measuring the correct ground state (0110) is notably lower than in the cases where $g = 0$ and $g \gg 1$. This reduction stems from the more complex spectral structure at intermediate field strengths: the energy levels are more densely packed and the ground state is less isolated. As a result, the initial state used has non-negligible overlap with excited states, leading to a broader distribution of measurement outcomes. This highlights the importance of selecting an ansatz informed by the physical characteristics of the system. In regimes where eigenstates are not well-separated, accurate phase estimation requires a deeper understanding of the Hamiltonian's structure to ensure sufficient overlap with the target eigenstate.
**8. Conclusion.** This work was guided by the lecture notes of Prof. Lin on Quantum Algorithms for Scientific Computing [1], which provided a conceptual foundation for our exploration of quantum phase estimation and its practical implementation. Building on this framework, we presented a pedagogical yet rigorous initial study of two essential quantum algorithms: the *Hadamard Test* (HT) and *Quantum Phase Estimation* (QPE). After introducing the mathematical structure and circuit implementations of both algorithms, we evaluated their performance in realistic computational settings.

A key contribution of this work is a detailed comparison of HT and QPE on both simulated and real quantum hardware, in our case, the IBM Quantum System One at RPI (ibm_rensselaer) [13], a 127-qubit Eagle-class quantum processor. Our main findings are:

* On ideal fault-tolerant simulators, HT and QPE perform in accordance with theoretical expectations: QPE shows exponential accuracy scaling in the number of ancilla qubits, while HT exhibits convergence consistent with Monte Carlo sampling error.

* On real hardware, HT remains robust across a wide range of shot counts, closely matching simulated behavior. QPE, however, deviates significantly from its ideal performance once the number of ancilla qubits exceeds $t = 3$, due to increased circuit depth and accumulation of two-qubit gate errors.

* Analysis of gate depth and transpiler optimization revealed a steep rise in the number of two-qubit gates with increasing QPE register size. Despite leveraging Qiskit's highest transpilation optimization level, QPE circuits with $t \ge 4$ proved too deep to maintain fidelity on Eagle-class hardware.

* In a case study involving the Transverse-Field Ising Model (TFIM), we demonstrated the role of physically informed ansatz selection in enhancing QPE performance, especially in regimes with nontrivial ground state structure.

Together, these results provide a realistic benchmark for basic quantum linear algebra methods on near-term quantum devices. While QPE remains a cornerstone of long-term algorithmic development, our experiments suggest that sampling-based methods like HT currently offer greater resilience to hardware noise, making them more suitable for early applications of quantum computing in scientific computing. As hardware matures and deeper circuits become feasible, these findings will serve as a reference point for re-evaluating algorithmic choices in the quantum simulation of physical systems.

**Acknowledgments.** The authors are thankful for the useful discussions with Dr. Cameron Cogburn, Dr. Osama Raisuddin, and Mr. Ibrahim Abdulazeez.
REFERENCES

[1] Lin Lin. Lecture notes: Quantum algorithms for scientific computing, 2024. Available at https://arxiv.org/abs/2201.08309.

[2] Mauro E S Morales, Lirandê Pira, Philipp Schleich, et al. Quantum linear system solvers: A survey of algorithms and applications. arXiv preprint arXiv:2411.02522v1, 2024.

[3] Michael A. Nielsen and Isaac L. Chuang. *Quantum Computation and Quantum Information*. Cambridge University Press, 10th anniversary edition, 2010.

[4] Thomas G. Wong. *Introduction to Classical and Quantum Computing*. Rooted Grove, Omaha, Nebraska, 2022.

[5] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd. Quantum algorithm for linear systems of equations. *Phys. Rev. Lett.*, 103:150502, 2009.

[6] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe. Quantum singular value transformation and beyond. In *Proc. 51st ACM Symposium on Theory of Computing (STOC)*, pages 193–204, 2019.

[7] Shi Jin and Nana Liu. Quantum simulation of discrete linear dynamical systems and simple iterative methods in linear algebra via schrödingerisation. arXiv preprint arXiv:2304.02865, 2023.

[8] John Preskill. Quantum computing in the NISQ era and beyond. *Quantum*, 2:79, 2018.

[9] IBM Quantum. Ibm unveils 127-qubit *eagle* quantum processor, 2021. IBM Research Blog.

[10] Alex Weiss and Paul Bruzzi. Intro to quantum linear algebra, 2025. https://github.com/Quantum-Linear-Algebra/Intro-to-Quantum-Linear-Algebra.

[11] Richard Cleve, Artur Ekert, Chiara Macchiavello, and Michele Mosca. Quantum algorithms revisited. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 454(1969):339–354, 1998.

[12] Kosuke Mitarai and Keisuke Fujii. Methodology for replacing indirect measurements with direct measurements. *Phys. Rev. Res.*, 1:013006, Aug 2019.

[13] Qiskit Development Team. Qiskit: An open-source framework for quantum computing, 2024.

[14] https://www.desmos.com/calculator/oebpz5tpyt.

[15] Joseph K. Blitzstein and Jessica Hwang. *Introduction to Probability*. CRC Press, Boca Raton, Florida, second edition, 2019.

[16] Ronald N. Bracewell. *The Fourier Transform and its Applications*. McGraw-Hill, 3rd edition, 2000.

[17] Peter W. Shor. Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer. *SIAM Journal on Computing*, 26(5):1484–1509, October 1997.

[18] Ryan LaRose. The phase estimation algorithm. https://www.ryanlarose.com/uploads/1/1/5/8/115879647/quic03-phase-estimation.pdf, n.d. Accessed: 2025-01-12.

[19] Guang Hao Low and Isaac L. Chuang. Optimal hamiltonian simulation by quantum signal processing. *Physical Review Letters*, 118(1):010501, 2017.

[20] various authors. *Qiskit Textbook*. Github, 2023.

[21] Ali Javadi-Abhari, Matthew Treinish, Kevin Krsulich, Christopher J. Wood, Jake Lishman, Julien Gacon, Simon Martiel, Paul D. Nation, Lev S. Bishop, Andrew W. Cross, Blake R. Johnson, and Jay M. Gambetta. Quantum computing with Qiskit, 2024.

[22] Robin Zhang. The baker–campbell–hausdorff formula, 2017. Available at https://math.mit.edu/~robinz/files/BCHFormula.pdf.

[23] Masuo Suzuki. Generalized trotter's formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems. *Communications in Mathematical Physics*, 51:183–190, 1976.

[24] Alba Cervera-Lierta. Exact Ising model simulation on a quantum computer. *Quantum*, 2:114, December 2018.
Appendix A. Qubits.

**A.1. The Bloch Sphere.** A common technique for visualizing qubit states is the Bloch sphere, where it takes the idea of a state vector and turns it into a literal vector. The *z* axis represents the probability of measuring the qubit to be |0⟩ or |1⟩. The closer the vector is to pointing directly at the top of the sphere, the higher the probability of measuring |0⟩, and the closer the vector is to pointing to the bottom of the sphere, the higher the probability of measuring |1⟩. The *x* axis represents the real component of that probability, and the *y* axis represents the imaginary component of that probability. With this, a vector on the Bloch sphere can represent any single state or superposition of states. This representation makes the idea of partial phase and bit rotations much easier to visualize.

<figure><img src="image_18.png" /><figcaption>Figure 16. Bloch sphere denoting a variety of states.</figcaption></figure>

**A.2. Base-10 to Binary.** Quantum states can be written as a number in base-10. To represent that number with multiple qubits, the base-10 number can be converted to binary. This binary number can be easily represented with qubits.

*Example A.1.* For example, the number 3 in base 10 written in binary is 011. With 5 qubits it is represented as

$$ (A.1) \qquad |3\rangle = |00011\rangle. $$

**Appendix B. Single Qubit Gates.** The following single qubit gates are commonly used in quantum circuits.

**B.1. I-Gate.** The Identity Gate does nothing to the state it is applied to. In matrix form, it is simply the identity matrix, i.e.,

$$ (B.1) \qquad I = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}. $$

It acts on the states |0⟩ and |1⟩ in the following way:

$$ (B.2) \qquad \begin{aligned} I|0\rangle &= \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix} = |0\rangle, \\ I|1\rangle &= \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \end{bmatrix} = |1\rangle. \end{aligned} $$
**B.2. X-Gate.** The X-gate is a bit-flip, or a 180° rotation about the x axis on the Bloch sphere written as

$$ (B.3) \qquad X = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} $$

in matrix form. It acts on the states $|0\rangle$ and $|1\rangle$ in the following way:

$$ (B.4) \qquad \begin{aligned} X|0\rangle &= \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \end{bmatrix} = |1\rangle, \\ X|1\rangle &= \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix} = |0\rangle. \end{aligned} $$

The X-basis states are $|+\rangle$ and $|-\rangle$.

**B.3. Z-gate.** The Z-gate is a phase-flip, or a 180° rotation about the z axis on the Bloch sphere written as

$$ (B.5) \qquad Z = \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} $$

in matrix form. It acts on the states $|0\rangle$ and $|1\rangle$ in the following way:

$$ (B.6) \qquad \begin{aligned} Z|0\rangle &= \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix} = |0\rangle, \\ Z|1\rangle &= \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \begin{bmatrix} 0 \\ -1 \end{bmatrix} = |-1\rangle. \end{aligned} $$

The Z-Basis states are $|0\rangle$ and $|1\rangle$.

**B.4. Y-gate.** The Y-gate is a phase-flip and a bit-flip, or a 180° rotation about the y-axis on the Bloch sphere written as

$$ (B.7) \qquad Y = \begin{bmatrix} 0 & -i \\ i & 0 \end{bmatrix} $$

in matrix form. It acts on the states $|0\rangle$ and $|1\rangle$ in the following way:

$$ (B.8) \qquad \begin{aligned} Y|0\rangle &= \begin{bmatrix} 0 & -i \\ i & 0 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 0 \\ i \end{bmatrix} = i|1\rangle, \\ Y|1\rangle &= \begin{bmatrix} 0 & -i \\ i & 0 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \begin{bmatrix} -i \\ 0 \end{bmatrix} = -i|0\rangle. \end{aligned} $$

The Y-Basis states are $|i\rangle$ and $|-i\rangle$ denoted as

$$ (B.9) \qquad |i\rangle := \frac{1}{\sqrt{2}}(|0\rangle + i|1\rangle) = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ i \end{bmatrix} \quad \text{and} \quad |-i\rangle := \frac{1}{\sqrt{2}}(|0\rangle - i|1\rangle) = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ -i \end{bmatrix}. $$
**B.5. Rz-Gate.** The $R_z$-gate is a $\theta$ rad rotation about the $z$ axis written as

$$ (B.10) \qquad R_z(\theta) = \begin{bmatrix} e^{-i\theta/2} & 0 \\ 0 & e^{i\theta/2} \end{bmatrix} $$

in matrix form. It acts on the states $|0\rangle$ and $|1\rangle$ in the following way:

$$ (B.11) \qquad \begin{aligned} R_z(\theta)|0\rangle &= \begin{bmatrix} e^{-i\theta/2} & 0 \\ 0 & e^{i\theta/2} \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = e^{-i\theta/2} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = e^{-i\theta/2}|0\rangle, \\ R_z(\theta)|1\rangle &= \begin{bmatrix} e^{-i\theta/2} & 0 \\ 0 & e^{i\theta/2} \end{bmatrix} \begin{bmatrix} 0 \\ 1 \end{bmatrix} = e^{i\theta/2} \begin{bmatrix} 0 \\ 1 \end{bmatrix} = e^{i\theta/2}|1\rangle. \end{aligned} $$

Recall that $|e^{ix}|^2 = 1$, so the states $|0\rangle$ and $e^{-i\theta/2}|0\rangle$ are indistinguishable upon measurement. This is also true for the states $|1\rangle$ and $e^{i\theta/2}|1\rangle$.

**B.6. H-Gate.** The Hadamard Gate puts an individual qubit into an even superposition of $|0\rangle$ and $|1\rangle$, written in matrix form as

$$ (B.12) \qquad H = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}. $$

It acts on the states $|0\rangle$ and $|1\rangle$ in the following way:

$$ (B.13) \qquad \begin{aligned} H|0\rangle &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 1 \end{bmatrix} = \begin{bmatrix} \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} = \frac{1}{\sqrt{2}} (|0\rangle + |1\rangle) = |+\rangle, \\ H|1\rangle &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ -1 \end{bmatrix} = \begin{bmatrix} -\frac{1}{\sqrt{2}} \\ -\frac{1}{\sqrt{2}} \end{bmatrix} = \frac{1}{\sqrt{2}} (|0\rangle - |1\rangle) = |-\rangle. \end{aligned} $$

**Appendix C. The Hadamard Test.** Let $U \in \mathbb{C}^{n \times n}$ be unitary. The real and imaginary Hadamard tests (HT) are defined as the gate sequences in Figure 17.

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Hadamard test for Re$\langle\psi|U|\psi\rangle$ (left) and Im$\langle\psi|U|\psi\rangle$ (right) where $H$ is the Hadamard gate, $S$ is the Phase gate, $U \in \mathbb{C}^{n \times n}$ is a unitary gate. The topmost qubit is the control qubit and the lower qubit(s) are the computation bits representing the eigenstate.</figcaption></figure>

**C.1. Uses of Phase Estimation.** Quantum phase estimation has many different applications to real world problems. Specifically, for a quantum system, the phase of a state is directly correlated with the energy of that state. By simulating a quantum system, the energy levels of a molecule can be estimated as it undergoes some chemical reaction. In addition, phase estimation is also key for factoring large prime numbers such as in Shor's algorithm. This level of efficient factoring will theoretically compromise current public-key cryptosystems in the near future.
**C.2. Step-by-step (real) Hadamard test proof.** Recall that measuring the control qubit in the real HT to be in the $|0\rangle$ or $|1\rangle$ state is represented with a random variable $X_{\text{Re}}$ being $|0\rangle$ or $|1\rangle$. Instead of referencing these random variables directly, we use the shorthand notation

$$ (C.1) \qquad \mathcal{P}_{\text{Re}|0\rangle} := \mathcal{P}(X_{\text{Re}} \in \{|0\rangle\}) \quad \text{and} \quad \mathcal{P}_{\text{Re}|1\rangle} := \mathcal{P}(X_{\text{Re}} \in \{|1\rangle\}), $$

respectively. Similarly, a random variable $X_{\text{Im}}$ being $|0\rangle$ or $|1\rangle$ describes the measurement of the control qubit in the imaginary HT, where we use the shorthand notation

$$ (C.2) \qquad \mathcal{P}_{\text{Im}|0\rangle} := \mathcal{P}(X_{\text{Im}} \in \{|0\rangle\}) \quad \text{and} \quad \mathcal{P}_{\text{Im}|1\rangle} := \mathcal{P}(X_{\text{Im}} \in \{|1\rangle\}). $$

The initial state of the circuit is

$$ (C.3) \qquad |\pi_0\rangle = |0\rangle \otimes |\psi\rangle = |0\rangle|\psi\rangle = \begin{bmatrix} 1 \\ 0 \end{bmatrix} |\psi\rangle. $$

Applying the Hadamard gate to the control qubit transforms it to the $|+\rangle$ state, i.e.,

$$ \begin{align*} |\pi_1\rangle &= (H \otimes I)|0\rangle|\psi\rangle \\ &= H|0\rangle I|\psi\rangle \\ &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} |\psi\rangle \\ &= \frac{1}{\sqrt{2}} \begin{bmatrix} (1 \times 1) + (1 \times 0) \\ (1 \times 1) + (-1 \times 0) \end{bmatrix} |\psi\rangle \\ &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 1 \end{bmatrix} |\psi\rangle \\ &= \frac{1}{\sqrt{2}} (|0\rangle + |1\rangle) |\psi\rangle \\ &= |+\rangle|\psi\rangle. \end{align*} $$

The controlled unitary operator is then performed rotating the control qubit equivalent to the eigenvalue of $|\psi\rangle$, i.e.,

$$ \begin{align*} |\pi_2\rangle &= CU|+\rangle|\psi\rangle \\ &= (|0\rangle\langle 0| \otimes I + |1\rangle\langle 1| \otimes U)|+\rangle|\psi\rangle \\ &= \frac{1}{\sqrt{2}}|0\rangle I|\psi\rangle + \frac{1}{\sqrt{2}}|1\rangle U|\psi\rangle \\ &= \frac{1}{\sqrt{2}}(|0\rangle I + |1\rangle U)|\psi\rangle. \end{align*} $$
Finally, the second Hadamard gate is applied to the circuit resulting in the state

$$
\begin{align*}
|\pi_3\rangle &= H \frac{1}{\sqrt{2}} (|0\rangle I + |1\rangle U) |\psi\rangle \\
&= \frac{1}{\sqrt{2}} (|+\rangle I + |-\rangle U) |\psi\rangle \\
(C.6) \qquad &= \frac{1}{2} [(|0\rangle + |1\rangle) I + (|0\rangle - |1\rangle) U] |\psi\rangle \\
&= \left( \frac{|0\rangle I + |0\rangle U}{2} \right) |\psi\rangle + \left( \frac{|1\rangle I - |1\rangle U}{2} \right) |\psi\rangle \\
&= |0\rangle \frac{I+U}{2} |\psi\rangle + |1\rangle \frac{I-U}{2} |\psi\rangle.
\end{align*}
$$

The final probabilities can be calculated by taking the probability amplitude times its complex conjugate. First looking at the probability of measuring the control qubit to be in the $|0\rangle$ state,

$$
\begin{align*}
\mathcal{P}_{\text{Re}|0\rangle} &= |\langle 0|\pi_3\rangle|^2 \\
&= \left| \langle 0|0\rangle \frac{I+U}{2} |\psi\rangle + \langle 0|1\rangle \frac{I-U}{2} |\psi\rangle \right|^2 \\
(C.7) \qquad &= \left| \frac{I+U}{2} |\psi\rangle \right|^2 \\
&= \frac{\langle \psi|(I+U^\dagger)(I+U)|\psi\rangle}{4} \\
&= \frac{\langle \psi|(I+U+U^\dagger+U^\dagger U)|\psi\rangle}{4}.
\end{align*}
$$

Now taking a similar process for the probability of measuring the control qubit to be in the $|1\rangle$ state,

$$
\begin{align*}
\mathcal{P}_{\text{Re}|1\rangle} &= |\langle 1|\pi_3\rangle|^2 \\
&= \left| \langle 1|0\rangle \frac{I+U}{2} |\psi\rangle + \langle 1|1\rangle \frac{I-U}{2} |\psi\rangle \right|^2 \\
(C.8) \qquad &= \left| \frac{I-U}{2} |\psi\rangle \right|^2 \\
&= \frac{\langle \psi|(I-U^\dagger)(I-U)|\psi\rangle}{4} \\
&= \frac{\langle \psi|(I-U-U^\dagger+U^\dagger U)|\psi\rangle}{4}.
\end{align*}
$$

Using the fact that $U^\dagger U = I$, these probabilities can be further simplified,

$$
(C.9) \qquad \mathcal{P}_{\text{Re}|0\rangle} = \frac{\langle \psi|(2I + U + U^\dagger)|\psi\rangle}{4} = \frac{2 + \langle \psi|(U + U^\dagger)|\psi\rangle}{4},
$$

$$
(C.10) \qquad \mathcal{P}_{\text{Re}|1\rangle} = \frac{\langle \psi|(2I - U - U^\dagger)|\psi\rangle}{4} = \frac{2 - \langle \psi|(U + U^\dagger)|\psi\rangle}{4}.
$$
For a unitary operator $U$, $\langle\alpha|U|\beta\rangle = \langle\beta|U^{\dagger}|\alpha\rangle^*$, therefore, if $\alpha = \beta$ as in this case, $\langle\psi|U|\psi\rangle = \langle\psi|U^{\dagger}|\psi\rangle^*$. Taking the complex conjugate of both sides then gives that $\langle\psi|U|\psi\rangle^* = \langle\psi|U^{\dagger}|\psi\rangle$. Also note that a complex number plus its complex conjugate is equal to two times its real part. Applying this to the probability equations, we obtain

$$ (C.11) \qquad \mathcal{P}_{\text{Re}|0\rangle} = \frac{2 + 2\text{Re}\langle\psi|U|\psi\rangle}{4} = \frac{1}{2}(1 + \text{Re}\langle\psi|U|\psi\rangle), $$

$$ (C.12) \qquad \mathcal{P}_{\text{Re}|1\rangle} = \frac{2 - 2\text{Re}\langle\psi|U|\psi\rangle}{4} = \frac{1}{2}(1 - \text{Re}\langle\psi|U|\psi\rangle). $$

With the probabilities of measuring $|0\rangle$ and $|1\rangle$ on the control qubit, the real component of the phase can be calculated, i.e.,

$$ (C.13) \qquad \mathcal{P}_{\text{Re}|0\rangle} - \mathcal{P}_{\text{Re}|1\rangle} = \frac{1}{2}(1 + \text{Re}\langle\psi|U|\psi\rangle) - \frac{1}{2}(1 - \text{Re}\langle\psi|U|\psi\rangle) = \text{Re}\langle\psi|U|\psi\rangle, $$

$$ (C.14) \qquad 2\mathcal{P}_{\text{Re}|0\rangle} - 1 = 2\left(\frac{1}{2}(1 + \text{Re}\langle\psi|U|\psi\rangle)\right) - 1 = (1 + \text{Re}\langle\psi|U|\psi\rangle) - 1 = \text{Re}\langle\psi|U|\psi\rangle, $$

$$ (C.15) \qquad 1 - 2\mathcal{P}_{\text{Re}|1\rangle} = 1 - 2\left(\frac{1}{2}(1 - \text{Re}\langle\psi|U|\psi\rangle)\right) = 1 - (\text{Re}\langle\psi|U|\psi\rangle - 1) = \text{Re}\langle\psi|U|\psi\rangle. $$

Therefore, the real expectation value can be found from these probabilities in the following three different ways:

$$ (C.16) \qquad \mathcal{P}_{\text{Re}|0\rangle} - \mathcal{P}_{\text{Re}|1\rangle} = 2\mathcal{P}_{\text{Re}|0\rangle} - 1 = 1 - 2\mathcal{P}_{\text{Re}|1\rangle} = \text{Re}\langle\psi|U|\psi\rangle $$

**C.3. Imaginary Hadamard test proof.** The first two states, $|\pi_0\rangle$ and $|\pi_1\rangle$ are the same as for the real test, i.e., $|\pi_0\rangle$ is

$$ (C.17) \qquad \begin{aligned} |\pi_0\rangle &= |0\rangle|\psi\rangle \\ &= \begin{bmatrix} 1 \\ 0 \end{bmatrix} |\psi\rangle, \end{aligned} $$

and $|\pi_1\rangle$ is

$$ (C.18) \qquad \begin{aligned} |\pi_1\rangle &= H|0\rangle|\psi\rangle \\ &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} |\psi\rangle \\ &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 1 \end{bmatrix} |\psi\rangle \\ &= \frac{1}{\sqrt{2}} (|0\rangle + |1\rangle) |\psi\rangle \\ &= |+\rangle |\psi\rangle. \end{aligned} $$
Now the hermitian of the phase gate, $S^\dagger$ is applied to this state, i.e.,

$$
\begin{align}
|\pi_2\rangle &= S^\dagger|+\rangle|\psi\rangle \nonumber \\
&= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 0 \\ 0 & -i \end{bmatrix} \begin{bmatrix} 1 \\ 1 \end{bmatrix} |\psi\rangle \nonumber \\
&= \frac{1}{\sqrt{2}} \begin{bmatrix} (1 \times 1) + (0 \times 1) \\ (0 \times 1) + (-i \times 1) \end{bmatrix} |\psi\rangle \tag{C.19} \\
&= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ -i \end{bmatrix} |\psi\rangle \nonumber \\
&= \frac{1}{\sqrt{2}} (|0\rangle - i|1\rangle) |\psi\rangle \nonumber \\
&= |-i\rangle|\psi\rangle. \nonumber
\end{align}
$$

Performing the controlled $U$,

$$
\begin{align}
|\pi_3\rangle &= CU|-i\rangle|\psi\rangle \nonumber \\
&= (|0\rangle\langle 0| \otimes I + |1\rangle\langle 1| \otimes U)|-i\rangle|\psi\rangle \tag{C.20} \\
&= \frac{1}{\sqrt{2}} (|0\rangle I|\psi\rangle - i|1\rangle U|\psi\rangle) \nonumber \\
&= \frac{1}{\sqrt{2}} (|0\rangle I - i|1\rangle U) |\psi\rangle. \nonumber
\end{align}
$$

Lastly, applying the second Hadamard gate to the state gives

$$
\begin{align}
|\pi_4\rangle &= H \frac{1}{\sqrt{2}} (|0\rangle I - i|1\rangle U) |\psi\rangle \nonumber \\
&= \frac{1}{\sqrt{2}} \frac{1}{\sqrt{2}} (|+\rangle I - i|-\rangle U) |\psi\rangle \nonumber \\
(C.21) \qquad &= \frac{1}{2} [(|0\rangle + |1\rangle)I - i(|0\rangle - |1\rangle)U] |\psi\rangle \nonumber \\
&= \left( \frac{|0\rangle I - i|0\rangle U}{2} + \frac{|1\rangle I + i|1\rangle U}{2} \right) |\psi\rangle \nonumber \\
&= |0\rangle \frac{I - iU}{2} |\psi\rangle + |1\rangle \frac{I + iU}{2} |\psi\rangle. \nonumber
\end{align}
$$
Using the same approach as before for calculating the probabilities, the probability of measuring the control qubit to be in the $|0\rangle$ state for the imaginary test is

$$ (C.22) \qquad \begin{aligned} \mathcal{P}_{\text{Im}|0\rangle} &= |\langle 0|\pi_3\rangle|^2 \\ &= \left|\langle 0|0\rangle \frac{I-iU}{2} |\psi\rangle\right|^2 \\ &= \left|\frac{I-iU}{2} |\psi\rangle\right|^2 \\ &= \frac{\langle\psi|(I+iU^\dagger)(I-iU)|\psi\rangle}{4} \\ &= \frac{\langle\psi|(I-iU+iU^\dagger+U^\dagger U)|\psi\rangle}{4}. \end{aligned} $$

And again for measuring the control qubit to be in the $|1\rangle$ state,

$$ (C.23) \qquad \begin{aligned} \mathcal{P}_{\text{Im}|1\rangle} &= |\langle 1|\pi_3\rangle|^2 \\ &= \left|\langle 1|1\rangle \frac{I+iU}{2} |\psi\rangle\right|^2 \\ &= \left|\frac{I+iU}{2} |\psi\rangle\right|^2 \\ &= \frac{\langle\psi|(I-iU^\dagger)(I+iU)|\psi\rangle}{4} \\ &= \frac{\langle\psi|(I+iU-iU^\dagger+U^\dagger U)|\psi\rangle}{4} \end{aligned} $$

These probabilities can be simplified further using similar steps to the real circuit:

$$ (C.24) \qquad \begin{aligned} \mathcal{P}_{\text{Im}|0\rangle} &= \frac{\langle\psi|(I - iU + iU^\dagger + U^\dagger U)|\psi\rangle}{4} & \mathcal{P}_{\text{Im}|1\rangle} &= \frac{\langle\psi|(I + iU - iU^\dagger + U^\dagger U)|\psi\rangle}{4} \\ &= \frac{2 - \langle\psi|iU|\psi\rangle + \langle\psi|iU^\dagger|\psi\rangle}{4} & &= \frac{2 + \langle\psi|iU|\psi\rangle - \langle\psi|iU^\dagger|\psi\rangle}{4} \\ &= \frac{1}{2} + i \frac{-\langle\psi|U|\psi\rangle + \langle\psi|U^\dagger|\psi\rangle}{4} & &= \frac{1}{2} + i \frac{\langle\psi|U|\psi\rangle - \langle\psi|U^\dagger|\psi\rangle}{4} \\ &= \frac{1}{2} + \frac{2\operatorname{Im}\langle\psi|U|\psi\rangle}{4} & &= \frac{1}{2} - \frac{2\operatorname{Im}\langle\psi|U|\psi\rangle}{4} \\ &= \frac{1}{2}(1 + \operatorname{Im}\langle\psi|U|\psi\rangle), & &= \frac{1}{2}(1 - \operatorname{Im}\langle\psi|U|\psi\rangle). \end{aligned} $$

Following the same rules mentioned for the real test, subtracting $\mathcal{P}_{\text{Im}|1\rangle}$ from $\mathcal{P}_{\text{Im}|0\rangle}$ should give the imaginary part of the expectation value $\operatorname{Im}\langle\psi|U|\psi\rangle$,

$$ (C.25) \qquad \mathcal{P}_{\text{Im}|0\rangle} - \mathcal{P}_{\text{Im}|1\rangle} = \frac{1}{2}(1 + \operatorname{Im}\langle\psi|U|\psi\rangle) - \frac{1}{2}(1 - \operatorname{Im}\langle\psi|U|\psi\rangle) = \operatorname{Im}\langle\psi|U|\psi\rangle. $$
Therefore, the imaginary part of the expectation value can be found in the following three ways:

$$ (C.26) \qquad \mathcal{P}_{\mathrm{Im}|0\rangle} - \mathcal{P}_{\mathrm{Im}|1\rangle} = 1 - 2\mathcal{P}_{\mathrm{Im}|1\rangle} = 2\mathcal{P}_{\mathrm{Im}|0\rangle} - 1 = \mathrm{Im}\langle\psi|U|\psi\rangle. $$

By combining the real and imaginary results, the expectation value $\langle\psi|U|\psi\rangle$ can be estimated.

Since the eigenvalue is of the form $e^{i\theta}$, Euler's formula can be used to convert this to $\cos\theta + i\sin\theta$ where $\cos\theta$ is the real part and $\sin\theta$ is the imaginary part. Thus, $\mathrm{Re}\langle\psi|U|\psi\rangle = \cos\theta$ and $\mathrm{Im}\langle\psi|U|\psi\rangle = \sin\theta$. So, the real probabilities can be written as

$$ (C.27) \qquad \mathcal{P}_{\mathrm{Re}|0\rangle} = \frac{1}{2}(1 + \mathrm{Re}\langle\psi|U|\psi\rangle) = \frac{1}{2}(1 + \cos\theta) $$

$$ (C.28) \qquad \mathcal{P}_{\mathrm{Re}|1\rangle} = \frac{1}{2}(1 - \mathrm{Re}\langle\psi|U|\psi\rangle) = \frac{1}{2}(1 - \cos\theta) $$

Therefore, $\theta$ can be estimated with

$$ (C.29) \qquad \theta = \arccos(2\mathcal{P}_{\mathrm{Re}|0\rangle} - 1) = \arccos(1 - 2\mathcal{P}_{\mathrm{Re}|1\rangle}). $$

The same can be done with the imaginary part of Euler's formula and the expectation value giving

$$ (C.30) \qquad \mathcal{P}_{\mathrm{Im}|0\rangle} = \frac{1}{2}(1 + \mathrm{Im}\langle\psi|U|\psi\rangle) = \frac{1}{2}(1 + \sin\theta), $$

$$ (C.31) \qquad \mathcal{P}_{\mathrm{Im}|1\rangle} = \frac{1}{2}(1 - \mathrm{Im}\langle\psi|U|\psi\rangle) = \frac{1}{2}(1 - \sin\theta). $$

Therefore, $\theta$ can also be estimated with,

$$ \theta = \arcsin(2\mathcal{P}_{\mathrm{Im}|0\rangle} - 1) = \arcsin(1 - 2\mathcal{P}_{\mathrm{Im}|1\rangle}). $$

Finally, $\theta$ can be estimated using both the real and imaginary expectation values. The real and imaginary expectation values can be thought of as vectors on the complex unit circle with the real part on the x axis and the imaginary on the y axis. Taking $\arcsin$ of the imaginary part or $\arccos$ of the real part effectively finds the vector of the entire complex expectation value. This means that arctan of the imaginary part over the real part also gives $\theta$. Therefore, the angle $\theta$ can be estimated with the following three methods:

$$ \theta = \arccos(\mathrm{Re}\langle\psi|U|\psi\rangle) = \arcsin(\mathrm{Im}\langle\psi|U|\psi\rangle) = \arctan\left(\frac{\mathrm{Im}\langle\psi|U|\psi\rangle}{\mathrm{Re}\langle\psi|U|\psi\rangle}\right). $$

Note that using both the real and imaginary parts to estimate $\theta$ increases the error of the estimation as this uses data from two random variables rather than just one.
Measurement probabilities for the real HT

<figure><img src="image_20.png" /><figcaption>Figure 18. Probabilities of measuring 0 and 1 for the real Hadamard test</figcaption></figure>

**C.4. HT Probabilities.** Cycling through values of $\theta$ from 0 rad to $2\pi$ rad, by plotting the probabilities of measuring $|0\rangle$ and $|1\rangle$ for the real HT in Figure 18, it is clear that they always sum up to 1.

**Appendix D. Building the QFT algorithm.** To try to construct this algorithm independently, we will try to find a pattern to build a circuit from the equation. Suppose that we wanted to transform a single qubit into the Fourier basis. Using the tensor product definition for the final state of a QFT gives us the states

$$ (D.1) \quad \text{QFT}_1|j\rangle = \frac{1}{\sqrt{2^n}} \bigotimes_{l=1}^{1} (|0\rangle + e^{\frac{2\pi i}{2^l}j}|1\rangle) = \frac{1}{\sqrt{2^1}} (|0\rangle + e^{\frac{2\pi i}{2^1}}|1\rangle) = \frac{1}{\sqrt{2}} (|0\rangle + e^{j\pi i}|1\rangle), $$

$$ (D.2) \quad \text{QFT}_1|0\rangle = \frac{1}{\sqrt{2}} (|0\rangle + e^0|1\rangle) = \frac{1}{\sqrt{2}} (|0\rangle + |1\rangle) = |+\rangle, $$

$$ (D.3) \quad \text{QFT}_1|1\rangle = \frac{1}{\sqrt{2}} (|0\rangle + e^{\pi i}|1\rangle) = \frac{1}{\sqrt{2}} (|0\rangle - |1\rangle) = |-\rangle. $$

Like we saw before, this means that QFT acting on the first qubit is equivalent to a Hadamard Gate (Figure 19).

$$ q_0 - \boxed{H} - $$

<caption>Figure 19. Quantum Fourier Transform Circuit for a single qubit.</caption>
Suppose now that we wanted to do a Fourier transform with two qubits instead of one. Let $S$ be the binary number that represents the state to which QFT is applied. This means that $S = q_02^0 + q_12^1 = q_0 + 2q_1$ for two qubits. So, QFT acting on 2 qubits is

$$ (D.4) \qquad \begin{aligned} \text{QFT}_2|q_0q_1\rangle &= \frac{1}{\sqrt{2}} \bigotimes_{l=1}^{2} (|0\rangle + e^{\frac{2\pi i}{2}S}|1\rangle) \\ &= \frac{1}{\sqrt{2}} (|0\rangle + e^{\frac{2\pi i}{2}S}|1\rangle) \otimes \frac{1}{\sqrt{2}} (|0\rangle + e^{\frac{2\pi i}{2}S}|1\rangle) \\ &= \frac{1}{\sqrt{2}} (|0\rangle + e^{\pi i S}|1\rangle) \otimes \frac{1}{\sqrt{2}} (|0\rangle + e^{\frac{\pi i}{2}S}|1\rangle). \end{aligned} $$

Substituting in $S = q_0 + 2q_1$,

$$ (D.5) \qquad \begin{aligned} \text{QFT}_2|q_0q_1\rangle &= \frac{1}{\sqrt{2}} (|0\rangle + e^{q_0\pi i} e^{q_12\pi i} |1\rangle) \otimes \frac{1}{\sqrt{2}} (|0\rangle + e^{q_0\frac{\pi i}{2}} e^{q_1\pi i} |1\rangle) \\ &= \frac{1}{\sqrt{2}} (|0\rangle + e^{q_0\pi i} |1\rangle) \otimes \frac{1}{\sqrt{2}} (|0\rangle + e^{q_0\frac{\pi i}{2}} e^{q_1\pi i} |1\rangle). \end{aligned} $$

The state can be rewritten using $\text{QFT}_1$ as

$$ (D.6) \qquad \text{QFT}_2|q_0q_1\rangle = \text{QFT}_1|q_0\rangle \otimes \frac{1}{\sqrt{2}} (|0\rangle + e^{q_0\frac{\pi i}{2}} e^{q_1\pi i} |1\rangle). $$

Suppose $q_0 = 0$, the state then becomes

$$ (D.7) \qquad \text{QFT}_2|0q_1\rangle = \text{QFT}_1|0\rangle \otimes \frac{1}{\sqrt{2}} (|0\rangle + e^{q_1\pi i} |1\rangle). $$

If $q_1 = 0$, then $q_1$ goes from the $|0\rangle$ state to the $|+\rangle$ state after the transformation. If $q_1 = 1$, then $q_1$ goes from the $|1\rangle$ state to the $|-\rangle$ state after the transformation. This means that we must apply a Hadamard gate to $q_1$.

Suppose now $q_1 = 0$, the state then becomes

$$ (D.8) \qquad \text{QFT}_2|q_00\rangle = \text{QFT}_1|q_0\rangle \otimes \frac{1}{\sqrt{2}} (|0\rangle + e^{q_0\frac{\pi i}{2}} |1\rangle). $$

If $q_0 = 0$, then $q_1$ goes from the $|0\rangle$ state to the $|+\rangle$ state after the transformation. If $q_0 = 1$, then the $|1\rangle$ part of the $q_1$ state rotates by $\frac{\pi}{2}$ rad. This can be added to the circuit with a controlled phase gate. The phase gate is defined as

$$ P(\theta) = \begin{bmatrix} 1 & 0 \\ 0 & e^{i\theta} \end{bmatrix}. $$

Therefore, applying this gate to a superposition of $|0\rangle$ and $|1\rangle$ will rotate the $|1\rangle$ part by $\theta$. By controlling this phase gate with $q_0$ targeting $q_1$ where $\theta = \frac{\pi}{2}$, the $|1\rangle$ part of the state of $q_1$ will rotate by $\frac{\pi}{2}$ rad only in the case that $q_0 = 1$. These rotations are applied as gates in a circuit as shown in Figure 20.
<figure><img src="image_21.png" /><figcaption>Figure 20. Quantum Fourier Transform Circuit for two qubits.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 21. 2 qubit QFT circuit (left) and expected (right) state vector outputs for $q_0 = 1$ and $q_1 = 0$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 22. Quantum Fourier Transform Circuit for two qubits with swaps.</figcaption></figure>

Now that the circuit has been constructed, let us look at the resulting state vector where $q_0 = 1$ and $q_1 = 0$ on the Bloch sphere and compare it to what we expect.

Looking at the two different sets of Bloch spheres in Figure 21, it appears that the order of the qubits is opposite to what we expect. This can be easily resolved by implementing a swap gate to the end of our circuit as shown in Figure 22.

Now that we know how to transform the first 2 qubits, what happens if we try to transform the third? Let $S = q_02^0 + q_12^1 + q_22^2 = q_0 + 2q_1 + 4q_2$ since there are now 3 qubits. QFT acting on 3 qubits is

$$ (D.9) \qquad \begin{aligned} \text{QFT}_3|q_0q_1q_2\rangle &= \frac{1}{\sqrt{2^3}} \bigotimes_{l=1}^{3} (|0\rangle + e^{\frac{2\pi i}{2}S}|1\rangle) \\ &= \frac{1}{\sqrt{2}} (|0\rangle + e^{\pi i S}|1\rangle) \otimes \frac{1}{\sqrt{2}} (|0\rangle + e^{\frac{\pi i}{2}S}|1\rangle) \otimes \frac{1}{\sqrt{2}} (|0\rangle + e^{\frac{\pi i}{4}S}|1\rangle). \end{aligned} $$

This can be rewritten using $\text{QFT}_2$ as

$$ (D.10) \qquad \text{QFT}_3|q_0q_1q_2\rangle = \text{QFT}_2|q_0q_1\rangle \otimes \frac{1}{\sqrt{2}} (|0\rangle + e^{\frac{\pi i}{4}S}|1\rangle). $$
Substituting in $S = q_0 + 2q_1 + 4q_2$, the state becomes

$$ (D.11) \qquad \begin{aligned} \mathrm{QFT}_3|q_0q_1q_2\rangle &= \mathrm{QFT}_2|q_0q_1\rangle \otimes \frac{1}{\sqrt{2}}(|0\rangle + e^{\frac{\pi i}{4}q_0}e^{\frac{\pi i}{4}2q_1}e^{\frac{\pi i}{4}4q_2}|1\rangle) \\ &= \mathrm{QFT}_2|q_0q_1\rangle \otimes \frac{1}{\sqrt{2}}(|0\rangle + e^{q_0\frac{\pi i}{4}}e^{q_1\frac{\pi i}{2}}e^{q_2\pi i}|1\rangle). \end{aligned} $$

So, to transform $q_2$ into the Fourier basis, we must apply a rotation of $\frac{\pi}{4}$ rad if and only if $q_0 = 1$, $\frac{\pi}{2}$ rad if and only if $q_1 = 1$, and $\pi$ rad if and only if $q_2 = 1$ to the $|1\rangle$ part of $q_2$. These transforms can be easily done with controlled phase operations and one Hadamard gate. The implementation of these transformations in a circuit is shown in Figure 23.

<figure><img src="image_24.png" /><figcaption>Figure 23. Quantum Fourier Transform Circuit for three qubits.</figcaption></figure>

Transforming the number 3, where $q_0 = 1$, $q_1 = 1$ and $q_2 = 0$, we can again observe the output state vector in the left panel of Figure 24 to see if it matches what we expect in the right panel of Figure 24.

<figure><img src="image_25.png" /><figcaption>Figure 24. 3 qubit QFT circuit (left) and expected (right) state vector outputs for $q_0 = 1$, $q_1 = 1$, and $q_2 = 0$.</figcaption></figure>

It appears that the qubits are again in reverse order; thus, we still need the swap gate as shown in Figure 25.

<figure><img src="image_26.png" /><figcaption>Figure 25. Quantum Fourier Transform Circuit for three qubits with swaps.</figcaption></figure>

With the swaps, it is clear that the output state vector in Figure 26 is the same as expected.
<figure><img src="image_27.png" /><figcaption>Figure 26. 3 qubit QFT circuit output state vector with swaps.</figcaption></figure>

The pattern that follows from these examples is

$$ (D.12) \qquad \begin{aligned} & \mathrm{QFT}_n |q_0 \cdots q_{n-1}\rangle \\ & = \mathrm{QFT}_{n-1} |q_0 \cdots q_{n-2}\rangle \otimes \frac{1}{\sqrt{2}} \left( |0\rangle + e^{\frac{\pi i}{2^{n-1}} q_0} e^{\frac{\pi i}{2^{n-2}} q_1} \cdots e^{\pi i q_{n-1}} |1\rangle \right). \end{aligned} $$

This means that to transform $q_n$, we can apply a Hadamard gate to a single qubit, then apply a controlled phase gate for each qubit with rotation $\frac{\pi}{2^{n-m}}$ where $n$ is the number of the target qubit and $m$ is the number of the control qubit. Then the order of the qubits must be reversed, which can be done with swap gates.

Now that we have constructed a procedure for an $n$- dimensional QFT, we can transfer this procedure into a for loop to make any size QFT circuit.

**Appendix E. Trotter Approximation.** We are trying to determine the error of approximating

$$ (E.1) \qquad e^{-itH} = \lim_{L \to \infty} \left( e^{-itA/L} e^{-itB/L} \right)^L, $$

using a finite order expansion constant, i.e.,

$$ (E.2) \qquad e^{-itH} = \left( e^{-itA/L} e^{-itB/L} \right)^L. $$

To determine the error of such an approximation, we start by finding the error of a simpler expression,

$$ (E.3) \qquad e^{-itH/L} = e^{-itA/L} e^{-itB/L}. $$

Defining $\delta$ as the difference between the terms in (E.3),

$$ (E.4) \qquad \delta = \left| e^{-itH/L} - e^{-itA/L} e^{-itB/L} \right|, $$

we can Taylor expand the exponentials,

$$ (E.5) \qquad \delta = \left| \sum_{k=0}^{\infty} \frac{1}{k!} \left( -\frac{it}{L} H \right)^k - \sum_{k=0}^{\infty} \frac{1}{k!} \left( -\frac{it}{L} (A+B) \right)^k \right|. $$
Writing out terms,

$$ (E.6) \qquad \delta = \left| \frac{1}{0!} (I - I) - \frac{1}{1!} \frac{it}{L} (H - H) + \mathcal{O} \left( \frac{1}{2!} \frac{t^2}{L^2} (H^2 - (A+B)^2) \right) \right| = \mathcal{O} \left( \frac{t^2}{L^2} \right). $$

Returning to determining the error of the finite expansion, let $\epsilon$ be the difference of the first two terms in (E.2), i.e.,

$$ (E.7) \qquad \epsilon = \left| e^{-itH} - \left( e^{-itA/L} e^{-itB/L} \right)^L \right| $$

Notice that $\epsilon$ follows a form similar to $\delta$, but is raised to the power of $L$. Thus, it is safe to assume that the order of $\epsilon$ is $L$ times the order of $\delta$ in (E.6). Therefore, to obtain precision $\epsilon$,

$$ (E.8) \qquad \epsilon = \mathcal{O} \left( \frac{t^2}{L} \right) \implies L = \mathcal{O} \left( \frac{t^2}{\epsilon} \right). $$

Note that even with a perfect decomposition to Pauli operators, there still exists error in the simulation of the Hamiltonian if the Pauli operators do not commute.