Simulating Models of Challenging Correlated Molecules and Materials on the Sycamore Quantum Processor

Simulating complex molecules and materials is an anticipated application of quantum devices. With the emergence of hardware designed to target strong quantum advantage in artiﬁcial tasks, we examine how the same hardware behaves in modeling physical problems of correlated electronic structure. We simulate static and dynamical electronic structure on a superconducting quantum processor derived from Google’s Sycamore architecture for two representative correlated electron problems: the nitrogenase iron-sulfur molecular clusters and α -ruthenium trichloride, a proximate spin-liquid material. To do so, we simplify the electronic structure into low-energy spin models that ﬁt on the device. With extensive error mitigation and assistance from classical recompilation and simulated data, we achieve quantitatively meaningful results deploying about one ﬁfth of the gate resources used in artiﬁcial quantum advantage experiments on a similar architecture. This increases to over half of the gate resources when choosing a model that suits the hardware. Our work serves to convert artiﬁcial measures of quantum advantage into a physically relevant setting.


I. INTRODUCTION
There has been much interest in simulating problems of chemistry and materials science on quantum computers [1,2].This is not least because the first claims of potential quantum "advantage" in a number of tasks have appeared [3][4][5].Such simulations have deployed an impressive level of quantum resources; the random-circuit-advantage experiments on the Sycamore processor used more than 50 qubits and approximately 500 two-qubit gates.However, the common feature of recent attempts to reach quantum advantage is their artificial nature.In particular, the experiments propose a special selection of tasks and metrics designed to facilitate large quantum advantage on current hardware.For example, Ref. [3] involved distinguishing a theoretical signal of 1 from 0, but a noisy experimental result of 0.002 was already been sufficient to claim success at the time due to the comparative complexity of the chosen classical algorithm to obtain the same result.The extent to which success in such a metric is transferable to generally accepted criteria for successfully simulating actual chemical and materials problems is unclear.
In this work, we study the simulation of two "realistic" (if still highly simplified) problems of strongly correlated molecular and materials electronic structure on a stateof-the-art superconducting quantum processor.Aside from the necessary simplifications, neither the systems nor the observables studied are chosen to favor the demonstration of advantage on the experimental hardware.The current work thus attempts to report on the ability of current quantum devices to tackle problems of real-world interest without careful preselection.In particular, we examine the relation between the deployable quantum resources to achieve artificial quantum advantage metrics and those that can be used in successful electronic structure simulations.
The first problem we consider is the low-energy electronic structure of the iron-sulfur clusters of nitrogenase, including the Fe-Mo cofactor, critical components of the natural nitrogen cycle.The second is the electronic structure of α-RuCl 3 , a candidate material for realizing spin-liquid physics.Both are strongly correlated electron structure problems that contain challenging features for quantitative heuristic quantum chemistry approximations.Although for practical reasons we need to reduce the electronic structure into low-energy spin models, the range, nature, and topology of the interactions provide some elements of real-world complexity that can be absent from more specially chosen simulation problems.We target finite-temperature static and dynamical properties representative of accessible physical chemistry measurements and study both systems using the Weber superconducting qubit processor derived from Google's Sycamore architecture, using the finite-temperature version of quantum imaginary time evolution with global classical recompilation of the circuits [6,7].To obtain meaningful data, many types of error mitigation are necessary; thus we discuss the impact of different protocols.Our results illustrate to what extent current superconducting quantum processors can be used to simulate real-world chemical and materials problems.

Fe-S clusters and nitrogenase
Nitrogenases are enzymes that convert atmospheric dinitrogen into ammonia.The process involves the coordinated transfer of multiple electrons and protons to dinitrogen and utilizes multiple metalloclusters found in the nitrogenase enzyme: the [4Fe-4S] Fe-cluster, the [8Fe-7S] P cluster, and the [7Fe-1Mo-9S] Fe-Mo cofactor (the latter can also be found with other metals replacing molybdenum) [8].The electronic structure of the Fe-S clusters remains incompletely understood; theoretical calculations on the Fe and P clusters have unveiled a large number of low-lying spin states [9,10], the role of which in the electron-transfer process is unclear.In the case of the Fe-Mo cofactor, the oxidation and spin states of the ions are unresolved.
The simplest electronic models of these compounds involve the Fe and S valence orbitals, with over 100 spin orbitals in the case of the Fe-Mo cofactor, too large for current quantum devices [11].However, both theoretical and spectroscopic studies indicate that the electron-hole-like excitations mostly lie at higher energies than a manifold of low-energy spin-coupled states.Thus the low-lying electronic spectrum can be captured using a variety of spin models, which are commonly used to interpret the electronic structure [9,12,13] as well as to model measured quantities, such as the magnetic susceptibility and heat capacity [14,15].In the most commonly studied Fe-S models, the metals occupy vertices and the S bridges form the edges, with the latter enforcing antiferromagnetic interactions parametrized by the exchange coupling J [9,12].
The three-dimensional models and the two-dimensional pattern of the spin couplings are shown in Fig. 1.Although the metal ions have large total spin (e.g., S = 2, S = 5/2), we can further simplify (reducing the quantum resources) by representing each metal by a S = 1/2 spin, retaining only the topology of the spin-spin couplings, and with an isotropic Heisenberg exchange interaction between each spin pair.Although the spectrum is changed by this reduction, it retains similar features in the lowest few states, e.g., the total spin value and degeneracy, as discussed in Appendix A. The simplified Hamiltonian can thus be written as where S i and S j are S = 1/2 operators and ij denotes each of the bonds in Fig. 1.We note that other instances of Heisenberg models have been simulated on quantum devices (see, e.g., Refs.[16,17]) but the interesting properties of the FeS clusters relate to the specific topology and magnitudes of the spin couplings.
The magnetic model parameters in the FeS clusters are usually derived through fitting to experimental data or through ab initio calculations [9,12].In the Heisenberg model, one parameter can be set as the scale of energy (J ), while the antiferromagnetic nature of the interaction induced by the Fe-S-Fe bridges fixes the sign of the couplings to be positive.In the case of [4Fe-4S], we use only two different J ij , namely J = 1 and J = 1.17J , as previously determined to characterize the all-ferric cubane [18] [J couples (1-2) and (3-4), while J couples the other pairs, as shown in Fig. 1(b)].Deriving accurate exchange couplings in the P cluster and the Fe-Mo cofactor is more challenging; thus, noting the modest variation in J ij in [4Fe-4S], we use a single J for all couplings in the P cluster and the Fe-Mo cofactor (all J positive).With this simplification, the P cluster and the Fe-Mo cofactor have the same spin Hamiltonian.(Note that the units chosen for energy define an inverse energy unit for time; all times will hence be assumed to be in such units.)

Ruthenium trichloride, α-RuCl 3
Ruthenium trichloride is a transition-metal material of interest as a "proximate" spin liquid [19].Similarly to in the Fe-S clusters, the low-energy excitations are spin excitations.In particular, the edge-sharing octahedral coordination around the Ru(III) ions (see Fig. 1), together with the spin-orbit coupling, leads to a low-energy Hamiltonian that approximates the exactly solvable Kitaev model in the spin-liquid regime.However, the degree of similarity to the Kitaev model, as well as the interpretation of spectroscopic and thermal measurements in terms of modifications to the Kitaev model, is much debated.Much work has been devoted to deriving and parametrizing low-energy Hamiltonians for α-RuCl 3 [20].The simplest family of Hamiltonians are the Kitaev-Heisenberg models, which take the form where γ = X , Y, Z.We also use a parametrization in the literature from Refs.[20][21][22], where the above form is augmented by additional couplings: where γ (α, β) = X , Y, Z, and α and β take indices different from γ .Specifically, we use the parameters: J = −1.53,K = K γ = −24.4,= 5.25, and = −0.95.The Kitaev point, which is exactly solvable, corresponds to taking only the K γ term in the above models.It is argued that in the parameter regime of α-RuCl 3 , both the excitations and heat capacity show echoes of the two kinds of Majorana fermions that exist at the exactly solvable point [23,24].

Formal algorithm
We simulate the finite-temperature energies and dynamical correlation functions of the systems.The quantum imaginary time evolution (QITE) algorithm is used to prepare a sample of the finite-temperature state ρ = e −βH /Tr(e −βH ).The basic idea in QITE is to prepare normalized imaginary time-evolved states is the unitary determined by the QITE procedure, and | i (0) is the initial state (here, the ith computational-basis state).We are interested in the finite-temperature static and dynamical observables A β = Tr[ρA] and A(t)B = Tr[ρA(t)B], computed as where ) .The expectation values A i are obtained by reading out Pauli operators (static observables) or from the Hadamard test (dynamic observables).
In the original QITE procedure, the unitaries are determined at each time step from a set of linear equations constructed from measurements on the quantum device.Similarly, the normalization weights P i are accumulated at each time step.Once the imaginary time states | i (β/2) are prepared, they can then be propagated in real time to yield states | i (β/2 + t) , with the real-time propagation unitary U(t) generated, e.g., by a Trotter evolution of the Hamiltonian.The schematic circuit for a dynamical simulation is shown in Fig. 1(a).In the current experiments, however, we use classical recompilation [7] to generate both the circuits U(β/2) and U(t) and the weights P i .The details are described further below.

Implementation and error mitigation
Hardware and qubit selection.The simulations are run on Google's 53-qubit Weber processor based on the Sycamore architecture [3].The Fe-S cluster simulations use up to five and nine qubits, respectively, while the α-RuCl 3 simulations use up to seven and 11 qubits, respectively.(In the above qubit counts, we include the one ancilla qubit used for the Hadamard test circuit.)The best-performing qubits are selected using a combination of single-qubit randomized benchmarking, two-qubit cross-entropy benchmarking, and Loschmidt-echo metrics on the hardware; an example of the embedding of the qubits onto the Weber architecture for the Fe-S clusters is shown in Fig. 1(e).
Circuit recompilation.The circuit realization of the formal QITE algorithm [including U(β/2) and U(t)] using a reasonable Trotter time step (e.g., 0.1 inverse energy units) is far too deep for current quantum devices.For example, using a first-order Trotter expansion with a time step of 0.1 in the six-site Kitaev-Heisenberg model at the Kitaev point, we already exceed the two-qubit gate count used in random-circuit advantage after three units of time.Thus, to reduce the circuit depth, for each desired imaginary-time or imaginary-plus real-time point at which observables are to be measured, we define the corresponding U i by classical recompilation.This is a variational procedure whereby the exact unitary is constructed classically and then a circuit of given depth ( Ūi ) is (classically) variationally optimized to maximize the fidelity [7].We note that classical optimization of the global fidelity is not a scalable procedure but it is necessary to produce the best results on quantum hardware; optimization to reproduce marginals (density matrices) could, in principle, be used if scalability was important [7]  Postselection.The Fe-S clusters and simplified α-RuCl 3 Hamiltonians possess Z 2 symmetry.We use this to perform postselection as discussed in Ref. [7].
Floquet calibration.As discussed in Refs.[25,26], one can calibrate an "excitation-number"-conserving gate in terms of five angles: θ, φ, ζ , γ , and χ .The ideal √ iSWAP † gate should have θ = π/4, with all other angles being zero, but in practice this is not the case.We perform "Floquet" calibration [25,26] before each experiment to calibrate the actual angles for each run.The modified √ iSWAP † gate is then used in the classical recompilation procedure to obtain the most faithful compressed circuit.
It is important to note that this calibration is only performed on isolated √ iSWAP † gates.However, due to the specifics of the hardware, the presence of PhXZ gates ("microwave" gates) can lead to additional errors in the two-qubit gates that are not accounted for in the calibration.
Dynamical decoupling.In the schematic shown in Fig. 2(a) for dynamical observables, the ancilla qubit in the Hadamard test is idle for large parts and can decohere.
To mitigate this effect, we employ a dynamical decoupling sequence [27,28] consisting of inserting identities generated by XX and YY gates on the ancilla.
Rescaling.To improve the data, we perform postprocessing, rescaling all dynamical correlation functions.To do so, we perform a classical noiseless simulation using a modified Hamiltonian Ĥ containing only commuting terms on pairs of qubits.For example, for the hexagonal model shown in Fig. 1(d), we use only the Hamiltonian terms on (1-2), (3)(4), and (5-6).The commuting form of Ĥ means that the exact classical simulation can be performed easily.We then perform the same simulation on the quantum device, generating the circuit by recompilation with the same ansatz as for the full Hamiltonian and including all the error-mitigation techniques described above.It is important to note that even though Ĥ contains only commuting terms, the ansatz used for the quantum device does not have this simplifying commuting structure-e.g., in Fig. 1(b), the √ iSWAP † gates couple all the qubits-thus the quantum simulation of Ĥ samples similar noise to that of Ĥ .For each time point t (in imaginary and in real time), we then define a rescaling factor where A ideal t ( Ĥ ) is obtained using classical noiseless emulation and A hw t ( Ĥ ) is obtained from the hardware.f (t) is then used to rescale the dynamical correlation functions of the real Hamiltonian.

Fe-S clusters
Figure 3 shows simulation results on the Fe-S clusters.Despite the simplifications of the model (e.g., to S = 1/2), the static Z i Z j correlation functions at the lowest temperature (β = 2) for the [4Fe-4S] and P and/or Fe-Mo-cofactor clusters reflect the known ground-state pairing patterns of the spins in the true [4Fe-4S] and P clusters [Figs.3(a) and 3(c)].The largest error in the correlation functions is approximately 22% in the P and/or Fe-Mo-cofactor cluster.The majority of this error is from the hardware rather than the classical recompilation; the error from classical recompilation can be seen in the difference between the red and black lines in Fig. 3(c).The total recompiled circuit complexity for the P and/or Fe-Mo-cofactor cluster dynamical correlation functions is 82 two-qubit gates.The larger circuit complexity is reflected in the somewhat lower quality of the dynamical correlation functions.Up to frequency ω = 3 (units of energy), there are seven identifiable peaks in the spectrum.Encouragingly, the hardware results even before rescaling capture the correct frequency of most of these peaks; however, the amplitudes are severely degraded: the largestamplitude peak is only about 40% of the expected height.Postprocessing of the data (rescaling as well as shifting to satisfy the imaginary part to zero mean) significantly improves the results, although the largest-amplitude peak still has an amplitude error of 20%.ancilla qubit, the latter system corresponds to the largest simulations that we perform.We start with six-site simulations of the heat capacity, obtained by numerically differentiating the finitetemperature energies.The energy and heat capacity are shown in Fig. 4(a).One of the main features of the proximate spin-liquid behavior of α-RuCl 3 is a two-peak structure in the heat capacity.While the energy has a significant error at lower temperatures, the two-peak structure of the heat capacity is visible at around T = 1 and T = 6 − 7, although it is extremely noisy, as the numerical derivative amplifies the noise.

Ruthenium trichloride, α-RuCl 3
Figures 4(b) and 5(a) show the dynamical correlation functions for the α-RuCl 3 six-and ten-site models (64 and 310 two-qubit gates, respectively), as well as a comparison with a second anisotropic parameter set (with smaller YY and ZZ terms) in Fig. 5 any reasonable physical spectra.Figure 5(b) shows the sensitivity of the quality of simulation to the choice of model, as the more anisotropic Kitaev-Heisenberg simulations are of significantly higher quality.

III. DISCUSSION AND CONCLUSIONS
In the current work, we discuss the quantum simulation of two representative real-world problems: the Fe-S clusters of nitrogenase, including the P cluster and the Fe-Mo cofactor, a problem of interest in correlated quantum chemistry, and the proximate spin liquid α-RuCl 3 , a correlated material.These systems are not selected for their suitability for implementation on the Sycamore platform and for practical computation, simplification of the Hamiltonians into effective low-energy spin models is required.Some reasonable results are obtained for such simplified models of these problems.In the Fe-S clusters and the smaller α-RuCl 3 instance, qualitatively correct features in the spin structure, excited-state spectrum, and heat capacity can be obtained.However, to achieve this, the  implemented circuits need to be obtained with the help of classical recompilation and the data require significant processing, including using data from exact classical simulations of related tractable problems.Unfortunately, these steps raise questions with regard to effectively simulating more classically difficult systems.The main limitation in the experiments is the two-qubit gate count, rather than the number of qubits.Simulations with up to 100 two-qubit gates, such as the larger Fe-S cluster simulations and the six-site α-RuCl 3 simulations, can be carried out with some confidence on the hardware.However, our largest simulations for α-RuCl 3 , which use 11 qubits, 310 two-qubit gates, and 782 single-qubit gates, are not successful.However, we can obtain meaningful simulation data with these quantum resources for Hamiltonian parameters that are tuned away from the α-RuCl 3 regime, indicating the impact of tuning the problem for the characteristics of the hardware.The successfully deployed circuit resources are less than those used in some recent simulation experiments using a similar chip [25].However, the discrepancy may also be understood in terms of how well the experiment matches the capabilities of the chip.For example, the simulations of the free-fermion dynamics in Ref. [25] do not require microwave gates, which helps with Floquet calibration and allows more gates to be used.
The questions of interest in realistic molecular and materials simulations are ones that require some degree of  quantitative precision.Compared to the resources deployable to achieve random-circuit advantage on the same chip architecture, our representative simulation problems can use about one fifth of the qubit and gate resources.If we adjust our models to be more tuned to the hardware, it is possible to use more than half of the gate resources of quantum advantage experiments, while retaining some level of physical accuracy.This provides an understanding of the relevance of artificial quantum advantage experiments to problems of physical simulation and reflects the current status of quantum hardware and quantum simulation.
The data sheet of the Weber device can be found in Ref. [32].The data and scripts for the density-matrix renormalization group (DMRG) simulation in Appendix A can be found in Ref. [33].Other data that support the findings of this study are available from the corresponding author upon reasonable request.The code used to generate the numerical results presented in this paper can be made available upon reasonable request.

FIG. 1 .
FIG. 1.(a) The molecular structures of the [4Fe-4S] and [8Fe-7S] P clusters.(b) The topology of the four-site and eight-site Heisenberg spin models of the (top) [4Fe-4S] and (bottom) [8Fe-7S] P clusters.(c) The crystal structure of α-RuCl 3 : trigonal space group P3 1 12.(d) The topology of the six-site and ten-site Kitaev-Heisenberg spin models of α-RuCl 3 .The adjacent ancilla qubits are indicated as "A."(e) Illustrative examples of qubit patches (blue) used for simulations of four-site and eight-site models.The ancilla qubits (dotted boxes) are shown in red.The color scheme of the qubits and connectivities represents the "single-qubit randomized benchmarking (RB) average error per gate" and "two-qubit √ iSWAP gate cross-entropy benchmarking (XEB) average error per cycle," respectively.

FIG. 2 .
FIG. 2. (a)The quantum circuit to calculate a finite-temperature dynamic correlation function Û(t) V β using QITE and one ancilla qubit.Measuring X and Y on the ancilla yields the real and imaginary parts of the correlation function, respectively.The thermal average over the initial states is obtained according to Eq. (4).(b) The circuit-recompilation ansatz.Each gate round after the base PhXZ layer includes a layer of two-qubit √ iSWAP † and single-qubit PhXZ gates, as shown in the box.

FIG. 3 .
FIG. 3. (a) The finite-temperature Z i Z j (β = 2) correlation functions computed for the [4Fe-4S] model.The simulation values are in good agreement with the exact results in parentheses.(b) The finite-temperature dynamical correlation function Z 1 (t)Z 2 (0) β (β = 2) and its Fourier transform for the [4Fe-4S] model.(c) The spin-coupling pattern (derived from Z i Z j correlation functions) in the P-cluster and/or Fe-Mo-cofactor model and the simulated finite-temperature Z 1 Z 8 and Z 2 Z 3 for the model, showing the correct magnetic coupling pattern in the low-temperature regime.The deviation between the exact and no-noise results shows the effects of recompilation.(d) The real and imaginary parts of the finite-temperature dynamical correlation function Z 1 (t)Z 2 (0) β (β = 1) and its Fourier transform for the P-cluster and/or Fe-Mo-cofactor model.Error-mitigation acronyms: SE, spin echoes; F, Floquet calibration; RO, readout-error mitigation; PS, postselection; HW, hardware.

FIG. 4 .
FIG. 4. (a) The thermal energy E and heat capacity (c/n = 1 n ∂ E /∂T) for the α-RuCl 3 six-site model (n = 6).The two-peak structure is indicative of the proximate spin-liquid character [29-31].(b) The dynamical correlation function Z 3 (t)Z 4 (0) β and its transform S(ω) at β = 1 for the α-RuCl 3 six-site model.The three plots for the real part show the successive improvement obtained by standard error mitigation (Floquet, readout error, postselection, orange); the introduction of spin echoes (light green); and rescaling (dark green).The same is seen in the plots for the imaginary part and for S(ω), where the results successively including these three classes of techniques are superposed.Error-mitigation acronyms: SE, spin echoes; F, Floquet calibration; RO, readout-error mitigation; PS, postselection; HW, hardware.

Figure 3 (
Figure3shows simulation results on the Fe-S clusters.Despite the simplifications of the model (e.g., to S = 1/2), the static Z i Z j correlation functions at the lowest temperature (β = 2) for the [4Fe-4S] and P and/or Fe-Mo-cofactor clusters reflect the known ground-state pairing patterns of the spins in the true [4Fe-4S] and P clusters [Figs.3(a)and 3(c)].The largest error in the correlation functions is approximately 22% in the P and/or Fe-Mo-cofactor cluster.The majority of this error is from the hardware rather than the classical recompilation; the error from classical recompilation can be seen in the difference between the red and black lines in Fig.3(c).Figure 3(b) shows the dynamical correlation functions and the effects of the different error-mitigation strategies.The total circuit complexity for the [4Fe-4S] cluster dynamical correlation function is 22 two-qubit gates after circuit recompilation of the imaginary-time evolution and real-time evolution blocks and decomposition of the controlled gates [Fig.2(a)] into √ iSWAP † gates; as discussed, the recompilation provides nearly exact results in the classical noiseless emulator, indistinguishable from

TABLE I .
The excitation energy levels of the S = 1/2 and S = 5/2 Heisenberg models for the [4Fe-4S] cluster.The numbers in parentheses indicate the degeneracy.

TABLE II .
The excitation energy levels of the S = 1/2 and S = 5/2 Heisenberg models for the P and/or Fe-Mo-cofactor cluster.

TABLE III .
The ground-state energies of the Kitaev-Heisenberg models.