Simulating spectroscopy experiments with a superconducting quantum computer

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 [1] 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 [2][3][4][5][6][7], even starting to surpass the accuracy of state-of-the-art classical approximations [8]. While quantum phase estimation (QPE) [9] offers an exponential speed-up in determining Hamiltonian spectra, it requires a depth of circuit unfeasible for near-term noisy quantum hardware. Shortdepth 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 [10-14, 16? ]. 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.
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 [9,18], this technique does not require controlled-coherent evolution of the register by auxiliary qubits. Instead, the probe qubit evolving according to an energy (via ω) 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, sub- Figure 1. Quantum circuits for simulating the time evolution U Pauli = e −iH Pauli t : (a) a cartoon analogy of the spectroscopy eigensolver method to scanning tunneling spectroscopy, where a tip (qubit) at potential (energy) ω 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 λ as a dip in probe expectation Z0 . (c) The spectroscopic eigensolver algorithm for a probe qubit interacting with a single qubit of the simulation register (for Nt Trotter steps Ures, 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 (θ) interaction generated by echoed cross resonance [17]. stantially 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 ω, 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 [19].
Explicitly, the spectroscopic eigensolver technique finds the energy differences of a system by coupling a probe qubit to a Pauli Hamiltonian H Pauli suitably encoded to enforce the commutation relations of the system Hamiltonian H 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 where the probe qubit q 0 with energy ω 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 is an XXrotation on sites j and k and U Pauli = e −iH Pauli dt is the Trotter step of time evolution dt = t/N t for the system.
We start in a state |0, ψ 0 with the probe qubit in the ground and the system |ψ 0 is in an arbitrary state. For simplicity we will take |ψ 0 to be the nqubit ground state. After N t applications of U res , we measure the probe qubit in the Z-basis, Z 0 (ω) = 0, When c is small compared to the energy scales in H Pauli and when ω is on resonance with an energy transition of H Pauli , the probability of the probe qubit flipping will peak, resulting in a dip in Z 0 (ω).
The following experiments were conducted within the Qiskit [20] 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, ω 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 π rotation in the former case [2]. These scaled pulses are implemented as calibrations for R zx (θ) gates following an analysis of the circuit using a novel template optimization technique [21]. Each spectroscopic experiment is then built as an array of circuits, one for each value of ω. 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 Z 0 (ω). 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, 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 Z 0 (ω) for several values of a. Notice that there is a dip at ω = 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 Z 0 (ω). The Root-Mean-Square (RMS) of the variance of all three transi- where N d is the number of data points, w * i is a minimum and ∆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-Widthat-Half-Minimum (FWHM), which includes both hardware and Trotter error, FWHM = 1 is half way between its minimum and maximum value to the left of ω * i and ω R i is the same but to the right of ω * i . This simple test system stands to demonstrate how accurate our algorithm can be when the decoherence time of the qubits (∼ 100 − 150 µs) allows for a long evolution time t dt, where t corresponds to ∼ 15 µ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 [22][23][24][25][26][27]. 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 → ∞ limit, a parity flip occurs even in a two site model.

The model Hamiltonian of the interacting Kitaev chain with Coulomb interactions is
where L is the number of sites, c † i (c i ) is the creation (annihilation) operator on site i, µ is the chemical potential, g is the hopping rate, ∆ 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 [28,29].
Using the Jordan-Wigner encoding [18], we can express Eq. 4 as 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 + ∆, 2y = g − ∆, 4z = V , 4m = 2µ + V , and 4m = 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 [30].
In the limit n → ∞, 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 → ∞ 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 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 ω for three values of c. The energy difference between the ground and first excited state is determined by the ω at which Z 0 (ω) 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 |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 ω = 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.
Scans of Z 0 (ω) 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 ω = 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 ≈ 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 ≈ 1.72 for the [0,1] transition which is well within the width of the peak.
Taking several scans of ω 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 ∆z = 0.19 which is well within the average FWHM ≈ 2.35 found above. Using the best fit z we find that the shape of the boundary is fairly accurate with RMS ≈ 0.15 which is again well within the average FWHM. The error is largely due to the width of the Z 0 (ω) dips and the influence of nearby transitions. As seen in Fig. 3 the [0,1] dip is pulled to higher ω 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 [31,32]. 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 bound-   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 ∆z = 0.19. ary 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 ± 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. Hell  H Figure 5. Comparison of (a) QPE and (b) the spectroscipic eigensolver algorithms .

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 [9] -the backbone of proposed universal quantum simulation [33]. 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 ≤ k ≤ 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 [34]. 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 [35]. 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 [36]. 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.
The accuracy of the spectroscopic eigensolver on the Landau-Zener model is approximately equivalent to three bits of precision (∼ 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 [37].
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. [2]) with a floor of the X-error that forms the echo. These are reported in Table I. 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.

II. FUNCTIONAL FORM OF THE RESONANCE EIGENSOLVER
The algorithm works by applying unitary time evolution to an arbitrary state (in our simulation we choose |ψ 0 = |00 ) 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 ω in the resonance Hamiltonian, not the actual transition energy of transmon as in [38]) 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 res (ω) = e −i(−ωZ0/2+H Pauli +cX0Xi)dt (7) where H 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 |n be the energies and eigenvectors of the system Hamiltonian H Pauli |n = E n |n so that the unperturbed eigenvalues are E a,n = (−1) a+1 ω/2 + E n and the unperturbed eigenvectors are |a, n where a ∈ {0, 1} labels the state of the probe qubit, i.e.
LetĒ an be the energy of the entire resonance Hamiltonian and |an be the eigenvectors so that To first order in c we have, 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, where the parameters α n are arbitrary reflecting that we start in an arbitrary state of the system. To evaluate this probability, we need to evaluate where l is any non-negative integer. To first order, we have where χ mn = m| X 1 |n . We multiply these together and get, Now we can go back to the probability We see that the biggest contributions to the probability are the energies which are on resonance with the probe qubit. The pole at ω = ±(E m − E n ) is unphysical since the perturbation is only valid for c < E m − E n − ω for all m and n.
To analyze the behavior near the resonance condition ω = ±(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 ω is near an energy transition. Let us write the two level system Hamiltonian as, where d = E 1 − E 0 . The full Hamiltonian is, This has eigenvalues and eigenvectors where Thus, we can write time evolution as, Using this form of the time evolution operator it is straight forward (although tedious) to derive the bracket, where A(ω, d) = i 2c Once again, we see that the probability is maximized when ω = ±d. 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(ω, 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(ω, d1)+A(ω, 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 Z 0 (ω) at ω = 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 We see that there is a dip in Z 0 (ω) when ω = E m −E n . Since the sum is over all m and n we should expect a zero energy dip for m = n unless χ nn = 0. Recall that which may not be zero depending on the eigenvectors |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 χ 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.  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 ω 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 Once the optimum c value was found we optimized time step dt and total time t in a similar way. Then we reoptimized 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 ω 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
We know that the energy spectrum is symmetric about m → −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]   We can approach the L → ∞ 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 =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. which has eigenvalues 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
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 smallenergy transitions increases causing the width of the dips [3,2] [2,3] [5,4] [4,5]   in Z 0 (ω) to often absorb several of the weaker transitions. In Fig. 11 we plot Z 0 (ω) 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 Z 0 (ω). 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.