# Simulating spectroscopy experiments with a superconducting quantum computer

John P. T. Stenger  NRC Postdoctoral Associate, U.S. Naval Research
Laboratory, Washington, DC 20375, United States  Gilad Ben-Shach  IBM Quantum,
IBM Canada, Markham, Ontario, L6G 1C7, Canada  David Pekker  Department of
Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA 15260, USA
Nicholas T. Bronn  IBM Quantum, IBM T.J. Watson Research Center, Yorktown
Heights, New York 10598, USA

###### Abstract

We present a novel method for solving eigenvalue problems on a quantum
computer based on spectroscopy. The method works by coupling a “probe” qubit
to a set of system simulation qubits and then time evolving both the probe and
the system under Hamiltonian dynamics. In this way, we simulate spectroscopy
on a quantum computer. We test our method on the IBM quantum hardware for a
simple single spin model and an interacting Kitaev chain model. For the Kitaev
chain, we trace out the pseudo-topological phase boundary for a two-site
model.

The simulation of quantum systems with a quantum computer was originally
proposed by Feynman Feynman (1982) which sparked the nascent field. While the
dynamics of quantum systems are fully determined (e.g. by the Schrödinger
equation), the dimension of the Hilbert space that describes their evolution
scales exponentially in the number of degrees of freedom which makes classical
calculations impossible for a large enough system size. The recent
availability of quantum computers on the cloud allows practitioners to begin
to map problems of quantum simulation and achieve proof-of-principle results
Stenger  _et al._ (2021); Xiao  _et al._ (2021); Rančić ; Sung  _et al._ ;
Harle  _et al._ ; Mi  _et al._ , even starting to surpass the accuracy of
state-of-the-art classical approximations Kim  _et al._ . While quantum phase
estimation (QPE) Kitaev  offers an exponential speed-up in determining
Hamiltonian spectra, it requires a depth of circuit unfeasible for near-term
noisy quantum hardware. Short-depth hybrid algorithms such as the variational
quantum eigensolver (VQE) are amenable to current hardware at the expense of
many iterations over measurements and optimization cycles Peruzzo  _et al._
(2014); Kandala  _et al._ (2017); McClean  _et al._ (2016); O’Malley  _et al._
(2016); Wang  _et al._ (2019); Gonthier  _et al._ (2020); McArdle  _et al._
(2020). Furthermore, VQE requires a problem-specific ansatz which is not
generalizable to all systems. Here we propose a spectroscopic eigensolver
method which makes an efficient use of available noisy near-term quantum
resources. Our method is both amenable to noisy quantum hardware and is
generally applicable to determining Hamiltonian spectra.

Figure 1: Quantum circuits for simulating the time evolution $U_{\rm
Pauli}=e^{-iH_{\rm Pauli}t}$: (a) a cartoon analogy of the spectroscopy
eigensolver method to scanning tunneling spectroscopy, where a tip (qubit) at
potential (energy) $\omega$ interacts with a single mode of a model system via
electron tunneling (transverse $cXX$ interaction), revealing (b) a response as
a function of energy which determines an eigenvalue $\lambda$ as a dip in
probe expectation $\langle Z_{0}\rangle$. (c) The spectroscopic eigensolver
algorithm for a probe qubit interacting with a single qubit of the simulation
register (for $N_{t}$ Trotter steps $U_{\rm res}$, first-order shown for
clarity) followed by measurement of the probe qubit, and (d) a circuit
equivalence that allows the efficient implementation of scaled pulses with the
native $R^{zx}(\theta)$ interaction generated by echoed cross resonance
Sheldon  _et al._ (2016).

We present a quantum eigensolver method which is akin to the experimental tool
of tunneling spectroscopy. Our method involves a “probe qubit” and a
“simulation” register that undergo Hamiltonian evolution for a time $t$. A
cartoon of our method is depicted in Fig. 1a. Measurement of the “probe” qubit
$q_{0}$ reveals the spectrum of the evolved Hamiltonian via its response when
the probe energy is on-resonance with an energy transition in the system,
schematically shown as a dip in Fig. 1b. While similar to other time-evolution
algorithms Kitaev ; Somma  _et al._ (2002), this technique does not require
controlled-coherent evolution of the register by auxiliary qubits. Instead,
the probe qubit evolving according to an energy (via $\omega$) interacts
tranversely with a single qubit of the simulation register. This process is
repeated over a number of Trotter steps as shown in Fig. 1c. We have the
freedom to select the specific qubit being probed. It is only necessary that
the wavefunction of the simulation register has strong weight on that qubit.
This method is conducive for execution on near-term noisy quantum hardware
because it only requires a local interaction between qubits instead of the
controlled-evolution of an entire simulation register, substantially reducing
circuit depth because no SWAP operations are necessary between probe qubit and
simulation register and naturally maps to the planar architecture of current
superconducting quantum hardware (see Supplement I). Furthermore, only
measurement of the probe qubit is required, reducing readout error. The
underlying trade-off as compared to other algorithms is the number of circuits
that must be executed: one for each value of $\omega$, a sweep of which is
required for accurate fitting and hence determination of the eigenvalues (see
Supplement II). However, recent improvements as determined by speed
benchmarking show this is not a bottleneck for superconducting quantum
hardware Wack  _et al._ (2021).

Explicitly, the spectroscopic eigensolver technique finds the energy
differences of a system by coupling a probe qubit to a Pauli Hamiltonian
$H_{\rm Pauli}$ suitably encoded to enforce the commutation relations of the
system Hamiltonian $H_{\rm sys}$ it represents. The full Hamiltonian, which we
will refer to as the “resonance” Hamiltonian because of its response when the
probe and system are on-resonance, is

$H_{\rm res}=-\frac{1}{2}\omega Z_{0}+cX_{0}X_{i}+H_{\rm Pauli},$ (1)

where the probe qubit $q_{0}$ with energy $\omega$ is at index 0, $c$ is the
coupling parameter between the probe and probed qubit $q_{i}$ of the simulated
system, and the tensor products with identity matrices are omitted for
succinctness. The first-order Suzuki-Trotter decomposition of the time
evolution unitary is

$U_{\rm res}(\omega,dt)=e^{-iH_{\rm res}dt}\approx
R^{xx}_{0i}\left(2c\,dt\right)R^{z}_{0}\left(-\omega\,dt\right)U_{\rm Pauli},$
(2)

where $R^{z}_{j}(\phi)$ is a $Z$-rotation on site $j$, $R^{xx}_{jk}(\theta)$
is an $XX$-rotation on sites $j$ and $k$ and $U_{\rm Pauli}=e^{-iH_{\rm
Pauli}dt}$ is the Trotter step of time evolution $dt=t/N_{t}$ for the system.

We start in a state $\ket{0,\psi_{0}}$ with the probe qubit in the ground and
the system $\ket{\psi_{0}}$ is in an arbitrary state. For simplicity we will
take $|\psi_{0}\rangle$ to be the $n$-qubit ground state. After $N_{t}$
applications of $U_{\rm res}$, we measure the probe qubit in the $Z$-basis,
$\braket{Z_{0}}(\omega)=\bra{0,\psi_{0}}(U_{\rm
res}^{\dagger})^{Nt}Z_{0}U_{\rm res}^{Nt}\ket{0,\psi_{0}}$. When $c$ is small
compared to the energy scales in $H_{\rm Pauli}$ and when $\omega$ is on
resonance with an energy transition of $H_{\rm Pauli}$, the probability of the
probe qubit flipping will peak, resulting in a dip in
$\braket{Z_{0}}(\omega)$.

Figure 2: Scans of the probe qubit expectation value as a function of $\omega$
for many values of $a$. The minimums in $\braket{Z_{0}}(\omega)$ for each
value of $a$ are plotted as red circles. The error bars represent the full-
width-at-half-minimum for the $\braket{Z_{0}}(\omega)$ dips. The black dashed
lines are the expected value of the energy transitions as obtained from exact
diagonalization. The parameter values are $c=0.1$, $t=10$, $dt=0.3\bar{3}$.

The following experiments were conducted within the Qiskit Qis (2021)
framework with jobs sent to the cloud-based ibm_lagos, a 7-qubit
superconducting IBM Quantum backend. First, the system Hamiltonian is mapped
to a qubit (Pauli) Hamiltonian by a suitable encoding. The Pauli Hamiltonian
is then mapped to unitary time evolution via the second-order Suzuki-Trotter
transformation. The resulting circuit consists of two-qubit rotations on the
order of $c\,dt$, $\omega\,dt$, etc., where the digital synthesis of each is
locally equivalent to a $Z$-rotation sandwiched by two CNOTs, which is more
efficiently implemented by scaled gates that rotate smaller angles in the two-
qubit Hilbert space natively (the equivalence of which is shown in Fig. 1d)
instead of the (net) two-qubit $\pi$ rotation in the former case Stenger  _et
al._ (2021). These scaled pulses are implemented as calibrations for
$R^{zx}(\theta)$ gates following an analysis of the circuit using a novel
template optimization technique Iten  _et al._ . Each spectroscopic experiment
is then built as an array of circuits, one for each value of $\omega$. For
each circuit, a result consisting of the measurement strings of 8192 shots is
returned and analyzed. This allows us to find the energy transitions by
locating dips in $\braket{Z_{0}}(\omega)$. Further details are provided in the
Supplementary Material and full data and source code can be found at
www.github.com/qiskit-research/qiskit-research.

As a first test case, we take the system to be a single spin under two
opposing magnetic fields,

$H_{\text{Pauli}}=H_{\text{sys}}=aZ_{1}+bY_{1}.$ (3)

We treat $b$ as fixed and vary $a$ so that the energy levels form an avoided
crossing. Both $a$ and $b$ are treated as unitless algorithmic quantities. In
Fig. 2 we show scans of the orientation of the probe qubit
$\braket{Z_{0}}(\omega)$ for several values of $a$. Notice that there is a dip
at $\omega=0$. This corresponds to the transition of an energy level to
itself, see the Supplement III for details.

We find excellent agreement between the exact energy transitions and the local
minimum in $\braket{Z_{0}}(\omega)$. The Root-Mean-Square (RMS) of the
variance of all three transitions together is
$\text{RMS}=\sqrt{\frac{1}{N_{d}}\sum_{i=0}^{N_{d}}(\omega^{*}_{i}-\Delta
E_{i})^{2}}\approx 0.083$ where $N_{d}$ is the number of data points,
$w^{*}_{i}$ is a minimum and $\Delta E_{i}$ is the corresponding energy
transition. The RMS variance is about an order of magnitude less than the
average uncertainty calculated from the Full-Width-at-Half-Minimum (FWHM),
which includes both hardware and Trotter error,
$\text{FWHM}=\frac{1}{N_{d}}\sum_{i=0}^{N_{d}}(\omega^{L}_{i}-\omega^{R}_{i})\approx
0.811$ where $\omega^{L}_{i}$ is the value of $\omega$ where
$\braket{Z_{0}}(\omega)$ is half way between its minimum and maximum value to
the left of $\omega^{*}_{i}$ and $\omega^{R}_{i}$ is the same but to the right
of $\omega^{*}_{i}$. This simple test system stands to demonstrate how
accurate our algorithm can be when the decoherence time of the qubits ($\sim
100-150$ $\mu$s) allows for a long evolution time $t\gg dt$, where $t$
corresponds to $\sim 15\leavevmode\nobreak\ \mu$s of hardware time for this
example. With this in mind, we will push our algorithm to the limits of the
quantum hardware by testing a more complicated system.

A system of particular interest in tunneling spectroscopy is the $n$-site
Kitaev chain which hosts Majorana Zero Modes (MZMs) in the topological regime.
The Kitaev chain is a particularly good test case for the spectroscopic
eigensolver technique since the MZMs exist at the ends of the chain and are,
therefore, easy to probe. Distinguishing the topological phase from the
trivial superconducting phase in such nanowires has been a recent theoretical
interest largely motivated by the fact that Majorana zero modes promise
application in topologically protected quantum computation Kitaev (2003);
Nayak  _et al._ (2008); Alicea  _et al._ (2011); Bonderson and Lutchyn (2011);
Stenger  _et al._ (2019); Ben-Shach  _et al._ (2015). The topological regime
is characterized by a change in parity of the ground state. While the
topological regime is only truly present in a $n\rightarrow\infty$ limit, a
parity flip occurs even in a two site model.

The model Hamiltonian of the interacting Kitaev chain with Coulomb
interactions is

$\begin{split}H_{\rm
sys}=&\mu\sum_{i=1}^{L}c^{\dagger}_{i}c_{i}+g\sum_{i=1}^{L-1}(c^{\dagger}_{i}c_{i+1}+c^{\dagger}_{i+1}c_{i})\\\
+&\Delta\sum_{i=1}^{L-1}(c^{\dagger}_{i}c^{\dagger}_{i+1}+c_{i+1}c_{i})+V\sum_{i=1}^{L-1}c^{\dagger}_{i}c_{i}c^{\dagger}_{i+1}c_{i+1}\end{split}$
(4)

where $L$ is the number of sites, $c^{\dagger}_{i}$ ($c_{i}$) is the creation
(annihilation) operator on site $i$, $\mu$ is the chemical potential, $g$ is
the hopping rate, $\Delta$ controls the superconductivity pairing strength,
and $V$ is the interaction strength. The inclusion of the interaction term
makes the system difficult to study classically but poses no additional
fundamental challenges for quantum processors Childs  _et al._ (2018);
Tacchino  _et al._ (2020).

Using the Jordan-Wigner encoding Somma  _et al._ (2002), we can express Eq. 4
as

$\begin{split}H_{\rm
Pauli}=&x\sum_{i=1}^{n-1}X_{i}X_{i+1}+y\sum_{i=1}^{n-1}Y_{i}Y_{i+1}+z\sum_{i=1}^{n-1}Z_{i}Z_{i+1}\\\
&-m\sum_{i=1}^{n}Z_{i}-\bar{m}\sum_{i=2}^{n-1}Z_{i}\end{split}$ (5)

where $X_{i},Y_{i},Z_{i}$ are Pauli matrices acting on qubit $i$, $n=L$ is the
number of qubits used to represent the system, and $2x=g+\Delta$,
$2y=g-\Delta$, $4z=V$, $4m=2\mu+V$, and $4\bar{m}=V$. These parameters are
considered to be unitless algorithmic quantities, however, the values we will
use roughly correspond to those of real systems if they are taken to be in the
meV range Chen  _et al._ (2017).

In the limit $n\rightarrow\infty$, this system has a topological phase
transition in parameter space. We will restrict ourselves to the two-site
$n=2$ model in the main text. See Supplement VI for an explanation of how to
relate the two site model and the $n\rightarrow\infty$ limit and see
Supplement VIII for experimental results for the $n=3$ model. In the $n=2$
case the gap between the ground and first excited state closes along the
surface

$m=\sqrt{z^{2}+z(x+y)+xy}$ (6)

see Supplement VII for details. This surface can be thought of as a remnant of
the topological phase boundary. We now explore the use of the novel
spectroscopic eigensolver technique to analyze this surface.

Because the circuit depth is long in this case, it is important that we
optimize the parameters. We optimized the time $t$, time step $dt$ and
coupling $c$ (see Supplement IV). In Fig. 3a we show a sweep of the probe
qubit’s energy $\omega$ for three values of $c$. The energy difference between
the ground and first excited state is determined by the $\omega$ at which
$\langle Z_{0}\rangle(\omega)$ is minimum. When $c$ is too low, the probe
qubit will not respond to the system but if c is too high it will perturb the
system. Note the state transition from ground to first excited (denoted by the
notation [0,1]) is visible in Fig. 3a but the reverse transition is not. This
is due to the initial state of the system qubits, set to $\ket{00}$ for
simplicity, which has no overlap with the odd parity excited state and so only
the forward transition is possible. Considering a single parity, in this way,
allows maximum contrast in our signal. Absent also is the $\omega=0$ dip which
is true for all fermionic Hamiltonians (see Supplement III), making them
particularly amenable to our method. At the edges of Fig. 3a we see the onset
of higher energy transitions.

Figure 3: Scans of the probe qubit expectation value as a function of
$\omega$. (a) sweeps of the probe qubit frequency for different values of $c$.
For these sweeps, $m=1.0$, $x=1.5$, $y=0.4$, $z=0.2$, $dt=0.7$ and $t=5$. The
vertical lines are the expected transitions from the classical solution. For
more $c$-sweeps see the Supplement IV. (b-c) sweeps of the probe qubit
frequency for different values of $y$. For these sweeps, $c=0.3$, $m=0.1$,
$x=1.5$, $z=0.4$, $dt=0.7$ and $t=5$. The minimums in $\braket{Z_{0}}(\omega)$
for each value of $y$ are plotted as red circles. The error bars represent the
half-width-at-half-minimum for the $\braket{Z_{0}}(\omega)$ dips. The blue
dashed line represents the [0,1] transition while other transitions are in
grey. In panel (b) we show all transitions and all local minimums of
$\braket{Z_{0}}(\omega)$. In panel (c) we show only the [0,1] and [1,0]
transition and we select the minimums by tracking as depicted by the yellow
arrow in (b).

Scans of $\braket{Z_{0}}(\omega)$ for several values of $y$ are shown in Fig.
3b-c. In panel 3b the location of all local minima are displayed. In 3c we
isolate the [0,1] transition by locating the minimum which is closest to
$\omega=0$ (circled in yellow in panel-b) and then track the dip by going to
the neighboring values of $y$ and finding the minimum which is closest to the
last. The average FWHM $\approx 2.35$ is larger than the separation of some
energy transitions. This can cause dips from two different energy transitions
to merge into one, shifting the location of the dip. Even with this concern,
we get fairly accurate results. The RMS variance $\approx 1.72$ for the [0,1]
transition which is well within the width of the peak.

Figure 4: Energy difference between the ground and first excited state as a
function of $m$ and $y$. (a) Simulation on a classical computer with no error.
(b) results from the ibm_lagos backend. In both plots, $c=0.3$, $x=1.5$,
$z=0.4$, $dt=0.7$, and $t=5.0$. In both panels, the orange line shows the
location of the data in Fig. 3b,c and the blue curve is the expected zero-
crossing from Eq. 6 with the above parameters. The green curve in (b) is the
zero-crossing again using Eq. 6 with $z$ as the free parameter. We find the
best-fit $z$ is shifted by $\Delta z=0.19$.

Taking several scans of $\omega$ for different values of $m$ and $y$, we are
able to trace out the pseudo phase boundary. Fig. 4 shows the absolute value
of the gap for (a) a classical simulation without noise and (b) the results
from the quantum device. The curve of points at which the gap goes to zero is
the phase boundary. Notice that in the results from the quantum device the
boundary is pushed to higher values of $y$ than in the perfect simulation. The
green line is a best fit of the boundary to Eq. 6 with z as the fit parameter.
We find a shift of $\Delta z=0.19$ which is well within the average FWHM
$\approx 2.35$ found above. Using the best fit $z$ we find that the shape of
the boundary is fairly accurate with RMS $\approx 0.15$ which is again well
within the average FWHM. The error is largely due to the width of the
$\braket{Z_{0}}(\omega)$ dips and the influence of nearby transitions. As seen
in Fig. 3 the [0,1] dip is pulled to higher $\omega$ by the [2,3] transition
near the phase boundary while it is pulled down by the [1,0] transition after
the phase boundary. This causes the boundary to appear to be shifted to higher
values of $y$ (or $z$ if $y$ is taken as fixed). There may also be an effect
due to an unintended $ZZ$ rotation, which is a well-known error for quantum
computers built from fixed-frequency transmons Kandala  _et al._ (2020);
Sundaresan  _et al._ (2020). While it is difficult to disentangle all of the
sources of error, our goal is not to remove the error completely but to
demonstrate the procedure for finding the phase boundary using our
spectroscopic eigensolver. The accuracy is determined by the width of the
dips, set by the time $t$ for which the algorithm can run, which is limited by
the number of quantum gates which can be applied before decoherence becomes
large. As quantum hardware improves, this number of gates will increase and
the accuracy of our algorithm will improve.

While the energy transitions are symmetric between positive and negative
values of $m$, it is often the case that a particular transition is more
apparent in the data for either $\pm$ values of $m$. Fig. 4 is generated by
taking the most accurate data between positive an negative $m$. For the raw
data and a detailed discussion of how the data was filtered see Supplement V.

The experiments discussed here demonstrate the ability of the novel
spectroscopic eigensolver technique to solve for eigenspectra on near term
quantum devices. Similar to the experimental method of tunneling spectroscopy,
the simulation register interacts with a local probe qubit to determine its
spectrum. This requires only the realization of two-qubit interactions between
the probe and register without the coherent control of the register at the
heart of many other simulation algorithms, making this technique highly
amenable to exploration on near-term noisy quantum computing hardware.

The authors acknowledge the use of IBM Quantum Services for this work. The
authors also thank C.S. Hellberg, L.D. Gunlycke, D.J. Egger, S.M. Frolov, P.
Jurcevic, and Y. Kim for insightful discussions, D.T. McClure, P. Nation, S.
Panda and R. Woo for assistance in performing the experiments, and K.J.
Ferris, O.T. Lanes, and K.J. Sung for careful readings of the manuscript. We
acknowledge the Air Force Research Laboratory (AFRL) for providing additional
quantum resources through its partnership with IBM.

## References

  * Feynman (1982) R. P. Feynman, International Journal of Theoretical Physics 21, 467 (1982).
  * Stenger  _et al._ (2021) J. P. T. Stenger, N. T. Bronn, D. J. Egger, and D. Pekker, Physical Review Research 3, 033171 (2021).
  * Xiao  _et al._ (2021) X. Xiao, J. K. Freericks, and A. F. Kemper, Quantum 5 (2021).
  * (4) M. J. Rančić,  arXiv:2108.07235 .
  * (5) K. J. Sung, M. J. Rančić, O. T. Lanes, and N. T. Bronn,  arXiv:2206.00563 .
  * (6) N. Harle, O. Shtanko, and R. Movassagh,  arXiv:2203.15083 .
  * (7) X. Mi, M. Sonner, M. Y. Niu, K. W. Lee, B. Foxen, R. Acharya, I. Aleiner, T. I. Andersen, F. Arute, K. Arya, A. Asfaw, J. Atalaya, R. Babbush, D. Bacon, J. C. Bardin, J. Basso, A. Bengtsson, G. Bortoli, A. Bourassa, L. Brill, M. Broughton, B. B. Buckley, D. A. Buell, B. Burkett, N. Bushnell, Z. Chen, B. Chiaro, R. Collins, P. Conner, W. Courtney, A. L. Crook, D. M. Debroy, S. Demura, A. Dunsworth, D. Eppens, C. Erickson, L. Faoro, E. Farhi, R. Fatemi, L. Flores, E. Forati, A. G. Fowler, W. Giang, C. Gidney, D. Gilboa, M. Giustina, A. G. Dau, J. A. Gross, S. Habegger, M. P. Harrigan, J. Hilton, M. Hoffmann, S. Hong, T. Huang, A. Huff, W. J. Huggins, L. B. Ioffe, S. V. Isakov, J. Iveland, E. Jeffrey, Z. Jiang, C. Jones, D. Kafri, K. Kechedzhi, T. Khattar, S. Kim, A. Y. Kitaev, P. V. Klimov, A. R. Klots, A. N. Korotkov, F. Kostritsa, J. M. Kreikebaum, D. Landhuis, P. Laptev, K.-M. Lau, J. Lee, L. Laws, W. Liu, A. Locharla, E. Lucero, O. Martin, J. R. McClean, M. McEwen, B. M. Costa, K. C. Miao, M. Mohseni, S. Montazeri, A. Morvan, E. Mount, W. Mruczkiewicz, O. Naaman, M. Neeley, C. Neill, M. Newman, T. E. O’Brien, A. Opremcak, A. Petukhov, R. Potter, C. Quintana, N. C. Rubin, N. Saei, D. Sank, K. Sankaragomathi, K. J. Satzinger, C. Schuster, M. J. Shearn, V. Shvarts, D. Strain, Y. Su, M. Szalay, G. Vidal, B. Villalonga, C. Vollgraff-Heidweiller, T. White, Z. J. Yao, P. Yeh, J. Yoo, A. Zalcman, Y. Zhang, N. Zhu, H. Neven, S. Boixo, A. Megrant, Y. Chen, J. Kelly, V. Smelyanskiy, D. A. Abanin, and P. Roushan,  arXiv:2204.11372 .
  * (8) Y. Kim, C. J. Wood, T. J. Yoder, S. T. Merkel, J. M. Gambetta, K. Temme, and A. Kandala,  arXiv:2108.09197 .
  * (9) A. Y. Kitaev,  arXiv:9511026 .
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. R. McClean, P. Shadbolt, M.-H. H. Yung, X.-Q. Q. Zhou, P. J. Love, A. Aspuru-Guzik, J. L. O’Brien, and J. L. O’Brien, Nature Communications 5, 4213 (2014).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New Journal of Physics 18, 023023 (2016).
  * O’Malley  _et al._ (2016) P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik, and J. M. Martinis, Phys. Rev. X 6, 031007 (2016).
  * Wang  _et al._ (2019) D. Wang, O. Higgott, and S. Brierley, Physical Review Letters 122 (2019), 10.1103/physrevlett.122.140504.
  * Gonthier  _et al._ (2020) J. F. Gonthier, M. D. Radin, C. Buda, E. J. Doskocil, C. M. Abuan, and J. Romero, “Identifying challenges towards practical quantum advantage through resource estimation: the measurement roadblock in the variational quantum eigensolver,”  (2020).
  * McArdle  _et al._ (2020) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Reviews of Modern Physics 92 (2020), 10.1103/revmodphys.92.015003.
  * Sheldon  _et al._ (2016) S. Sheldon, E. Magesan, J. M. Chow, and J. M. Gambetta, Physical Review A 93, 060302 (2016).
  * Somma  _et al._ (2002) R. Somma, G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, Physical Review A 65, 042323 (2002).
  * Wack  _et al._ (2021) A. Wack, H. Paik, A. Javadi-Abhari, P. Jurcevic, I. Faro, J. M. Gambetta, and B. R. Johnson,  (2021), arXiv:2110.14108 .
  * Qis (2021) “Qiskit: An open-source framework for quantum computing,”  (2021).
  * (21) R. Iten, R. Moyard, T. Metger, D. Sutter, and S. Woerner,  arXiv:1909.05270 .
  * Kitaev (2003) A. Y. Kitaev, Annals of Physics 303, 2 (2003).
  * Nayak  _et al._ (2008) C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. S. Das Sarma, Reviews of Modern Physics 80, 1083 (2008).
  * Alicea  _et al._ (2011) J. Alicea, Y. Oreg, G. Refael, F. von Oppen, and M. P. A. Fisher, Nature Physics 7, 412 (2011).
  * Bonderson and Lutchyn (2011) P. Bonderson and R. M. Lutchyn, Physical Review Letters 106, 130505 (2011).
  * Stenger  _et al._ (2019) J. P. T. Stenger, M. Hatridge, S. M. Frolov, and D. Pekker, Physical Review B 99, 035307 (2019).
  * Ben-Shach  _et al._ (2015) G. Ben-Shach, A. Haim, I. Appelbaum, Y. Oreg, A. Yacoby, and B. I. Halperin, Physical Review B 91, 045403 (2015).
  * Childs  _et al._ (2018) A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su, Proceedings of the National Academy of Sciences of the United States of America 115, 9456 (2018).
  * Tacchino  _et al._ (2020) F. Tacchino, A. Chiesa, S. Carretta, and D. Gerace, Advanced Quantum Technologies 3, 1900052 (2020).
  * Chen  _et al._ (2017) J. Chen, P. Yu, J. Stenger, M. Hocevar, D. Car, S. R. Plissard, E. P. A. M. Bakkers, T. D. Stanescu, and S. M. Frolov, Science Advances 3 (2017), 10.1126/sciadv.1701476.
  * Kandala  _et al._ (2020) A. Kandala, K. X. Wei, S. Srinivasan, E. Magesan, S. Carnevale, G. A. Keefe, D. Klaus, O. Dial, and D. C. McKay, Physical Review Letters 127, 130501 (2020).
  * Sundaresan  _et al._ (2020) N. Sundaresan, I. Lauer, E. Pritchett, E. Magesan, P. Jurcevic, and J. M. Gambetta, PRX Quantum 1, 020318 (2020).
  * Lloyd (1996) S. Lloyd, Science 273, 1073 (1996).
  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Córcoles  _et al._ (2021) A. D. Córcoles, M. Takita, K. Inoue, S. Lekuch, Z. K. Minev, J. M. Chow, and J. M. Gambetta, Physical Review Letters 127, 100501 (2021).
  * Somma (2019) R. D. Somma, New Journal of Physics 21, 123025 (2019).
  * (37) “mapomatic: Automatic mapping of compiled circuits to low-noise sub-graphs,” https://github.com/Qiskit-Partners/mapomatic, accessed: 2022-03-30.
  * Roushan  _et al._ (2017) P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas, A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Fowler, B. Foxen, M. Giustina, E. Jeffrey, J. Kelly, E. Lucero, J. Mutus, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. White, H. Neven, D. G. Angelakis, and J. Martinis, Science 358, 1175 (2017).

##  I Hamiltonian Simulation Algorithms

Quantum computers will theoretically give an exponential speed-up to the
problem of solving for the eigenvalues of Hamiltonians by the celebrated
quantum phase estimation (QPE) algorithm Kitaev  – the backbone of proposed
universal quantum simulation Lloyd (1996). QPE estimates the phase of a
unitary operator $U$ acting on an $n$-dimensional Hilbert space encoded by a
“simulation” register to within an arbitrary precision specified by the number
$m$ of qubits in an auxiliary register. Performing QPE involves realizing
controlled-$U^{k}$ gates (Trotterized into $N_{t}$ steps each) on $n$ qubits
(with $1\leq k\leq 2^{m-1}$) followed by an inverse quantum Fourier transform
on the auxiliary register (as depicted in Fig. 5a), and is prohibitive for
near-term noisy quantum hardware Preskill (2018). Current research thrusts
involve reducing the quantum resources necessary for performing quantum
simulations while still seeking exponential speed-up. One such method known as
iterative phase estimateion (IPE) replaces the $m$-qubit auxiliary register
with a recycled “pointer” qubit via mid-circuit measurements and feed-forward
Córcoles  _et al._ (2021). Another method involves performing a classical
Fourier transform on the time series generated by expectation values of the
phase kickback on the auxiliary qubit from a controlled-$U$ on the simulation
register Somma (2019). The spectroscopic eigensolver (Fig. 5b for comparison),
however, uses only a local interaction and does not require coherent control
of the register to elucidate the eigenvalues. As a demonstration of the
expected performance, we will estimate a fidelity score for the resource
extremal case of QPE.

(a)

$\Qcircuit<EMAIL_ADDRESS>\cdots\cdots\nghost{QFT^{\dagger}}\cdots\\\
\lstick{\ket{0}}\gate{H}\qw\qw\cdots\ctrl{1}\ghost{QFT^{\dagger}}\meter\\\
\lstick{\ket{\psi_{0}}}{/}\qw\gate{U_{\rm
Pauli}^{N_{t}2^{m-1}}}\qw\cdots\gate{U_{\rm Pauli}^{N_{t}2^{0}}}\qw\qwa}$

(b)

$\Qcircuit@C=1em@R=0.8em{\lstick{\ket{0}}&\qw\qw\gate{R_{0}^{z}(-\omega\,dt)}\qw\multigate{1}{R_{0i}^{xx}(2c\,dt)}\qw\qw\meter\\\
\lstick{}\qw\qw\multigate{2}{U_{\rm
Pauli}}\qw\ghost{R_{0i}^{xx}(2c\,dt)}\qw\qw\qwa\\\ \cdots\nghost{U_{\rm
Pauli}}\cdots\cdots\\\ \lstick{}\qw\qw\ghost{U_{\rm
Pauli}}\qw\qw\qw\qw\qwa\inputgroupv{2}{4}{0.7em}{1.69em}{\ket{\psi_{0}}}\gategroup{1}{3}{4}{6}{2.4em}{(}\gategroup{1}{3}{4}{6}{2.4em}{)}}\raisebox{14.63881pt}{\makebox[-21.0pt][r]{$N_{t}$}}$

Figure 5: Comparison of (a) QPE and (b) the spectroscipic eigensolver
algorithms .

The accuracy of the spectroscopic eigensolver on the Landau-Zener model is
approximately equivalent to three bits of precision ($\sim 10\%$) and is
obtained with a single qubit simulation register. The same precision in QPE
requires a three-qubit measurement register. Using the PhaseEstimation circuit
from Qiskit with the same Trotterization, SWAP-mapping determined by the
(stochastic) SABRE SWAP method of the Qiskit transpiler, the resulting circuit
is mapped to the least noisy qubits on ibm_lagos with the mapomatic package
map . The mapomatic package provides an infidelity score calculated from
single- and two-qubit gate error and measurement error, with each $R^{zx}$
gate error calculated as linear in angle (as found in Ref. Stenger  _et al._
(2021)) with a floor of the $X$-error that forms the echo. These are reported
in Table 1. We see that QPE is estimated to have three-9’s of infidelity while
the spectroscopic eigensolver is around 65%, which improves to 6% when pulse-
scaled $R^{zx}$ gates are used. Note that while these exact values are based
on current error rates determined by calibration and are subject to frequent
change, they provide an estimate to the required quantum resources.

| QPE | Spec. Eigen. | Spec. Eigen. w/ $R^{zx}$
---|---|---|---
$N_{\rm qubits}$ | 4 | 2 | 2
#CNOTs | 1272 | 120 | $60^{*}$
Infid. Score | 0.99984 | 0.65363 | 0.06018
Table 1: A comparison of parameters between QPE and the spectroscopic
eigensolver for the Landau-Zener model, and that implemented with
$R^{zx}(\theta)$ operations, the count of which is denoted in the table by an
asterisk (*).

##  II Functional form of the resonance eigensolver

The algorithm works by applying unitary time evolution to an arbitrary state
(in our simulation we choose $|\psi_{0}\rangle=|00\rangle$) of the system with
the probe qubit in the ground state. We then observe the probability that the
probe qubit becomes excited. If the energy on probe qubit (as swept by
$\omega$ in the resonance Hamiltonian, not the actual transition energy of
transmon as in Roushan  _et al._ (2017)) matches an energy transition in the
system then we will find a peak in the probability. To demonstrate that the
probability indeed peaks at energy transitions we will use two approaches.
First we will derive the first order perturbation of the probability as a
function of probe qubit frequency. This will show us that the probability
decays as the frequency of the probe qubit is taken off resonance. Second we
will analyze the behavior of the probability when the system is composed of
only two energy levels. Since the probability decays away from resonance, a
two level system is a good approximation for a larger system where a
particular transition is close to resonance with the probe qubit.

Let us start with the perturbation approach. In the algorithm, we apply
unitary time evolution

$U_{\rm res}(\omega)=e^{-i(-\omega Z_{0}/2+H_{\rm Pauli}+cX_{0}X_{i})dt}$ (7)

where $H_{\rm Pauli}$ is the qubit-encoded system Hamiltonian, qubit zero is
the probe qubit, and qubit $i$ is the qubit being probed. We desire $c$ to be
small so that we can treat the $cX_{0}X_{i}$ term as a perturbation. Let
$E_{n}$ and $\ket{n}$ be the energies and eigenvectors of the system
Hamiltonian $H_{\rm Pauli}\ket{n}=E_{n}\ket{n}$ so that the unperturbed
eigenvalues are $E_{a,n}=(-1)^{a+1}\omega/2+E_{n}$ and the unperturbed
eigenvectors are $\ket{a,n}$ where $a\in\\{0,1\\}$ labels the state of the
probe qubit, i.e.

$\begin{split}\left(-\frac{1}{2}\omega Z_{0}+H_{\rm
Pauli}\right)\ket{a,n}&=E_{a,n}\ket{a,n}\\\
&=\left(\frac{(-1)^{a+1}}{2}\omega+E_{n}\right)\ket{a,n}.\end{split}$ (8)

Let $\bar{E}_{an}$ be the energy of the entire resonance Hamiltonian and
$\ket{\overline{an}}$ be the eigenvectors so that

$\begin{split}\left(-\frac{1}{2}\omega Z_{0}+H_{\rm
Pauli}+cX_{0}X_{1}\right)\ket{\overline{an}}=\bar{E}_{an}\ket{\overline{an}}\end{split}$
(9)

To first order in $c$ we have,

$\begin{split}&\bar{E}_{an}\approx
E_{an}+c\bra{a,n}X_{0}X_{1}\ket{a,n}=E_{an}\\\
&\ket{\overline{an}}\approx\ket{an}+c\sum_{bm\neq
an}\ket{bm}\frac{\bra{bm}X_{0}X_{1}\ket{an}}{E_{an}-E_{bm}}.\end{split}$ (10)

The initial state of the algorithm has the probe qubit in the ground state and
the system qubits in an arbitrary state. We are interested in the probability
that the auxiliary qubit will flip regardless of the state of the system after
time evolution. We can write the entire algorithm in one line,

$\braket{Z_{0}}(\omega)=1-2\sum_{m}\left|\sum_{n}\alpha_{n}\bra{1m}U_{\rm
res}(\omega)^{N_{t}}\ket{0n}\right|^{2}$ (11)

where the parameters $\alpha_{n}$ are arbitrary reflecting that we start in an
arbitrary state of the system. To evaluate this probability, we need to
evaluate

$\begin{split}\bra{1m}\left(-\frac{1}{2}\omega
Z_{0}+H_{S}+cX_{0}X_{1}\right)^{l}\ket{0n}=\\\
\sum_{ak}\braket{1m}{\overline{ak}}\braket{\overline{ak}}{0n}\bar{E}_{ak}^{l}\end{split}$
(12)

where $l$ is any non-negative integer. To first order, we have

$\begin{split}\braket{1m}{\overline{ak}}=\delta_{a1}\delta_{km}+c\frac{\chi_{mk}}{E_{0m}-E_{1k}}\delta_{a0}\end{split}$
(13)
$\begin{split}\braket{\overline{ak}}{0n}=\delta_{a0}\delta_{kn}+c\frac{\chi_{kn}}{E_{1k}-E_{0n}}\delta_{a1}\end{split}$
(14)

where $\chi_{mn}=\bra{m}X_{1}\ket{n}$. We multiply these together and get,

$\begin{split}\bra{1m}\left(-\frac{1}{2}\omega
Z_{0}+H_{S}+cX_{0}X_{1}\right)^{l}\ket{0n}=\\\
c\chi_{mn}\left(\frac{E_{1m}^{l}}{E_{1m}-E_{0n}}+\frac{E_{0n}^{l}}{E_{0m}-E_{1n}}\right).\end{split}$
(15)

Now we can go back to the probability

$\begin{split}&\braket{Z_{0}}(\omega)=\\\
&1-2c^{2}\left|\sum_{mn}\alpha_{n}\chi_{mn}\left(\frac{e^{iE_{1m}t}}{E_{m}-E_{n}+\omega}+\frac{e^{iE_{0n}t}}{E_{m}-E_{n}-\omega}\right)\right|^{2}\end{split}$
(16)

We see that the biggest contributions to the probability are the energies
which are on resonance with the probe qubit. The pole at
$\omega=\pm(E_{m}-E_{n})$ is unphysical since the perturbation is only valid
for $c<E_{m}-E_{n}-\omega$ for all $m$ and $n$.

To analyze the behavior near the resonance condition
$\omega=\pm(E_{1}-E_{0})$, we can ignore the other energy levels since they do
not contribute as strongly. In this case we are left with a two level system.
We can fully analytically analyze our algorithm for this two level system
which will tell us generally the behaviour of the algorithm when $\omega$ is
near an energy transition. Let us write the two level system Hamiltonian as,

$H_{0}=\frac{1}{2}dZ_{1}$ (17)

where $d=E_{1}-E_{0}$. The full Hamiltonian is,

$H=\frac{1}{2}\omega Z_{0}+\frac{1}{2}dZ_{1}+cX_{0}X_{1}$ (18)

This has eigenvalues

$\begin{split}E_{a}^{\pm}=\pm\frac{1}{2}\sqrt{4c^{2}+(d-\omega)^{2}}\\\
E_{b}^{\pm}=\pm\frac{1}{2}\sqrt{4c^{2}+(d+\omega)^{2}}\end{split}$ (19)

and eigenvectors

$\begin{split}&|a\pm\rangle=\left(\frac{d-\omega+2E_{a}^{\pm}}{2N_{a}^{\pm}}|01\rangle-\frac{c}{N_{a}^{\pm}}|10\rangle\right)\\\
&|b\pm\rangle=\left(\frac{d+\omega+2E_{b}^{\pm}}{2N_{b}^{\pm}}|00\rangle+\frac{c}{N_{b}^{\pm}}|11\rangle\right)\end{split}$
(20)

where

$\begin{split}&N_{a}^{\pm}=\frac{1}{\sqrt{2}}\sqrt{(d-\omega)^{2}+4c^{2}-2(d-\omega)E_{a}^{\pm}}\\\
&N_{b}^{\pm}=\frac{1}{\sqrt{2}}\sqrt{(d+\omega)^{2}+4c^{2}+2(d+\omega)E_{b}^{\pm}}\end{split}$
(21)

Thus, we can write time evolution as,

$\begin{split}e^{-iHt}&=e^{-iE_{a}^{-}}|a-\rangle\langle
a-|+e^{-iE_{a}^{+}}|a+\rangle\langle a+|\\\ &+e^{-iE_{b}^{-}}|b-\rangle\langle
b-|+e^{-iE_{b}^{+}}|b+\rangle\langle b+|\end{split}$ (22)

Using this form of the time evolution operator it is straight forward
(although tedious) to derive the bracket,

$\langle
1m|e^{iH(\omega)t}|0n\rangle=A(\omega)\delta_{n1}\delta_{m0}+B(\omega)\delta_{n0}\delta_{m1}$
(23)

where

$\begin{split}A(\omega,d)=i\frac{2c}{\sqrt{4c^{2}+(d-\omega)^{2}}}\sin(E_{a}t)\\\
B(\omega,d)=i\frac{2c}{\sqrt{4c^{2}+(d+\omega)^{2}}}\sin(E_{b}t)\end{split}$
(24)

Once again, we see that the probability is maximized when $\omega=\pm d$.

Figure 6: Fit of the ibm_lagos data to the functional form. The orange curves
are the experimental data, the blue curves are the fit, and the green dashed
lines show the location of the transitions according to the fit.

We use the functional form to assist in finding the minima in the data. Figure
6 shows the functional form fitted to the data from ibm_lagos. We begin the
fitting by smoothing the data (averaging over the nearest four data points)
and then finding the zeros in the derivative. We expect two energy transitions
within the data range, specifically the [0,1] transition and the [2,3]
transition which are both small. If we find two zeros (left panel in Fig. 6))
then we take the ten nearest points to each minimum found by the derivative
and separately fit them to $|A(\omega,d)|^{2}$. If we find only one minimum
(right panel in Fig. 6) then we take the 20 nearest points and fit it to
$|A(\omega,d1)+A(\omega,d2)|^{2}$. In other words, we expect that the two
transitions are close enough that they form a single dip.

##  III Zero energy resonance dip

In Fig. 2 of the main text, we see a dip in $\braket{Z_{0}}(\omega)$ at
$\omega=0$. However, this peak is absent in the Kitaev chain data. To
understand the origin of the zero-peak let us return to the function form of
the probe qubit orientation as derived in appendix II

$\begin{split}&\braket{Z_{0}}(\omega)=\\\
&1-2c^{2}\left|\sum_{mn}\alpha_{n}\chi_{mn}\left(\frac{e^{iE_{1m}t}}{E_{m}-E_{n}+\omega}+\frac{e^{iE_{0n}t}}{E_{m}-E_{n}-\omega}\right)\right|^{2}.\end{split}$
(25)

We see that there is a dip in $\braket{Z_{0}}(\omega)$ when
$\omega=E_{m}-E_{n}$. Since the sum is over all $m$ and $n$ we should expect a
zero energy dip for $m=n$ unless $\chi_{nn}=0$. Recall that

$\chi_{mn}=\bra{m}X_{1}\ket{n}$ (26)

which may not be zero depending on the eigenvectors $\ket{n}$ of the
Hamiltonian. If, however, the Hamiltonian describes fermions and we are using
the Jordan-Wigner transformation as in the main text then $X_{1}$ creates or
destroys a single fermion. Since fermion Hamiltonians conserve parity, they
always have $\chi_{nn}=0$. Therefore, we should not expect a zero dip for
fermion Hamiltonians. If, however, there is some small error in term which
couples the system to the probe when implemented on the quantum hardware then
a zero dip can show up. For this work we were careful to find qubits that
which did show a zero peak for the fermion model.

##  IV Optimization

Figure 7: Expectation values of the $Z$-component of the probe qubit as a
function of $\omega$ and $c$. The other parameters are set to $x=1.5$,
$y=0.4$, $z=-0.3$, $t=5.0$, $dt=1.2$. The dashed red lines show the location
of the energy transitions. The lines are labeled on the right of the figure by
the two states involved in the transition, e.g. $[0,1]$ is the transition from
state $0$ to state $1$.  Figure 8: Expectation values of the $Z$-component of
the probe qubit as a function of $\omega$ and $t$ for $dt=0.5$ (left) and
$dt=1.0$ (right). The other parameters are set to $x=1.5$, $y=0.4$, $z=-0.3$,
$t=0.4$. The dashed red lines show the location of the energy transitions. The
lines are labeled on the right of the figure by the two states involved in the
transition, e.g. $[0,1]$ is the transition from state $0$ to state $1$.

The coupling strength $c$ between the probe qubit and the system, the total
evolution time $t$, and the time step $dt$ all need to be optimized for the
algorithm. In Fig. 7 we see the orientation of the probe qubit as a function
of both $\omega$ and $c$. We want to focus on the ground to first excited
state transition labeled $[0,1]$ in the figure. We see that when $c>0.5$ the
dip corresponding to the $[0,1]$ transition is shifted due to the large
interaction between the probe and the system, and in fact the reverse
transition $[1,0]$ and higher-order transitions $[2,3]$ and $[3,2]$ also
appear (yet also shifted) at these strong coupling values. However, for
$c<0.1$ the dip is washed out. We choose $c=0.3$ where the dip corresponding
to the correct transition ($[0,1]$) is apparent and at the correct value, and
no other nearby transitions are observed.

Once the optimum $c$ value was found we optimized time step $dt$ and total
time $t$ in a similar way. Then we re-optimized $c$ at the new values of $t$
and $dt$ and found that the optimum $c$ value had not changed. Figure 8 shows
the $Z$-expectation of the probe qubit as a function of $\omega$ and $t$ for
two values of $dt$. We see that if $t$ is too short the dip is washed out just
like for small $c$. Additionally, if $dt$ is too small then we run into errors
due to the increased number of gates that are applied to the quantum register.
We want the smallest $dt$ we can manage so that our Suzuki-Trotter
decomposition is as accurate as possible.

Informed by these optimization experiments and others, we chose the parameter
set $t=5.0$, $dt=0.7$, $c=0.3$ for this work. More data was analyzed than
those presented here, and we have made our full data set accessible at
www.github.com/qiskit-research/mzm-phase-boundary.

##  V Applying the chemical potential symmetry

Figure 9: The raw data of the energy gap as a function of $m$ and $y$. The
parameter values match those in Fig. 4 of the main text. The dashed blue curve
shows the location of the phase transition calculated from Eq. (6) of the main
text while the green curve is a best fit of the data to Eq. (6) using $z$ as
the fit parameter.

We know that the energy spectrum is symmetric about $m\rightarrow-m$.
Therefore, if we know the zero-crossing phase boundary for $+m$ we can infer
the phase boundary for $-m$ and visa-versa. Figure 9 shows the raw data for
the energy gap as a function of $m$ and $y$. In the main text we selected the
smallest gap between positive and negative $m$ to be the gap for both. In this
way, we filter out the data in which the [0,1] transition has been washed out
by the [2,3] transition. Sometimes the $[2,3]$ transition has a much larger
response than the $[0,1]$ transition. When this happens the $[2,3]$ transition
can wash out the $[0,1]$ transition so that the fitting technique fails to
provide an accurate estimate of the $[0,1]$ transition. Thus, we get more
accurate data by selecting between the transitions at positive and negative
$m$. We find that the zero crossing of the $[0,1]$ transition is easily
extractable from the data for negative $m$ in the range of $0.7\leq y\leq 1.1$
while for $1.1<y$ the positive $m$ values show the $[0,1]$ zero-crossing
clearly.

##  VI Approaching the $n\rightarrow\infty$ limit in the Kitaev chain

We can approach the $L\rightarrow\infty$ limit by increasing the number of
sites. In Fig. 10 we plot maps of the energy gap for different values of $L$.
The parameters are set to $x=1.5$, $z=\bar{m}=0.4$. For $L>2$, we see that the
gap opens and closes multiple times. The larger $L$ the more the gap closes
and the smaller are the peaks between the lines where the gap closes. In this
way, one could approach the case where the gap stays closed.

Figure 10: Energy gap for increasing chain length. The plots are generated
from a classical simulation. In panel (a) we plot the (absolute) energy of the
lowest even and lowest odd parity states versus $m$ for two sites and for six
sites at $y=1.3$ ($m$ and $y$ are defined in Eq. (6)). In the other panels,
the energy gap between the lowest even and lowest odd parity state is plotted
against the parameters $m$ and $y$. The points where $|\Delta E|$ goes to zero
corresponds to the crossing of the two lowest energy levels. Each plot has a
different number of sites: (b) three sites, (c) four sites, (d) five sites,
(e) six sites, the dotted yellow line shows the location ($y=1.3$) of the six
site energy cuts in panel (a). For a plot of the simulated 2-site model, see
Fig. 4a.

##  VII Two site party flipping boundary

In the main text we stated that the parity of the ground state changes along
the boundary

$m=\sqrt{z^{2}+z(x+y)+xy}$ (27)

for the two site Kitaev chain. This comes simply from diagonalizing the two
site Hamiltonian

$H=mZ_{1}+mZ_{2}+xX_{1}X_{2}+yY_{1}Y_{2}+zZ_{1}Z_{2}$ (28)

which has eigenvalues

$\begin{split}&E_{a}^{\pm}=-z\pm(x+y)\\\
&E_{b}^{\pm}=z\pm\sqrt{4m^{2}+(x-y)^{2}}\end{split}$ (29)

The parameter regime of interest, $E_{a}^{-}$ and $E_{b}^{-}$ are the two
lowest energy levels. We want to know when the ground state switches from
$E_{a}^{-}$ to $E_{b}^{-}$. Setting $E_{a}^{-}=E_{b}^{-}$ and solving for $m$
gives the desired equation.

##  VIII Three site model data

Figure 11: Expectation of the $Z$-component of the probe qubit as a function
of $\omega$ for a three-site Kitaev chain. Left: for $m=0.5$. Right: for
$m=-0.5$. The other parameters are set to $x=1.5$, $y=1.1$, $z=0.4$, $c=0.3$,
$dt=0.7$, $t=5$. The green dashed lines show the energy transitions which are
labeled at the top/bottom of the plot. Figure 12:  The $[5,4]$ energy
transition as a function of $y$ and $m$. Left: from directly solving for the
eigenvalues of the Hamiltonian. Right: from the spectroscopic algorithm ran on
Lagos. The other paramaters are set to $x=1.5$, $z=0.4$, $c=0.3$, $dt=0.7$,
$t=5$.

We ran the spectroscopic algorithm for a three site Kitaev chain. While the
data is very accurate at certain parameter values, we where not able to map
out the full phase diagram. One issue is that the number of small-energy
transitions increases causing the width of the dips in
$\braket{Z_{0}}(\omega)$ to often absorb several of the weaker transitions. In
Fig. 11 we plot $\braket{Z_{0}}(\omega)$ for $m=0.5$ (left) and $m=-0.5$
(right). While the energy transitions are symmetric in $m$, the algorithm
often favors different transitions for positive and negative $m$. In this
case, the $[5,4]$ transition is favored for positive $m$ while $[3,2]$ is
favored for negative $m$. While this asymmetry in the favored transition is
present for the two site case as well, there we where able to resolve the
secondary transition for many parameter values. For the three site model, the
secondary transitions are very difficult to resolve if they are present at
all.

Still, we are able to trace out the zero crossings for certain energy
transitions for a range of parameters. Figure 12 shows the zero-crossing for
the $[5,4]$ transition in a small range of $m$ and $y$. While the transition
is shifted just like in the two site case, the main features are present. The
data could be improved by narrowing the dips in $\braket{Z_{0}}(\omega)$. A
possible solution for narrowing the dips would be to increase $t$, which for
this current work was limited by qubit coherence and classical waveform
generation bandwidth for ibm_lagos. It is expected that future generations of
superconducting quantum backends will allow for such exploration.


###### Abstract

Classification of $n-$th $(n\geq 2)$ order linear ODEs is considered. The
equation reduced to Laguerre Forsyth form by a point transformation then, the
other calculations would have done on this form. This method is due to V.A.
Yumaguzhin.
Keywords: Linear ODE, symmetry, Lie algebra, projective transformations.
© 200x Published by Islamic AZAD University-Karaj Branch.

_Mathematical Sciences_  _Vol. 1, No. 1,2 (2007) 01-12_

Classification of $n-$th order linear ODEs up to projective transformations

Mehdi Nadjafikhaha,111Corresponding Author. E-mail Address:
<EMAIL_ADDRESS>, Seyed Reza Hejazib

aFact. of Math., Dept. of Pure Math., Iran University of Science and
Technology, Narmak, Tehran, I.R. Iran.
bSame address.

## Introduction

The local classification of linear ODEs up to projective transformations is
obtained in this article. For $n\leq 2$, it is well known that any $n-$th
order linear ODE can be transformed locally to the form $y^{(n)}=0$ by a point
transformation. For $n\geq 0$, this statement is incorrect: there is finite
number of different equivalence classes of linear ODEs.

First this problem was posed by classics of the 15 century E. Laguerre, G.H.
Halphen and others. They obtain results concerning classification of third and
fourth orders linear ODE. Here, this problem is solved for $n\geq 0$ in a
neighborhood of regular germs.

Consider a general $n-$th order ODE which is solved by the higher order
derivative

$\displaystyle y^{(n)}=\sum_{i=1}^{n}a_{n-i}(x)y^{(n-i)},$ (1)

where $y(x)$ is a smooth function of $x$.

Lie shows that the point symmetry group of a second ordinary linear
differential equation has dimension at most eight, conversely the equation
admits an eight-dimensional symmetry group if and only if it can be mapped, by
a point transformation, to the linear equation $y^{\prime\prime}=0$. Thus, the
main result is any linear second ordinary differential equation can mapped to
the equation $y^{\prime\prime}=0$. So, the condition of second ordinary linear
differential equation is specified.

A same result shows that for $n\geq 3$, any linear ODE admits at most an
$(n+4)-$dimensional symmetry group of point transformation, therefore, the
symmetry group is $(n+4)-$dimensional if and only if the equation is
equivalent to the linear equation $y^{(n)}=0$. In continuation we will work on
the general form of linear ODE in the form of (1) once $n\geq 3$.

##  1 Laguerre-Forsyth form

The classification of linear differential equations is a special case of the
general problem of classifying differential operators, which has a variety of
important applications. Consider an $n-$th order ordinary differential
operator corresponding to (1)

$\displaystyle{\cal
D}=a_{n}(x)D_{x}^{n}+a_{n-1}(x)D_{x}^{n-1}+\cdots+a_{1}(x)D_{x}+a_{0}.$ (2)

The aim is finding out when two operators, or two linear ODE, of type (2), can
be mapped to each other by a suitable change of variables. To preserve
linearity, we restrict to those of the form

$\displaystyle\bar{x}=\varphi(x),\qquad\bar{y}=\psi(x)y,$ (3)

the chain rule action shows that $D_{\bar{x}}=(\varphi(x))^{-1}D_{x}$, and
with a rescaling of the dependent variable by
$\displaystyle\psi(x)=e^{\varphi(x)}$ we obtain the gauge factor. So, two
differential operator $\overline{\cal D}$ and $\cal D$ is called gauge
equivalent if they satisfy

$\displaystyle{\overline{\cal D}}=\psi\cdot{\cal
D}\cdot\frac{1}{\psi},\qquad\bar{x}=\varphi(x).$ (4)

A straightforward calculation shows that the change of variables (3) is given
by

$\displaystyle\bar{x}=\varphi(x)=\int\frac{dx}{\sqrt[n]{|a_{n}(x)|}},\;\;\;\;\psi(x)=|a_{n}(x)|^{\frac{1-n}{2n}}\exp\Big{\\{}\int^{x}\frac{a_{n-1}(y)}{na_{n}(y)}dy\Big{\\}},$

thus (1) is gauge equivalent to an operator of the form

$\displaystyle{\cal D}=\pm D_{x}^{n}+a_{n-2}(x)D_{x}^{n-2}+\cdots+a_{0}(x).$
(5)

If $\rho(x)$ be a nonvanishing smooth function, two differential operator
$\overline{\cal D}$ and $\cal D$ is called projective equivalence if they
satisfy

$\displaystyle{\overline{\cal D}}=\rho\cdot\psi\cdot{\cal
D}\cdot\frac{1}{\psi},\qquad\bar{x}=\varphi(x).$ (6)

A nonsingular $n-$th order linear operator of type (5) is projectively
equivalent to one in Laguerre-Forsyth form

$\displaystyle{\cal D}=D_{x}^{n}+a_{n-3}(x)D_{x}^{n-3}+\cdots+a_{0}(x),$ (7)

with change of variable (6) in the form of

$\displaystyle\bar{x}=\varphi(x),\qquad\bar{y}=\varphi_{x}^{\frac{n-1}{2}}y,\qquad\rho=\varphi_{x}^{-n},$

where $\varphi(x)$ is a solution of the Schwarzian equation

$\displaystyle\frac{n(n^{2}-1)}{12}\frac{\varphi_{x}\varphi_{xxx}-\frac{3}{2}\varphi_{xx}^{2}}{\varphi_{x}^{2}}=a_{n-2}(x).$

##  2 Classification of linear ODEs of Laguerre-Forsyth form

A useful theorem help us to reduce the classification of ODEs up to a special
transformation.

######  Theorem 2.1

Let $\Delta_{1}$ and $\Delta_{2}$ be ODEs of the form (7). If there is a point
transformation that takes $\Delta_{1}$ to $\Delta_{2}$, that is

$\displaystyle
f(x)=\frac{ax+b}{cx+d},\quad\hat{f}(x,y)=|f^{\prime}|^{\frac{n-1}{2}}\cdot
y,\quad a,b,c,d\in\bf{R}.$ (8)

A transformation $(f,\hat{f})$ of the form (8) is generated by a projective
transformation $f$ on $\bf R$. The isomorphisms $f\rightarrow(f,\hat{f})$
makes a group of point transformations in the form of (8). Consider these
projective transformations in a group $G$ and denoted by all projective
transformations of $\bf R$,i.e.,

$\displaystyle G=\Big{\\{}f(x)=\frac{ax+b}{cx+d}\Big{|}a,b,c,d\in{\bf
R}\;\mbox{and}\;ad\neq bc\Big{\\}}.$

It is easy to check that $G$ has two connected component $G_{1}=\\{f\in
G|f^{\prime}>0\\}$ and $G_{2}=\\{f\in G|f^{\prime}<0\\}$, thus, $G=G_{1}\cup
G_{2}$.

###  2.1 Bundles of Laguerre-Forsyth form

Consider $x$ as a coordinate on $\bf R$ and $a_{n-3},a_{n-2},...,a_{0}$
coordinates on ${\bf R}^{n-2}$. Then, we can construct a fiber bundle
corresponding to (7) in the form of

$\displaystyle p:{\bf R}\times{\bf R}^{n-2}\rightarrow{\bf R}.$ (9)

Any ODE of type (7) identifies with
$\Delta=\\{p_{n}=a_{n-3}(x)p_{n-3}+\cdots+a_{0}(x)p_{0}\\}$ is a section of
(9) denoted by $S_{\Delta}:x\rightarrow(x,a_{n-3}(x),...,a_{0}(x))$, where the
identification $\Delta\rightarrow S_{\Delta}$ is a bijection.

Let
$\Delta_{2}=\\{\tilde{p}_{n}=\tilde{a}_{n-3}(\tilde{x})\tilde{p}_{n-3}+\cdots+\tilde{a}_{0}(\tilde{x})\tilde{p}_{0}\\}$
be an ODE of the form (7). Subjecting $\Delta_{2}$ to an transformation (8),
the, we obtain linear ODE
$\Delta_{1}=\\{p_{n}=a_{n-3}(x)p_{n-3}+\cdots+a_{0}(x)p_{0}\\}$. The
coefficients $\Delta_{2}$ are expressed in terms of coefficients of
$\Delta_{1}$ and projective transformation $f^{-1}$ by the equation

$\displaystyle\tilde{a}_{n-i}=F_{n-i}\Bigg{(}a_{n-3},...,a_{n-i};\frac{df^{-1}}{d\tilde{x}},...,\frac{d^{i+1}f^{-1}}{d\tilde{x}^{i+1}}\Bigg{)},\quad
i=3,4,...,n.$ (10)

The equation (10) is a lifting of a projective transformation $f$ to
diffeomorphism $\bar{f}:{\bf R}\times{\bf R}^{n-2}\rightarrow{\bf R}\times{\bf
R}^{n-2}$ such that $p\circ\bar{f}=f\circ p$.

For any $f\in G$, a transformation of sections of $p$ defined by the formula

$\displaystyle S\rightarrow f(S)=\bar{f}\circ S\circ f^{-1},$

then, equation (10) can be represented as $S_{\Delta_{2}}=f(S_{\Delta_{1}})$.

######  Lemma 2.2

Consider two equation of the form (7). Then a transformation $(f,\hat{f})$ of
the form (8) maps $\Delta_{1}$ to $\Delta_{2}$ if and only if
$S_{\Delta_{2}}=f(S_{\Delta_{1}})$.

The main result of the lemma (2.2) is the classification of ODEs of the form
(7) up to transformation (8) reduces to classification of germs of sections of
$p$ up to projective transformation on $\bf R$.

###  2.2 Classification of regular germs

Let $S$ be a section of $p$ and $a$ be a point in domain of $S$. Denoted by
$\\{S\\}_{a}$ the germ of $p$ at $a$. Let $\\{S\\}_{a_{1}}$ and
$\\{S\\}_{a_{2}}$ be germs of sections $S_{1}$ and $S_{2}$ respectively. We
say that $\\{S\\}_{a_{1}}$ and $\\{S\\}_{a_{2}}$ are $G_{+}$-equivalent if
there exist $f\in G_{+}$ such that
$\\{f(S_{1})\\}_{f(a_{1})}=\\{S\\}_{a_{2}}$. A germ $\\{S\\}_{a}$ is regular
of class i if there exist a neighborhood $\cal O$ of $a$ and subbundle $E_{i}$
such that Im$S|_{\cal O}\subset E_{i}$. If $\\{S\\}_{a}$ is a regular germ of
class $i\geq 0$, then in a neighborhood of $a$ we have
$S(x)=(x,0,...,0,a_{i}(x),...,a_{0}(x))$. In the rest of the paper we will
often denoted $\\{S\\}_{a}$ by $\\{a_{i},...,a_{0}\\}_{a}$. If $\\{S\\}_{a}$
is a regular germ, then $a$ is a regular point of $S$.

######  Definition 2.3

Let $S$ be a section of $p$ and v be a vector field of the Lie algebra of
group $G$, if $\theta_{t}$ be the flow of v, we say v is a projective symmetry
of S if one of the following statements satisfied:

  * 1)

$\theta_{t}(S)=\overline{\theta_{t}}\circ S\circ\theta_{t}^{-1}=S$,

  * 2)

$\displaystyle\frac{d}{dt}\theta_{t}(S)\Big{|}_{t=0}=0.$

Denote by ${\cal P}(S)$ the Lie algebra of all projective symmetries of $S$.

Let $\Upsilon$ be the set of all regular germs at $0\in\bf R$ of sections of
$p$. Define

$\displaystyle\Upsilon_{i}=\Big{\\{}\\{S\\}_{a}|\mbox{dim}{{\cal
P}(S)}=i\Big{\\}},\quad i=0,1,3,$

and denote $\Upsilon=\Upsilon_{0}\cup\Upsilon_{1}\cup\Upsilon_{3}$. If $G_{0}$
be the isotropic subgroup of $G$ in 0, then, $\Upsilon_{i}$’s are
$G_{0}$-invariant.

Define $\Upsilon_{r,i}\subset\Upsilon_{r}$ be the subset of all regular germs
of class $i$. It follows from the invariance of subbundle $E_{i}$’s under
$G_{0}$, $\Upsilon_{r,i}$ is $G_{0}-$invariant. Consequently we have

$\displaystyle\Upsilon_{r}=\bigcup_{i=0}^{n-3}\Upsilon_{r,i},$

where this union is separated invariant subsets.

Let ${\bf R}_{+}$ and ${\bf R}_{-}$ be the set off positive and negative real
numbers respectively. If $\ell_{r,i}:\Upsilon_{r,i}\rightarrow({\bf
R}\backslash\\{0\\})\times{\bf R}$ be a map by the formula
$\\{a_{i},...,a_{0}\\}\mapsto(a_{i}(0),a_{i}^{\prime}(0))$ and

$\displaystyle G_{0+}\times\Upsilon_{r,i}\rightarrow\Upsilon_{r,i}$ (11)
$\displaystyle(f,\\{S\\}_{0})\mapsto\\{f(S)\\}_{0},$

be the action of $G_{0+}$ on $\Upsilon_{r,i}$ then,

######  Lemma 2.4

The map $\ell_{r,i}|_{\Theta}$ is a bijection from the orbit $\Theta$ of the
action (11) either to $({\bf R}_{+})\times{\bf R}$ or to $({\bf
R}_{-})\times{\bf R}$.

Let $\Omega_{r,i}^{+}=\ell_{r,i}^{-1}((1,0))$ and
$\Omega_{r,i}^{-}=\ell_{r,i}^{-1}((-1,0))$. Denote by $\Gamma_{r,i}$ the
subset of $\Omega_{r,i}^{+}\cup\Omega_{r,i}^{-}$ defined in the following way:

  * 1)

$\Gamma_{r,0}=\Omega_{r,0}^{+}\cup\Omega_{r,0}^{-}$ for $i=0$,

  * 2)

if $i>0$, then, $\Omega_{r,i}$ consists of all germs $\\{a_{i},...,a_{0}\\}$
from $\Omega_{r,i}^{+}\cup\Omega_{r,i}^{-}$ satisfying one of the following
conditions:

    * i)

$a_{i-j}=0$ for all odd numbers $j$ with $1\leq j\leq i$,

    * ii)

there exist an odd number $r$ with $1\leq r\leq i$ such that $a_{i-r}(0)>0$
and if $r>1$, then $a_{i-j}(0)=0$ for all odd numbers $j$ with $1\leq j<r$.

###  2.3 Classification of regular germs from the family
$\displaystyle\bf\Omega_{r,i}$

Let $\mu\in G_{-}$ defined by $\mu(x)=-x$ for all $x\in\bf R$, then, due to
lemma (2.4) and attentive to $\mu(\Omega_{r,i}^{-})=\Omega_{r,i}^{+}$ we have:

######  Theorem 2.5

  * 1)

The set $\Omega_{r,i}^{+}\cup\Omega_{r,i}^{-}$ is a family of all germs from
$\Upsilon_{r,i}$ nonequivalent with respect to $G_{0+}$.

  * 2)

If $n-i$ is odd, then $\Omega_{r,i}^{+}$ is a family of all germs from
$\Upsilon_{r,i}$ nonequivalent with respect to $G_{0}$.

  * 3)

If $n-i$ is even, $\Gamma_{r,i}$ is a family of all germs from
$\Upsilon_{r,i}$ nonequivalent with respect to $G_{0}$.

an important corollary concludes this section as follows:

######  Corollary 2.6

Classification of regular germs of sections of (7) is:

  * 1)

The family of germs of the form

$\displaystyle\\{\pm 1+a(x)x^{2},a_{i-1}(x),...,a_{0}\\}_{0}$

is a family of all regular germs of class $i$ nonequivalent with respect to
$G_{0+}.$

  * 2)

If $n-i$ is odd, then the family of germs of the form

$\displaystyle\\{1+a(x)x^{2},a_{i-1}(x),...,a_{0}\\}_{0}$

is a family of all regular germs of class $i$ nonequivalent with respect to
$G_{0}.$

  * 3)

If $n-i$ is even, then the family of germs of the form

$\displaystyle\\{\pm 1+a(x)x^{2},a_{i-1}(x),...,a_{0}\\}_{0},$

satisfying one of the following conditions:

    * a)

$a_{i-j}(0)=0$ for all odd numbers $j$ with $1\leq j\leq i$,

    * b)

there exist an odd number $r$ with $1\leq r\leq i$ such that $a_{i-r}(0)>0$
and if $r>1$, then $a_{i-j}(0)=0$ for all odd number $j$ with $1\leq j\leq r$,

is the family of germs of class $i$ nonequivalent with respect to $G_{0}.$

##  3 Conclusion

This article was a qualification of classification of linear ODEs due to V.A.
Yumaguzhin. First we transform the general form of ODEs to Laguerre-Forsyth
form, then by a suitable change of variable up to projective transformation we
reduce this classification to classification of the sections of bundles, next
by construction germs and specially regular germs of this sections near
identity, the classification reduced to classifying of regular germs by
providing some invariant subsets of the bundles.

## References

  * [1] Lee, John,M., Introduction to Smooth Manifolds, Springer Verlage, New York, 2002.
  * [2] Olver P.J., Equivalence, Invariant and Symmetry, Cambridge University Press, Cambridge 1995.
  * [3] Olver P.J., Applications of Lie Groups to Differential equations, Seconed Edition, GTM, Vol. 107, Springer Verlage, New York, 1993.
  * [4] Olver P.J., Differential Invariant and Differential Invariant Equations, University of Min-nesuta, 1994.
  * [5] Olver P.J., Differential Invariants: Algebraic and Geometric Structure in Differential Equa-tions, P.H.M. Kersten and I.S. Krasil’shchik, eds., Proceeding, University of Twente, 1993, to appear.
  * [6] Ovsiannikov, L.V., Group Analysis of Differential Equations, Academic press, New York, 1982.
  * [7] Yumaguzhin V.A., Contact classification of linear differential equations. I., Program System Institute, M. Botik, Preslavl-Zalessky, 152020, Russia.
  * [8] Yumaguzhin V.A., Contact classification of linear differential equations. II., Program System Institute, M. Botik, Preslavl-Zalessky, 152020, Russia.
  * [9] Yumaguzhin V.A., Point transformation and classification of 3rd-order linear ODEs, Russian journal of Mathematical Physics, 4 (1996) No. 3, 403-410.
  * [10] Yumaguzhin V.A., Classification of 3rd-order linear ODEs up to equivalence, Journal of Differential Geometry and its applications, 6 (1996) No. 4, 343-350.