Chemistry Beyond the Scale of Exact Diagonalization on a Quantum-Centric Supercomputer

• Exact finite-basis solutions grow combinatorially, limiting FCI to about (22e,22o) (1) and (26e,23o) (2); larger systems require approximate methods (3,4).Ground-state energy computation is identified as a common task in theoretical quantum chemistry.Source: Introduction · pp. 1–2
• Pre-fault-tolerant chemistry experiments are limited by impractical energy-estimation runtimes and deep circuits whose gate errors and decoherence weaken signals; circuit costs can scale as M^4 (19).The passage also cites efforts to improve measurement (15–17) and techniques to improve scaling (20).Source: Introduction · p. 2
• SQD is used to study ground-state properties of N2 and the [2Fe-2S] and [4Fe-4S] clusters on 58, 45, and 77 qubits, with up to 3.5 K two-qubit gates.The authors describe the workflow as enabling realistic electronic-structure problems beyond exact diagonalization on pre-fault-tolerant processors.Source: Introduction · p. 2

Results

• The workflow assigns large-circuit execution to quantum processors and distributes classical post-processing, including configuration recovery, projection, and diagonalization, across classical nodes.The method maps molecular spin-orbitals to qubits and measures circuit outputs as bitstrings representing electronic configurations.Source: Results · pp. 2–3, Figure 1 · p. 3

Configuration recovery

• In simulated N2, configuration recovery reduced the quantum signal needed for errors below 10 mEh from approximately 20% to approximately 2%.The self-consistent routine targets configurations with incorrect particle numbers and uses average orbital occupancies to guide bit flips.Source: Configuration recovery · pp. 4–5, Figure 2 · p. 4

Quantum circuits

• For N2 in a 6-31G basis, optimized LUCJ produced a qualitatively correct dissociation curve, whereas restricted CCSD failed under strong static correlation.The experiments used a truncated LUCJ circuit initialized from CCSD parameters without quantum-classical optimization; further optimization could improve sample quality.Source: Figure 3 · p. 5, Quantum circuits · p. 6

Implementation on a quantum computer

• The largest experiment used 77 qubits and distributed [4Fe-4S] projection and diagonalization across 64 Fugaku nodes per batch, with 100 batches.At subspace dimension 100 million, the classical runtime was about 1.5 hours; sampled determinants were extended under spin inversion to approximately restore spin symmetry.Source: Implementation on a quantum computer · pp. 6–7

Triple bond breaking in N2

• Across N2 dissociation, SQD energies stayed within tens of mEh of HCI (31), while restricted CCSD failed to describe the dissociation.The experiment used 58 qubits and 16 million configurations; its combined quantum runtime across bond lengths was approximately 45 minutes.Source: Triple bond breaking in N2 · p. 7, Figure 4 · p. 8

[2Fe-2S] cluster: precision many-body physics

• For [2Fe-2S], energy-variance analysis assessed convergence of low-energy eigenstates as the diagonalization subspace dimension increased.The experiment used 10 batches and approximately 2.46 million measurement outcomes, with a quantum runtime of approximately 45 minutes.Source: Figure 4 · p. 8, [2Fe-2S] cluster: precision many-body… · p. 9

[4Fe-4S] cluster: a stress test for methodology and quantum processors

• At approximately 3.5k two-qubit gates, [4Fe-4S] quantum-derived subspaces had significantly lower energy and variance than uniform-sampled subspaces of equal size, though solutions remained worse than classical methods.SQD energies provide an accuracy metric for comparison with selected configuration interaction and uniform sampling without exact solutions; the sparse-wavefunction approach may also extend to other tasks.Source: Figure 4 · p. 8, [4Fe-4S] cluster: a stress test for… · pp. 9–10
• LUCJ samples produced subspaces of better quality than HCI selection, including dimensions both smaller and larger than the 58 M-dimensional target used for optimization.The LUCJ subspaces agreed well with subspaces labeled “upper bound to optimal.”Source: [4Fe-4S] cluster: a stress test for… · p. 12
• HCI’s more restrictive selection cutoffs converged to the first excited-state energy, while looser cutoffs improved the energy at the cost of much larger subspaces.Removing tens of millions of configurations from the largest HCI wave function did not significantly deteriorate the ground-state approximation, suggesting HCI explores more configurations than necessary.Source: [4Fe-4S] cluster: a stress test for… · p. 12
• The LUCJ circuit was optimized first to match amplitudes from a 58 M-dimensional ground-state estimate, then fine-tuned to minimize SQD energy.The first stage minimized Kullback–Leibler divergence; the second used differential evolution.Source: [4Fe-4S] cluster: a stress test for… · p. 12
• For the [2Fe-2S] active space, HCI could find relevant bitstrings despite exploring and diagonalizing subspaces larger than optimal, showing diagonalization need not be the runtime bottleneck.The passage identifies sub-optimal configuration proposals as the bottleneck in this instance and notes that larger, more strongly correlated systems may challenge classical heuristics, where quantum sampling could offer an advantage.Source: [4Fe-4S] cluster: a stress test for… · p. 12
• Characterizing the applicability of classical and quantum heuristics can guide their combination in quantum-centric supercomputing and shift the search for quantum advantage toward practical problems.The experiments suggest LUCJ circuits can produce better-quality subspaces than the HCI selection heuristic at the same depth and connectivity.Source: [4Fe-4S] cluster: a stress test for… · p. 12

Hamiltonian

• The Born–Oppenheimer Hamiltonian is written in second quantization over NMO orthonormal orbitals, using restricted closed-shell Hartree–Fock molecular orbitals as basis functions.Nuclear positions and atomic numbers, and the number of nuclei, specify the nuclear contribution; atomic units use Bohr for length and Hartree for energy.Source: Hamiltonian · pp. 12–13
• For relativistic or active-space calculations, orbital indices refer to active-space orbitals and Hamiltonian quantities are modified for relativistic effects or inactive-electron density.The passage identifies E0, h_pr, and (pr|qs) as quantities subject to modification.Source: Hamiltonian · p. 13

Molecular species and active spaces

• N2 potential-energy curves used all non-core electrons and orbitals on hardware, with RHF, MP2, CCSD, and CASCI at 6-31G and HCI at cc-pVDZ as classical methods.The calculations enforced D∞h symmetry.Source: Molecular species and active spaces · p. 13
• Table 1 lists the studied molecules, their hardware electron and orbital counts, basis sets, and Hilbert-space dimensions for equal spin populations, without molecular point-group symmetries.Its dimension is D = (NMO choose N/2)^2.Source: Table 1 · p. 13
• RHF, MP2, CCSD, and HCI were used to approximate ground-state energies, with sf-X2C Hamiltonians included for scalar relativistic effects.The RHF, MP2, CCSD, CISD, CASCI, and HCI calculations used PySCF.Source: Molecular species and active spaces · p. 14

Electron configurations and qubit mapping

• The Jordan–Wigner mapping represents each electronic configuration as a computational-basis bitstring using M = 2NMO qubits, with separate spin-up and spin-down bitstring halves.The Fock space contains configurations with all possible orbital filling factors, and the manuscript uses “determinant” and “electronic configuration” interchangeably.Source: Electron configurations and qubit… · p. 14
• A configuration belongs to the right particle sector when each spin-sector bitstring’s Hamming weight equals the exact ground-state electron count for that spin.The RHF bitstring satisfies this condition by construction.Source: Hamiltonian · p. 13, Electron configurations and qubit… · p. 14

Eigenstate solver

• Spin-aware sampling facilitates constructing singlet-compatible states, while a penalty with λ = 0.2 mitigates—but does not eliminate—spin contamination.The sampling procedure facilitates total-spin conservation but does not enforce it; the solver penalizes contributions with the wrong S^2 eigenvalue.Source: Eigenstate solver · p. 15

Sample-based Quantum Diagonalization

• SQD projects sampled configurations into subspaces, diagonalizes the projected Hamiltonian, and uses its ground state as an approximate ground state.Samples are drawn in batches from configurations with the correct particle number, including configurations recovered from noisy samples.Source: Sample-based Quantum Diagonalization · p. 14, Eigenstate solver · p. 16

Self-consistent configuration recovery

• Configuration recovery probabilistically flips bits in noisy samples to restore each spin species’ particle number, using orbital occupancies refined through self-consistent iterations.Samples already in the correct particle sector are left unchanged; approximate eigenstates update the occupancy estimates that guide subsequent corrections.Source: Self-consistent configuration recovery · pp. 16–17

Energy-variance extrapolation

• When approximate states are sufficiently close to an eigenstate, energy-variance points vary linearly, and a least-squares intercept estimates the exact energy (45,67–71).SQD points come from different configuration counts and sample batches; the analysis may also reveal multiple eigenstates near the ground state.Source: Energy-variance extrapolation · pp. 17–18

Quantum circuits

• The local LUCJ approximation supports depth O(1) Jastrow circuits with O(|S| + |S′|) adjacent number-number gates, balancing hardware friendliness and accuracy.Orbital rotations use O(NMO) depth; the ansatz’s accuracy is connected to coupled cluster theory and adiabatic state preparation.Source: Quantum circuits · p. 18
• For dynamically correlated species with restricted closed-shell CCSD parameters, adding exp(K2) to the truncated LUCJ circuit reduces excessive concentration around the Hartree–Fock configuration.The truncated two-layer circuit has the same depth as a single-layer LUCJ circuit.Source: Quantum circuits · p. 19

Initialization of LUCJ parameters

• Initialization uses restricted closed-shell CCSD amplitudes, orders t2 eigenvectors by absolute value, retains the first L matrices, sparsifies non-adjacent gates, and uses the resulting circuit without further optimization.The t2 tensor is reshaped and diagonalized before selecting the matrices and sparsifying the J parameters for the heavy-hex layout.Source: Initialization of LUCJ parameters · pp. 19–20

Mapping to heavy-hex processors

• For systems with NMO > 21 on processors with up to 133 qubits, the mapping limits opposite-spin auxiliary couplings to six orbitals and can discard large interaction terms, reducing accuracy.This occurs when the processor layout cannot connect the spin-orbital segments through enough auxiliary qubits.Source: Mapping to heavy-hex processors · p. 20
• A permutation can place the largest opposite-spin interactions at heavy-hex-compatible orbital positions before sparsification, retaining dominant terms with minimal SWAP overhead.This addresses the potential loss of large interaction terms when the standard mapping cannot connect enough orbitals through auxiliary qubits.Source: Mapping to heavy-hex processors · pp. 20–21

Supplementary materials

• The supplement includes Figs. S1 to S22, Tables S1 to S2, and References (7-106).Source: Supplementary materials · p. 30

Supercomputer

• The listed affiliations include IBM Quantum and its research labs, the University of Colorado Boulder, and RIKEN centers and programs.Source: Supercomputer · p. 31

This PDF file includes:

• The PDF includes Supplementary Text, Figures S1–S22, and Tables S1–S2.Source: This PDF file includes · p. 31

1.1.1 Sample complexity

• Under sufficient concentration, the estimator’s energy converges exponentially with sample count because the probability of missing any of the top m configurations decreases exponentially.The stated failure bound simplifies to p_fail ≤ e^(-N_sβ_m/2) when m ≤ e^(N_sβ_m/2).Source: §1.1.1 · pp. 32–33
• Polynomially many samples in qubit count suffice for inverse-power energy accuracy when configuration probabilities decay exponentially or algebraically.For exponential decay, N_s = Ω(M^2c) guarantees error O(∥Ĥ∥/M^c); for algebraic decay, N_s = O(M^cγ) with m ∝ M^c. These bounds assume the ground-state distribution is sufficiently concentrated.Source: §1.1.1 · p. 34
• Subsampling preserves the chance of capturing all m configurations with exponentially decreasing failure in K, while enabling parallel diagonalization of K matrices of dimension d instead of one N_s matrix.For exponentially decaying probabilities, polynomially many subsamples suffice for constant failure probability; choosing a larger N_s/d ratio reduces matrix dimension.Source: §1.1.1 · p. 34

1.1.2 Assessing wavefunction concentration

• Concentration can be assessed by energy convergence as basis size grows, variability across sample batches, and observables near extremal values, but these tests are not definitive.If energy has not converged at the largest computationally feasible basis size, the wavefunction support exceeds that size.Source: §1.1.2 · p. 35

1.1.3 A numerical study of the accuracy: Hubbard model

• Lower entropy yields lower relative ground-state energy error in the Hubbard study, with better accuracy at U = 16 than at U = 1 for d = 10^3 configurations.For U = 1, the selected configurations cannot represent the many relevant amplitudes; at U = 16, the SQD amplitudes better match the exact amplitudes.Source: §1.1.3 · pp. 36–37, Figure S1 · p. 36

1.2 An analytical lower bound to the probability of recovering families of configurations

• The recovery probability is lower-bounded by α + (1 − α)F(M, N, b, ε), while hardware recovery fractions exceed that bound by orders of magnitude and show no strong size dependence.The analytical bound decays exponentially with qubit count, but recovery is more likely for configurations near the reference bitstring; the bound approaches one as α approaches one.Source: §1.2 · p. 39, Figure S2 · p. 40

2 Optimization of the circuit parameters

• The proposed circuit-parameter formalism minimizes subspace energy averaged over sampled configuration batches, using Monte Carlo estimation instead of exhaustive double-combinatorial evaluation.The estimator uses K ≪ Kmax batches; Kmax can be double-combinatorial in the number of spin-orbitals and electrons.Source: §2 · p. 41
• With K = 1, the optimization reduces to stochastic gradient descent; gradient estimates can otherwise use Monte Carlo batches sampled from the circuit distribution.For many variational circuits, state gradients can be implemented with parameter-shift rules or similar techniques.Source: §2 · p. 42
• Gradient-based optimization has practical limitations: its estimator can be biased when configuration supports differ, and required Hadamard tests use controlled circuits too deep for current processors.The section also proposes gradient-free methods such as COBYLA or simulated annealing, which are used in its numerical study.Source: §2 · p. 42

2.1 Orbital optimization

• Orbital optimization searches for a rotated single-particle basis by jointly optimizing orbital rotations and circuit parameters, while the subspace amplitudes are obtained by diagonalization.The loss averages the Rayleigh quotient over configuration batches, and Monte Carlo estimation handles the double-combinatorial sum.Source: §2.1 · p. 43
• In this work, circuit parameters were not optimized during orbital optimization, so the circuit could not respond to the rotated Hamiltonian and the effect was not maximized.The section says optimization of circuit parameters on quantum experiments is future work, while fixed-parameter circuits already achieve good estimator accuracy.Source: §2 · p. 42, §2.1 · p. 44

3.0.1 Compiled circuits

• The compiled circuits use single-qubit and CNOT gates, with barriers omitted in final compilation; the first orbital rotation prepares the reference Slater determinant directly.This replaces implementing the Bogolyubov transformation on the full state, reducing circuit size without information loss; later orbital rotations do not admit the same simplification.Source: §3.0.1 · pp. 44–45, Figure S3 · p. 45

3.0.2 Effect of the Jastrow terms in the accuracy

• For most N2 bond lengths, SQD has its largest energy error with the determinant circuit, while adding density-density interactions and increasing LUCJ layers generally reduces error.The noiseless study compares LUCJ circuits with L=1 and L=4 against a determinant circuit; circuit optimization takes more iterations for LUCJ than for the determinant circuit.Source: §3.0.2 · pp. 45, 47, Figure S7 · p. 48

3.1.1 processor layouts

• Quantum-run circuits are truncated from two LUCJ layers and use classically obtained CCSD parameters without quantum-run optimization, while remaining classically challenging to simulate.Table S1 reports qubit counts, circuit depths, gate counts, processors, layouts, confidence intervals, and noisy sample counts for the hardware experiments.Source: §3.1.1 · p. 47, Table S1
• N2 runs collected 1 · 10^5 and 9.8304 · 10^4 outcomes per dissociation point at 6-31G and cc-pVDZ, respectively; iron-sulfur runs collected 2.4576 · 10^6 outcomes.The experiments used twirled readout error mitigation and dynamical decoupling, including two X pulses on idle qubits.Source: §3.1.1 · pp. 47–48

3.1.2 Classical resource details

• The [4Fe-4S] workflow uses DICE for distributed SCI, whereas PySCF is used for N2 and [2Fe-2S] systems.The [4Fe-4S] calculation runs on Fugaku and the largest calculation uses 6400 nodes.Source: §3.1.2 · p. 50

4.1 Visualization of the wave function amplitudes

• Ground-state energy errors stayed below 10 mEh across the N2 dissociation curve, and estimated amplitudes matched the largest exact amplitudes best.Agreement in the wavefunction tails deteriorated at larger bond lengths, where the wavefunction is less concentrated and softer orbital-occupancy profiles make configuration recovery less effective.Source: §4.1 · p. 51, Figure S9 · p. 52

4.2 Scaling and runtime

• For [2Fe-2S], eigenstate-solver scaling was optimal at 30+ CPUs per node; for [4Fe-4S], diagonalization speedup from added nodes saturated due to communication overhead.For [4Fe-4S], batch groups were also compared by average minimum energy while holding each diagonalization runtime at 1.44 hours.Source: §4.2 · pp. 52, 54, Figure S11 · p. 53, Figure S12 · p. 54
• For N2, increasing the subspace dimension d lowered estimated energies and reduced unphysical large-bond-length oscillations.The comparison used processor-generated measurement outcomes and the lowest energy among K batches.Source: §4.2 · p. 52, Figure S10 · p. 53

4.4 Orbital occupation numbers for the low-energy spectrum of [2Fe-2S]

• In [2Fe-2S], eigenstates B and C showed greater depletion of molecular orbitals 3 and 4 and greater occupancy of orbitals 16 and 19 than eigenstate A.Localized-orbital occupancies differed most in Fe 3d orbitals; the approximate states broke reflection symmetry, consistent with the importance of Fe 3d antiferromagnetic coupling (41,43).Source: §4.4 · pp. 55, 58, Figure S15 · p. 57

4.3 Effect of orbital optimizations in the dissociation of N2 (cc-pVDZ)

• For N2, orbital optimization had negligible effect near equilibrium but lowered dissociation energies and reduced oscillations for smaller subspaces.Unrestricted orbital optimization can break symmetries, so its use depends on whether those symmetries should be preserved.Source: §4.3 · p. 55, Figure S13 · p. 56

4.5 Quantum signal in the experiments

• Across the 40-, 52-, and 72-qubit experiments, correct-particle-sector fractions were compared with uniform samples, and the results indicated processor-sample signal distinguishable from white noise.The signal assessment also compared SQD energies with and without configuration recovery and compared processor samples against uniform distributions.Source: §4.5 · p. 58

4.5.1 N2 experiments

• For N2, configuration recovery produced more accurate energies and a qualitatively correct dissociation curve, while unrecovered samples yielded energies comparable to RHF.Without recovery, about 1 in 1000 measured configurations had the correct particle number; recovery improved accuracy without requiring more quantum measurements.Source: Figure S16 · p. 59, §4.5.1 · p. 59
• Quantum measurement samples produced more accurate N2 energy estimates than uniform samples from either the right particle sector or full Fock space.The right-sector uniform-sample potential energy surface was nonsmooth and sometimes worse than RHF; the full-Fock comparison favored uniform sampling by using 50·10^6 configurations.Source: §4.5.1 · p. 60, Figure S17 · p. 61

4.5.2 [2Fe-2S] experiments

• With configuration recovery, [2Fe-2S] best energies fell below CCSD, while without it they remained between RHF and CISD.The energy distributions shifted lower as subspace dimension d increased in both cases, but recovery produced lower energies.Source: §4.5.2 · p. 60, Figure S18 · p. 62
• For [2Fe-2S], quantum-sample energy-variance points formed two clusters with extrapolated energy in good agreement with HCI, while right-sector uniform samples did not reveal the lower-energy eigenstate.At the largest d values, the uniform-sample analysis began to form a cluster for the higher-energy eigenstate but did not reveal the second eigenstate.Source: §4.5.2 · p. 61, Figure S19 · p. 63
• At most d = 12.25M configurations let quantum-sample SQD resolve [2Fe-2S] Eigenstate B, versus requiring d > 16M for uniformly sampled Fock-space data.The difference in dimensions was at least 3.75M configurations.Source: §4.5.2 · p. 64, Figure S20 · p. 65