Multi-Level Variational Spectroscopy using a Programmable Quantum Simulator

Energy spectroscopy is a powerful tool with diverse applications across various disciplines. The advent of programmable digital quantum simulators opens new possibilities for conducting spectroscopy on various models using a single device. Variational quantum-classical algorithms have emerged as a promising approach for achieving such tasks on near-term quantum simulators, despite facing significant quantum and classical resource overheads. Here, we experimentally demonstrate multi-level variational spectroscopy for fundamental many-body Hamiltonians using a superconducting programmable digital quantum simulator. By exploiting symmetries, we effectively reduce circuit depth and optimization parameters allowing us to go beyond the ground state. Combined with the subspace search method, we achieve full spectroscopy for a 4-qubit Heisenberg spin chain, yielding an average deviation of 0.13 between experimental and theoretical energies, assuming unity coupling strength. Our method, when extended to 8-qubit Heisenberg and transverse-field Ising Hamiltonians, successfully determines the three lowest energy levels. In achieving the above, we introduce a circuit-agnostic waveform compilation method that enhances the robustness of our simulator against signal crosstalk. Our study highlights symmetry-assisted resource efficiency in variational quantum algorithms and lays the foundation for practical spectroscopy on near-term quantum simulators, with potential applications in quantum chemistry and condensed matter physics.

Digital quantum simulators, based on programmable quantum circuits, provide remarkable versatility and enable the study of arbitrary Hamiltonians [22].Variational quantum algorithms, such as the Variational Quantum Eigensolver (VQE) [23], present a promising avenue for exploring practical quantum applications using nearterm quantum simulators.By employing a hybrid approach that combines a quantum simulator with a classical optimizer, VQE simulations have successfully tackled ground state problems of many-body systems in quantum chemistry [24][25][26][27][28][29] and condensed matter physics [30,31].However, the considerable resource overheads associated with both quantum and classical components have restricted VQE spectroscopy to small system sizes, predominantly targeting only two eigenstates [32][33][34], with the exception of Ref. [28] which approximates four eigenenergies in a two-qubit system.In order to make quantum variational spectroscopy extensible, it is crucial to address resource overheads by simplifying the ansatz and reducing variational parameters.In fact, symmetry is a fundamental property that is found in every building block of our universe.The interaction between particles in strongly correlated many-body systems reveals several forms of symmetries.Theoretical proposals have suggested utilizing symmetries as a potential approach [35][36][37][38].Nonetheless, conclusive experimental evidence remains to be provided.
In this article, we present an experimental demonstration of multi-level energy spectroscopy for fundamental many-body Hamiltonians using VQE on a superconducting digital quantum simulator.By employing a mix of symmetry-preserving ansatzes, we significantly re-arXiv:2306.02110v1[quant-ph] 3 Jun 2023 duce circuit depth and the number of variational parameters while targeting multiple energy levels in conjunction with the Subspace Search VQE (SSVQE) method [39].We also introduce a circuit compiling method that conveniently corrects for both spatial and temporal signal crosstalk, improving the consistency of the gate performance in our device.We achieve full variational quantum spectroscopy for all 16 states of a 4-qubit Heisenberg spin chain, with the measured energies deviating from their theoretical values by 0.13 on average, assuming a unity coupling strength.We further extend our approach to 8-qubit Heisenberg and transverse-field Ising Hamiltonians, extracting the three lowest energy eigenstates and demonstrating the potential of our method for multi-level spectroscopy in many-body systems.Intriguingly, we find that in the cases where the cost function of the VQE algorithm contains several terms, the best convergence result is obtained using a hybrid optimizer that combines the gradient-free Nelder-Mead method with the gradient-based Adam method.
As illustrated in Fig. 1a, the VQE algorithm contains a quantum hardware unit that executes a quantum circuit consisting of an initialization part, a multi-layer parameterized circuit, and measurement.The initialization circuit prepares an initial state of a N -qubit system through unitary operation |Φ⟩=U I |0⟩ ⊗N .The parameterized circuit, however, is described by a unitary operation U ( ⃗ θ), where ⃗ θ = (θ 1 , θ 2 , . . ., θ l ) are variational parameters.The action of U ( ⃗ θ) on the initial state creates the output state |Ψ( ⃗ θ)⟩=U ( ⃗ θ)|Φ⟩ for which an observable can be efficiently measured.A proper cost function C( ⃗ θ) is calculated from the measured data which is then fed into a classical optimizer to be minimized through updating the circuit parameters ⃗ θ, a training process until an optimal set of parameters ⃗ θ * is obtained.
In comparison to typical VQE algorithms which find the ground state of a given Hamiltonian with the cost function being its expectation value [23-27, 30, 40], we employ a combination of several techniques to extend the spectrum and target more excited states.For a system with known symmetries, we can prepare an initial state with a specific set of symmetry numbers and then apply the corresponding symmetry-preserving circuit such that the search is restricted within a certain manifold.In this way, one can easily obtain the lowest energy eigenstate in that manifold -but an excited state to the whole system.For cases when the states are indistinguishable by symmetry or when it is hard to prepare the initial state, simulating the excited states requires either adding additional penalizing terms in the cost function [34,41] or simultaneous learning of multiple eigenstates through special training methods, such as SSVQE [39].In our experiment, we use a general form for the cost function where {|Φ i ⟩} P i=1 are a set of P orthogonal initial states and w i (real positive) are their weights with condition w 1 >w 2 > • • • >w P .Note that the output states |Ψ i ⟩ =U ( ⃗ θ) |Φ i ⟩ are naturally orthogonal due to the unitarity of U ( ⃗ θ).With proper constraints imposed on these weights, one obtains |E i ⟩ ≈|Ψ i ( ⃗ θ * )⟩ through classical minimization of the cost function, provided that the parameterized circuit exhibits adequate expressibility.
To understand how symmetry is integrated into the algorithm, let us first examine the well-known Heisenberg spin chain model that contains N spin-1/2 particles.The Hamiltonian is expressed by where J>0 is the coupling strength and σ i =(σ i x , σ i y , σ i z ) is the vector of Pauli operators at site i.The Heisenberg Hamiltonian in Eq. ( 2) supports several symmetries including: (i) the magnetization in every direction, i.e. [H, S α ]=0, where S α = or |1⟩ is the qubit state at site i.These symmetries imply that all the eigenstates |E j ⟩ have a specific total spin s, namely S 2 tot |E j ⟩ =s(s + 1) |E j ⟩, zmagnetization s z , given by S z |E j ⟩ =s z |E j ⟩ in which −s ≤ s z ≤ s, and the mirror symmetry m=±1, given by M |E j ⟩ =m |E j ⟩.
In general, each layer of our parameterized circuit contains two consecutive sets of two-qubit operations N (θ), acting on odd and even bonds, followed by a series of single-qubit rotations (see Fig. 1a).In the case of Heisenberg model, N (θ) acting on qubits i and i+1 is chosen to be An explicit design of N H (θ) as well as symmetrypreserving circuits are shown in Fig. 1b.Since [N H (θ), S 2 tot ]=0, to conserve the total spin s no singlequbit rotation is needed.To conserve s z , the singlequbit operations are chosen to be z-rotations, namely Z θ = e −iθσz/2 .Note that the s-conserving circuit readily conserves s z .The mirror symmetry can be easily preserved by choosing the gate parameters symmetrically with respect to the middle of the chain, and it may be used in combination with the other symmetries.
To diminish the effect of noise and improve the accuracy of the extracted energies, we employ the Zero-Noise N H (θ 3 ) N H,I (θ 1 ) N H,I (θ 2 ) Quantum simulator Classical computer N I (θ 3 ) Variational quantum algorithm with symmetry-preserving circuits.a, Schematic diagram of a variational quantum algorithm which is realized on a combination of a quantum simulator and a classical computer.The quantum circuit includes an initialization circuit UI , a parameterized symmetry-preserving circuit U ( ⃗ θ), and a measurement unit.The measurement results are fed into a classical optimizer for minimizing a cost function C( ⃗ θ) through iteratively updating the parameters ⃗ θ. b, Ansatz circuits for preserving different symmetries.In the top boxes are the two-qubit modular circuits used for the Heisenberg and transverse-field Ising problems.In the bottom are ansatz circuits for preserving z-polarization (sz), total spin (s), mirror symmetry (m), and zz-product (πz), respectively.c, Error mitigation using zero-noise extrapolation.UI and U † I are the initialization circuit and its reversal.Note that the N (θ/k) circuit is identical to the original N (θ) circuit except for single-qubit rotations according to the decomposition in (b).d, False-colour micrograph of the device.
Extrapolation (ZNE) method [42][43][44] for error mitigation.To perform ZNE, we rescale the circuit depths by a factor of k while maintaining the same unitary.As shown in Fig. 1c, for the parameterized circuit U ( ⃗ θ), every N (θ) operation is divided into k consecutive repetitions of N (θ/k).Such a treatment allows the noise scaling factor k to be any positive integer.For the initialization circuit U I , we append the U I U † I pairs, which means k can only be odd-numbered.To access more choices in k and improve the statistics, for an even-numbered k, we perform both cases of k=2n−1 and k=2n+1 for the initialization circuit and take the averaged outcome.
Our experiment is performed on a superconducting quantum processor with a ring of 16 transmon qubits (Fig. 1d) mounted inside a dilution refrigerator.A tunable coupler -also a transmon qubit -connects each pair of neighboring qubits and serves as an independent control knob for adjusting their coupling strength.A shared control line is used for delivering both the single-qubit gate signal to the qubit and the two-qubit gate signal to an adjacent coupler, reducing the wiring efforts both on the chip and inside the refrigerator.More details about the device can be found in the Methods section and its 8-qubit predecessor [45].In this experiment, we use the coupler-assisted adiabatic controlled-Z (CZ) gate [46] in conjunction with the net-zero pulsing technique to suppress dephasing from low-frequency flux noise [47].The CZ gate fidelities when simultaneously applied are nearly 99% on average.
The VQE algorithms and the ZNE error mitigation often involve long circuits with parallel gate operations.In practice, it is commonly seen that the performance of modular single-and two-qubit gates is inconsistent between calibration and final implementation, as the control signals are distorted due to the spatial and temporal crosstalk effect that varies with circuits.To combat this,

2-q
2-q 1-q 1-q we implement a generally applicable compilation strategy named isomorphic waveform.The key idea is to maximally match the waveforms used during calibration and subsequent implementations by rearranging an arbitrary quantum circuit into interleaved layers of singleand two-qubit gates and translating them into well-timed pulses (Figs.2a-c).First, for two-qubit gates that cannot be parallelized, one can always split them and insert a layer of identity gates between them (identity filling).Then, multiple concatenated single-qubit gates are squeezed into a single SU(2) operation.Note that we also perform identity filling to the "no gate" case.Now with one SU(2) operation for every qubit during each cycle, we apply the U3 decomposition [48] to compile the SU(2) operation into two physical π/2 pulses and additional virtual Z gates, resulting in restless single-qubit gate cycles.Since the virtual Z gates only affect the pulse phases, the entire waveform shape is almost identical given an arbitrary original circuit.
A direct consequence of this compiling scheme is that only the π/2 pulses are required to calibrate.More importantly, simultaneous calibration of these π/2 pulses immediately accounts for microwave signal crosstalk because all pulses have fixed shapes and hence fixed crosstalk influence; this substantially reduces the efforts for performing crosstalk compensation.In addition, the timing of both single-and two-qubit pulses, is fixed during calibration and subsequent experiments.In this way, pulse distortions -if not completely removable -can be more consistent in different circuits; this makes the gate performance more predictable for long circuits.Figure 2d shows the circuits used for calibrating the controlled phase and single-qubit phase correction.With the isomorphous waveform method, their compiled waveforms are almost identical to each other and to the VQE circuits.In addition, the method is compatible with randomized compiling [49], another powerful error mitigation technique.The method may require more singlequbit gates to be added to the circuits.However, this is unlikely to degrade the whole circuit fidelity due to the fact that the SU(2) squeeze operation can reduce the circuit depth and the fact that the single-qubit gate fidelities are typically much better than the two-qubit ones.
After the processor is calibrated, we implement the full energy spectroscopy of a 4-qubit Heisenberg spin chain on the selected qubits (Fig. 3a).Since all the eigenstates have definitive s, s z , and m numbers, the ground state |E 1 ⟩ must be a global singlet (i.e.s=s z =0) which we denote as |S 1 ⟩=|E 1 ⟩.Next in the spectrum are two triplets (i.e.s=1), each with degeneracy three corresponding to s z =0, ±1, which we denote as |T To solve the singlet |S 1 ⟩ which is the ground state, one can simply use the s-conserving circuit (3 parameters) in combination with mirror symmetry (to subtract 1 parameter) illustrated in Fig. 1b.Therefore, there are only 2 parameters per layer.We prepare the initial state Since the initial state has a total spin s=0, the final state remains in the singlet manifold.A simple average energy minimization with single weight (P =1) in the cost function (1) can estimate E S1 .The result converges to a precise value for a circuit with L=2 layers.Figure 3b shows the convergence of the cost function versus optimization iterations and the energy estimation as a function of noise scaling factor k. The ground state energy extrapolated at k=0 yields E S1 =(−6.570 ± 0.059) (theory: −6.464) assuming J=1.Note that our approach requires only 4 parameters to train, while the widely used hardware efficient ansatz [26] requires a total of 3 layers and 32 parameters.As shown in Ref. [36], the situation gets worse as the system size or the number of target states in SSVQE increases.
For singlet |S 2 ⟩, we resort to the SSVQE method by training |S 1 ⟩ and |S 2 ⟩ simultaneously with double initial states (P =2 in the cost function).Unlike the ground state, producing two orthogonal initial states with s=0 requires long-distance controlled-NOT gates which are not available in our device.We instead use the s zconserving circuit with L=3 layers and two initial states To ensure that the final states fall into the s=0 manifold, an additional penalty term ⟨S 2 tot ⟩ is included in the cost function (see Supplementary Information for details).Due to the complex nature of the cost function, the training becomes much harder than the ground state and takes about 24 hours until convergence.We only obtain one successful convergence (Fig. 3b) out of a few trials as the processor performance fluctuates over time.Interestingly, we find that in this case, the hybrid optimization of Nelder-Mead (first 14 iterations) followed by Adam performs better than either of them separately.
For the three triplet groups T 1,2,3 , the eigenstates can be efficiently solved with fewer training tasks by identifying symmetries in these states and choosing proper initial states.To understand this, let us focus on the case of s z =0.We use the total spin s-conserving circuit with mirror symmetry with L=2 layers.For the initial state 2 which is determined by symmetry numbers s=1, s z =0, m= − 1, minimizing the average energy yields the lowest eigenstate in the same-symmetry manifold, which is exactly 2 ⟩.In this manner, all three triplet states s z =0 share a common trained circuit, and error mitigation is conducted separately (Fig. 3c).The same treatment also works for the manifold of s z =±1.In Supplementary Information, we give the exact form of all the initial states and their corresponding circuits.
The last manifold of the spectrum is a quintuple with s=2.Two eigenstates with s z =±2 are easy to obtain as all the spins are aligned.Therefore, we focus on generating the three eigenstates with s z =0, ±1.Generating initial states with total spin s=2 and s z =0, ±1 is not easy without long-distance controlled-NOT gates.Instead, we add the penalty term ⟨(S 2 tot − 6) 2 ⟩ to the cost function with P =1 to guarantee the convergence to the s=2 manifold.We use the s z -conserving circuit with mirror symmetry and the initial states |0101⟩ (s z =0), |0111⟩ (s z =−1), |0001⟩ (s z =1).The results are shown in Fig. 3d.Again, we find that similar to the SSVQE simulation of |S 2 ⟩ (Fig. 3b), the hybrid optimization of the gradient-free Nelder-Mead (first 14 iterations) and gradient-based Adam method gives the best convergence.In fact, when the cost function contains several terms, the potential landscape becomes complex with multiple local minima.It is likely that a few iterations by a gradientfree optimizer, such as Nelder-Mead, bring the cost function to the vicinity of its global minimum whereas a following gradient-based optimizer, such as Adam, facilitates the convergence.This can be seen from the immediate drop in the cost function when switching from Nelder-Mead to Adam at the 15th iteration in all the relevant cases.The estimated energy eigenvalues of all the 16 states are compared with the theory in Fig. 3e.The average deviation is about 0.13, validating our approach for obtaining the full variational spectroscopy.
To show the scalability of our approach, we perform low-energy spectroscopy for an 8-qubit Heisenberg Hamiltonian.We consider estimating the three lowest   trace is dropped for clarity.The three estimated eigenenergies and their theoretical values are compared in Fig. 4b.While E S1 and E T1 can be estimated very precisely, the energy E T2 shows a larger error.This is because a circuit with L=4 layers is not capable to target two eigenstates simultaneously with high precision.
The key advantage of programmable digital quantum simulators is their capability in simulating different Hamiltonian models.To show the versatility, we also perform on the same device the variational spectroscopy for the transverse-field Ising Hamiltonian where J is the coupling strength and h is the transverse field strength.This system has a critical point at h=J where the ground state becomes highly entangled, and thus its variational simulation is expected to be harder [40,50].The Hamiltonian commutes with two operators, namely the mirror operator M and the parity operator Π z = N i=1 σ i z , which imply that every eigenstate has two symmetry numbers, m=±1 (i.e.M |E j ⟩ =m |E j ⟩) and π z =±1 (i.e.Π z |E j ⟩ =π z |E j ⟩).Hence, We can denote the state as |E (m,πz) j ⟩.For simulating the transverse Ising Hamiltonian, the two-qubit operation N (θ) is replaced with The circuit used to construct N I (θ) is depicted in Fig. 1b.To ensure the preservation of the Π z symmetry, the single-qubit operations should all be z-rotations.Furthermore, mirror symmetry is readily preserved by choosing parameters symmetrically with respect to the middle of the chain.The three lowest energy eigenstates are denoted as |E (1,1) 1

⟩ and |E
(1,−1) 3 ⟩.These three eigenstates possess symmetry numbers that allow their independent generation through the appropriate choice of the initial state and a cost function with only P =1.We simulate these three eigenstates for both 4-qubit and 8-qubit Ising chains, presenting the result for the 8-qubit case here.The 4-qubit results and other experimental details can be found in Supplementary Information.
In our study, we present an experimental demonstration of multi-level variational spectroscopy using a programmable superconducting quantum simulator.Our findings highlight the resource efficiency afforded by symmetry and the ability to target a greater number of states through the implementation of the subspace search method and the inclusion of penalty terms.However, a significant challenge is the circuit depth overhead, which is inherent to the SSVQE algorithm when the number of target states increases.Furthermore, the limited qubit connectivity of our one-dimensional simulator imposes constraints on generating appropriate initial states.Both increasing qubit connectivity by utilizing simulators based on two-dimensional qubit arrays or implementing the protocol on alternative platforms such as ion traps should address these limitations.Additionally, we discover that a hybrid optimization strategy, which combines gradient-free Nelder-Mead with gradient-based Adam, offers faster convergence when penalty terms are present in the cost function.This suggests that in a disrupted potential landscape caused by relatively strong penalty terms, Nelder-Mead is less sensitive to local minima, guiding the cost function toward its global minimum.At this point, the Adam optimizer can rapidly converge.

Methods
VQE training.To ensure that our convergence is accurate, we use a restricted range of the initial parameters, which are uniformly sampled from -π/2 to π/2.These parameters are used as a starting point for the VQE training procedure.The variational circuit with initial parameters is then compiled into isomorphic forms as described in the main text and further discussed in Supplementary Information.At each iteration, the cost function is measured by measuring relevant correlation functions.To ensure accuracy, each measurement is repeated around 6000 times.Then by feeding the cost function into our classical optimizer, we obtain a new set of parameters ⃗ θ.We consider three different ways for classical optimization, namely gradient-free Nelder-Mead, gradient-based Adam, and a hybrid of both.In the Nelder-Mead algorithm, the initial simplex is uniformly sampled from -π/2 to π/2.In the Adam optimizer, the initial parameters are also uniformly sampled from -π/2 to π/2.The configuration of Adam optimizer is chosen to be α=0.1 for the learning rate, β 1 =0.9 for the exponential decay rate of the 1st-moment, β 2 =0.999 for the exponential decay rate of the 2nd-moment, and finally ϵ=10 −8 for numerical stability.In the hybrid optimization, we utilize the Nelder-Mead for ∼15 iterations to obtain a rough training result, and then use Adam until the training converges.
In most cases, the Nelder-Mead algorithm performed better than the other two, See Supplementary Information for specifying the best optimizer for every single case.However, when the cost function contains several terms, e.g. in the SSVQE with P >1 or in the presence of penalizing terms, the best optimization method is often the hybrid method.This might be due to the presence of multiple local minima in such a complex potential landscape.In such cases, Nelder-Mead brings the cost function to the vicinity of its global minimum but it is slow to make the final convergence towards the exact location of the global minimum.However, the gradient-based Adam optimizer operates well when the cost function is close to its global minimum.As a result, for complex multi-term cost functions, hybrid optimizations seem to provide faster convergence.
Error mitigation.In our experiment, we utilize the initialization circuits to prepare initial states that satisfy the corresponding S z value and ansatz circuits to conserve S z .To ensure that the conservation of S z is maintained, we remove the data that fails to meet this criterion.
To mitigate errors, we employ the zero-noise extrapolation (ZNE) technique, which involves proportionally increasing the entire circuit error.For the initialization circuit U I , we increase the circuit error by a factor of k (where k=2n+1 is an odd number) by replacing the cir-cuit U I with (U I U † I ) n U I .However, if k=2n, we take the average of (U I U † I ) n−1 U I and (U I U † I ) n U I to obtain the desired results.For the ansatz circuits, we can increase the circuit error by a factor of k by substituting the original N (θ) with N (θ/k) k .Such a protocol allows us to utilize even noise scaling factors, improving the statistics during extrapolation given a limited circuit length.
Device and setup.The sample (10mm-by-10mm in size) is made of aluminum on a sapphire substrate and packaged in an aluminum sample enclosure.The sample is mounted inside a Bluefors LD400 dilution refrigerator at a base temperature of 12 mK and shielded with two layers of cryo-perm cans.
The qubits have alternating frequencies between 3.8 GHz and 4.2 GHz, designed for suppressing residual ZZ interactions.The couplers idle at their maximum frequencies of 6.4 GHz and their couplings to the qubits are about 100 MHz.The readout resonators designed at ∼7 GHz have photon loss rates of 1 MHz and a dispersive shift of 0.6 MHz.Eight resonators share a common transmission line to which a 20 fF capacitor is added at the input end to prevent photons from entering this port.
The XY microwave signal for single-qubit gates is synthesized by up-converting intermediate-frequency signals with a single-tone microwave carrier using an IQ mixer.To ensure high signal quality without reflections and spurious frequencies, an isolator and a band-pass filter are used.The fast Z signal for two-qubit gates is combined with the XY signal using a diplexer at room temperature.The combined signal is then attenuated and filtered at multiple cold stages of the DR for noise suppression.A custom-made low-pass IR filter is added to each control line at the mixing-chamber stage which attenuates signals at the qubit frequency by approximately 25 dB while allowing low-frequency signals to pass through.
The multiplexed readout signal is generated from upconversion using an IQ mixer, filtered by a 7-GHz bandpass filter, and attenuated at different temperature stages before entering the device.The output signal is sequentially amplified by a HEMT amplifier at 4K and an amplifier at room temperature before being down-converted and demodulated.To block noise from higher temperature stages, circulators and filters are used at the output port.More device and setup information can be found in Supplementary Information.

INITIAL STATE CONSTRUCTION
In order to approximate the desired eigenstate with the symmetry-preserving ansatzes in VQE, the initialization needs to be carefully chosen for the corresponding symmetry.For the Heisenberg Hamiltonian with system size N =4, the two singlet eigenvectors |S 1,2 ⟩ share the symmetry of total spin s=0, z-magnetization s z =0, and mirror symmetry m=1.In order to target |S 1,2 ⟩, two initial states with these symmetry numbers are required.These states can be obtained by Clebsch-Gordan decomposition as |ψ − ⟩⊗|ψ − ⟩ and ).The initialization circuit for the former can be easily obtained which is depicted in Table S1.However, the latter requires a non-trivial initialization circuit with long-distance controlled-not gates which is not available in our device.Therefore, we use initial states |0101⟩ and |1010⟩, see their corresponding quantum circuits in Table S1, with S zconserving circuit to approximate S 2 using SSVQE.
The three triplets |T ⟩.Therefore, for each given s z , we use the following initial states: For s z =0 : For s z =+1 : For s z =−1 : The quantum circuits for generating these states are shown in Table S1 and Table S2.
For the quintuple |Q (0,±1,±2) ⟩, which is determined by total spin s=2, there are five degenerate eigenstates given by s z =−2, −1, 0, +1, +2.To generate these eigenstates, we use a s z -conserving circuit together with a penalty term added to the cost function, see the main text for more details.Therefore, the initial states for these five eigenstates are simply set to be |1111⟩ , |0111⟩ , |0101⟩ , |0001⟩ , |0000⟩.The corresponding initialization circuits are shown in Table S2.
In the case of Heisenberg Hamiltonian with system size N =8, we simulate the 3 lowest eigenenergies.Similar to the N =4 case, the ground state |S 1 ⟩ is a global singlet with total spin s=0, z-magnetization s z =0, and mirror symmetry m=1.The initial state for the ground state is |ψ − ⟩ ⊗4 .For the first excited state 2 ⟩ share the same symmetry numbers as the first excited state, its initial state can similarly choose as |ψ − ⟩ |ψ − ⟩ |ψ + ⟩ |ψ − ⟩.The initialization circuits for these three initial states are shown in Table S3.For other eigenstates with s z =±1, one can simply replace |ψ + ⟩ with |00⟩ or |11⟩.The circuits are not shown, as they can be easily obtained by replacing the circuit for |ψ + ⟩ with the circuit for |00⟩ or |11⟩.
The symmetry numbers for the low-energy spectrum of Ising Hamiltonian with 4-qubit and 8-qubit systems are the same.Namely, the ground state |E 1 ⟩ is determined by z=1 for the parity operator Π z , and m=+1 for mirror symmetry M. Hence, its initial state for both system sizes N =4 and N =8 is set to be |0⟩ 2 .The second eigenstate, |E 2 ⟩, exhibits z=−1 and m=−1 symmetries, with the initial state given by |0⟩ The third eigenstate, |E 3 ⟩, has z=−1 and m=+1 symmetries, and its initial state is expressed as |0⟩ The initialization circuits for these three eigenstates are shown in Table S4 and Table S5.

DEVICE AND MEASUREMENT SETUP
The processor is made of aluminum on sapphire using a recipe similar to that in [1] and bonded in an aluminum sample box which is then mounted inside a Blufors LD400 dilution refrigerator with a base temperature of about 10 1 S3.Circuits for the 8-qubit Heisenberg model

State Initialization Circuit
Ansatz layer # of layers mK.The sample box is shielded by two layers of cryo-perm cylinders.Figure S1 depicts the schematic diagram of the measurement setup.
The XY microwave signals used for single-qubit operations are generated using the IQ mixing.A bandpass filter at the qubit frequency and an isolator are used for suppressing the spurious and reflection.Fast Z pulses used for two-qubit gates are combined with the XY pulses at room temperature using a diplexer, and then attenuated and filtered at multiple temperature stages for further noise reduction.A customized low-pass IR filter is added to the control lines at the mixing-chamber stage which attenuates signals at the qubit frequency by approximately 25 dB while allowing low-frequency signals to pass through, serving as an equalizer.
The multiplexed readout signals are attenuated and filtered before being fed into the readout transmission lines.The output signals are then sequentially amplified by a HEMT amplifier at 4 K and a room temperature amplifier  The average flux fluctuation in our device is estimated to be approximately 60 µΦ 0 .The measured ZZ interaction between neighboring qubits is around 10 kHz and therefore considered negligible.More device details are listed in Table S6.

SINGLE-QUBIT GATES
The microwave signal or XY crosstalk is a major source of errors when performing parallel single-qubit gates.In our experiment, we use a hybrid strategy to combat the XY crosstalk.For qubits with significant crosstalk and with frequencies close to each other (usually the next-nearest-neighbor qubits), we calibrate and add the compensation signal for crosstalk cancellation [2].Aside from those, we rely on the isomorphous waveform strategy presented in the main text which compiles the circuit into restless single-qubit cycles and automatically corrects for crosstalk.The flux or Z crosstalk can be effectively cancelled by adding compensating signals.The statistics of both XY and Z crosstalk -4% for XY and 0.2% for Z (median value) -in our device are shown in Fig. S3.
For single-qubit gates, we employ the U3 decomposition [3] with the need for calibrating X π/2 pulses only.We use 30-ns-long X π/2 .We add an additional 10-dB attenuation at room temperature to rescale the pulse amplitudes to approximately half of the maximum DAC output range for better signal quality.To simultaneously calibrate these X π/2 gates, we use tbe pulse train technique which progressively tunes up the pulse amplitude and DRAG coefficients with increasing cycle number.We characterize our calibrated single-qubit gate errors using simultaneous cross-entropy benchmarking (XEB), with the result shown in Table S6.The simultaneous single-qubit gate fidelities are limited by T 2 times and residual crosstalk effect.

TWO-QUBIT GATES
The CNOT gates used in our experiment are synthesized from native CZ gates and single-qubit gates.We use the coupler-assisted adiabatic CZ gate scheme demonstrated in our previous work [4].The bipolar net-zero pulse  technique is used here for improving robustness against low-frequency flux noise and minimizing signal distortion.For a unipolar half, we employ a hyperbolic tangent pulse shape described by where A and τ are the pulse amplitude and width respectively.ϵ is a parameter controlling pulse shape.Pulse profiles with different ϵ parameters are shown in Fig. S4.Larger ϵ value tends to make the pulse rise more rapidly and hence less adiabatic.In our setup, we set all the CZ pulse duration to 60+60 ns (τ = 60 ns for each pole) as a trade-off between leakage and decoherence errors.our pulses.Then, we calibrate the flux (Z) crosstalk among different control lines and cancel it with compensation pulse.The shape parameter ϵ is decided by performig the echoed CZ pulse train as a function of ϵ and pulse amplitude A (Fig. S5(a-b)).The control qubit population reveals the information about the controlled phase while the target qubit population reveals the information about the leakage.Monitoring both helps us determine the approriate value for ϵ.Here we choose ϵ ≈ 2 in our experiment.
Once the shape of the CZ pulse is determined, the only left pulse parameter to calibrate is the pulse amplitude A, which can be progressively tuned using the same echoed pulse train sequence (Fig. S6(a)).A similar procedure is followed to obtain the local Z phase correction for both qubits (Fig. S6(b)).Note that we devide the CZ gates in the 8-qubit chain into two groups (Table S6).These calibrations are simultaneously performed for qubit pairs in a same group.The calibrated CZ gates are simultaneously benchmarked using parallel two-qubit XEB.The average gate fidelities are approximately 0.99 after subtracting the single-qubit gate errors (Table S6).

READOUT
To improve the fidelity of the readout, we employ the shelving technique as described in Ref. [5].This method involves the addition of extra π ef pulses to all qubits before executing the measurement pulse.When the demodulated

POTENTIAL LANDSCAPE
To validate the accuracy of our convergence, we use the 4-qubit Heisenberg state S 1 as an example.The measured potential profiles are compared to theory for 6 different pairs of θ (4 parameters) in Fig. S7a-f.On these figures, we also plot trajectories of two convergence traces, showing how the they reach the minimum in all these dimensions.S4.

FIG. 2 .
FIG. 2. Circuit compilation using the isomorphous waveform method.a-c, Step-by-step compilation.An arbitrary quantum circuit composed of modular single-and two-qubit gates (a) can be rearranged into periodic single-qubit (1-q) and two-qubit (2-q) gate cycles (b).Red arrows indicate the identity filling and single-qubit gate squeeze operations.All SU(2) operations including identity operations are subsequently compiled into two π/2 pulses according to the U3 decomposition (c).d, Gate sequences used in the controlled-phase calibration, the local Z phase calibration, and the VQE experiment for the Heisenberg model in which the initialization circuit (red solid box) and the N (θ) circuit (brown solid box) are indicated.The black dashed box indicates a section of the three circuits that share a common waveform as illustrated in (c).

+1) 2 EnergyFIG. 3 .
FIG. 3. Full spectroscopy of a 4-qubit Heisenberg spin chain.a, The left plot shows the device connectivity and the four active qubits (solid circles).The right plot sketches the spectrum of the 16 states of the 4-qubit Heisenberg spin chain which are grouped by degeneracy.b-d, VQE training and error mitigation for each of the 16 states including the singlets S1 and S2 (b), the triplets T1, T2 and T3 (c), and the quintuple Q (d).The top panels (gray background) in each figure plot the training curves of the cost function C( ⃗ θ).The blue dots are the mean values of the training curves with different initial parameters.The shades indicate twice the standard deviation.The bottom panels (colored background) plot the error mitigation results for different training outcomes (green dots) and the exponential fit to their mean (green line).The black arrows indicate the theoretical values.Note that the triplet states with the same sz values share a common trained circuit and their energies can be obtained simply by switching the initialization circuit.e, Comparison of the extrapolated zero-noise energies (red bars) and the theoretical values (dashed boxes) for all 16 eigenstates.The error bars (blue bars) denote plus/minus twice the standard deviation.

FIG. 4 .
FIG. 4. Low-energy spectroscopy of an 8-qubit Heisenberg and transverse Ising chain.a, VQE training (left panels) and error mitigation (right panels) for the three lowest eigenstates for the 8-qubit chain with the Heisenberg interactions.b, Comparison of the extrapolated zero-noise energies (red bars) and the theoretical values (dashed box) of the Heisenberg model.The error bars (blue bars) are plus/minus twice the standard deviation from the measured values.c, Same as (a) but with the transverse Ising interactions.d, Same as (b) but for the Ising model.The center panel at bottom shows the eight active qubits (solid circles).

3 ⟩ 3 ⟩
share the symmetry of total spin s=1.Each of these eigenstates is triply degenerate with s z =0, ±1.The corresponding mirror symmetry are m=−1 for |T (0,±1) 1,, and m=+1 for |T (0,±1) 2 conversion and acquisition by the digitizer.Repeated characterization result of the T 1 relaxation time is shown in Fig. S2.The T 1 times of the eight selected qubits used in our experiment generally ranges from 40 µs to 80 µs, each showing a typical fluctuation of 10-20 µs.
FIG. S1.Schematic diagram of the measurement setup.One set of the readout resonator, qubit, tunable coupler, the shared control line, and corresponding wiring is shown.

FIG
FIG. S2.Qubit T1 time.Cumulative distribution of the repeated measured T1 relaxation times for the selected eight qubits (inset).
FIG. S3.Crosstalk characterization.Cumulative distribution of measured microwave (XY) and flux (Z) crosstalk ratios with the median values indicated.

THE 4 -
FigureS8shows the experimental results for the 4-qubit transverse Ising model.With similar training and error mitigation procedures, we obtain the three lowest energy states E 1 , E 2 and E 3 .The circuits used in the experiment are listed in TableS4.
FIG. S5.CZ pulse shape and leakage.a, The measured populations of the target qubit q0 versus the shape parameter ϵ and the pulse amplitude δA referenced to a certain value under the echoed CZ pulse train sequence shown above.Here the number of repeated cycles is 4. The red fringe in the middle corresponds to a controlled phase of π. b, The measured populations of the control qubit q1 under the same sequence described in a.The region with more population (red) indicates more leakage.
FIG. S7.Energy landscape of the S1 state in the 4-qubit Heisenberg model.a-f, The measured and theoretical values of the cost function in the S1 case as a function of different pairs of variational parameters selected from the set {θ0, θ1, θ2, θ3}.On the figures are plotted two convengence trajectories (orange and magenta).

E 1 E 2 E 3 E 1 E 2 E 3 ab
FIG. S8.Low-energy spectroscopy of a 4-qubit Ising spin chain.a, VQE training (left panels) and error mitigation (right panels) for the three lowest eigenstates for the 4-qubit chain with the transverse Ising interactions.b, Comparison of the extrapolated zero-noise energies (red bars) and the theoretical values (dashed box) of the transverse Ising model.The error bars (blue bars) are plus/minus twice the standard deviation from the measured values.

TABLE S1 .
Circuits for the 4-qubit Heisenberg model

TABLE S4 .
Circuits for the 4-qubit transverse Ising model

TABLE S5 .
Circuits for the 8-qubit transverse Ising model Qubits used in the 4-qubit experiment.See Fig.S2for the indexing information of the qubits.Both the |2⟩ and |1⟩ state have been treated as the |1⟩ state with the shelving protocol.The ZZ interaction strength between neighboring qubits is measured when couplers are biased at their idling point.

TABLE S6 .
Device parameters and gate fidelities.

TABLE S8 .
State properties and VQE information for the 4-qubit and 8-qubit transverse Ising model EnergyNoise scaling factor k