Quantum chemistry simulation of ground- and excited-state properties of the sulfonium cation on a superconducting quantum processor

Abstract

The computational description of correlated electronic structure, and particularly of excited states of many-electron systems, is an anticipated application for quantum devices. An important ramification is to determine the dominant molecular fragmentation pathways in photo-dissociation experiments of light-sensitive compounds, like sulfonium-based photo-acid generators used in photolithography. Here we simulate the static and dynamical electronic structure of the H₃S⁺ molecule, taken as a minimal model of a triply-bonded sulfur cation, on a superconducting quantum processor of the IBM Falcon architecture. To this end, we generalize a qubit reduction technique termed entanglement forging or EF [A. Eddins et al., Phys. Rev. X Quantum, 2022, 3, 010309], currently restricted to the evaluation of ground-state energies, to the treatment of molecular properties. While in a conventional quantum simulation a qubit represents a spin-orbital, within EF a qubit represents a spatial orbital, reducing the number of required qubits by half. We combine the generalized EF with quantum subspace expansion [W. Colless et al., Phys. Rev. X, 2018, 8, 011021], a technique used to project the time-independent Schrodinger equation for ground- and excited-states in a subspace. To enable experimental demonstration of this algorithmic workflow, we deploy a sequence of error-mitigation techniques. We compute dipole structure factors and partial atomic charges along ground- and excited-state potential energy curves, revealing the occurrence of homo- and heterolytic fragmentation. This study is an important step towards the computational description of photo-dissociation on near-term quantum devices, as it can be generalized to other photodissociation processes and naturally extended in different ways to achieve more realistic simulations.

---

1 Introduction

Solving the Schrödinger equation for ground- and excited-states of many-electron quantum systems is one of the grand challenges of contemporary science.^1,2 In particular, the accurate computation of excited-state properties by numerical simulations stands to impact many problems in pure and applied quantum chemistry, exemplified by photochemical processes that result from absorption of photons and promotion of electrons to excited states.

The semiconductor industry has employed these processes to use photolithographic materials in solid-state chip fabrication.^3,4 Fundamentally, fabrication involves coating a silicon wafer with a thin film of a precisely engineered block co-polymer with different functional side chains on blocks that self-assemble into lamellae when processed. In this example, blocks have distinct reaction profiles, and some engineered systems include acid-sensitive side chains that, when reacted, change the block solubility coefficient. A photo-acid generator (PAG) can be embedded in the polymer film and photochemically reacts at specific wavelengths of light to release a free proton in the solid-state that can subsequently react with acid-sensitive side chains.^5–8 Patterns form with the use of a mask that blocks or exposes different parts of the film to ultraviolet (UV) light, modifying the solubility of the polymer such that it can be selectively washed away from the wafer in aqueous solvents. One such effective industrial PAG contains the triphenylsulfonium (Ph₃S⁺) cation.^9–14

The computational description of these photochemical processes poses a number of challenges, for example: characterizing the light–matter interaction to determine transitions to electronic (and, in general, vibronic) excited states; and, in photo-dissociation reactions, assessing the electronic structure of excited states to determine the nature of the dissociation path (e.g. homolytic^15–17 or heterolytic¹⁸). In addition, qualitatively correct and quantitatively accurate calculations require incorporating solvation and thermal effects, and reliably assessing the electronic structure of the studied species by accounting for static and dynamical correlations in realistic basis sets. Capturing such effects accurately is essential for understanding molecular properties, for predictive computations, and ultimately for introducing new PAGs, because the current semiconductor process requires molecular-size order controlling.

Over the last decades, research in computational many-electron quantum mechanics has generated algorithms for conventional classical computers that yield approximate, though often very accurate, estimates of ground- and excited-state molecular properties at polynomial cost.^19–22

Digital quantum computers are an alternative and complementary framework to simulate many-body quantum systems.^23–26 Assuming high-quality qubits in a sufficiently large number, they allow simulation of the time-dependent Schrödinger equation at polynomial cost introducing controllable approximations only,^27,28 thus being capable of accessing a vast class of excited-state properties. Recent advances in hardware manufacturing has produced quantum computers that can carry out computations on a limited scale. Despite the rapid development of quantum hardware, modern quantum computation platforms are immature. As a consequence, simulations of excited states on near-term devices are typically restricted to heuristic quantum subspace algorithms,^29–35 that yield approximations to excited-state wavefunctions and properties within the budget of these devices by projecting the Schrödinger equation onto a suitably constructed subspace. It is therefore a real possibility, and of central importance at this time, to assess the potential usefulness of near-term quantum devices on problems of conceptual and practical interest, e.g. the computation of molecular excited states.

Here, we report the development of a heuristic methodology that leverages structured entanglement in many-electron wavefunctions to calculate ground- and excited-state molecular properties, and its experimental demonstration on a superconducting quantum processor. More specifically, we generalize a qubit reduction technique called entanglement forging (EF),³⁶ initially proposed for variational simulations of ground-state energies, to the computation of generic many-body observables. While in a conventional quantum simulation a qubit represents a spin-orbital, within EF a qubit represents a spatial orbital, reducing the number of required qubits by half.

To improve the accuracy of this technique, and to approximate excited-state energies and properties, we combine EF with quantum subspace expansion (QSE), an example of a heuristic quantum subspace algorithm^29,35,37 which, in its simplest form, projects the Schrödinger equation onto a subspace spanned by single and double excitations on top of a reference wavefunction. The proposed methodology extends the applicability of EF, allowing the computation of a significant set of observables, and that of QSE, facilitating its demonstration on contemporary quantum hardware due to the qubit reduction operated by EF.

We apply the proposed technique, in combination with multiple error mitigation methods, to investigate the gas-phase photo-dissociation of H₃S⁺, taken as the simplest molecular model for Ph₃S⁺. Common to both compounds is the presence of a triply-bonded sulfur cation. The most accurate description of the computational model for Ph₃S⁺ requires the inclusion of the π-conjugated phenyl groups, which determines the energy and character of its excited states, but this feature corresponds to a higher computational cost which makes H₃S⁺ a more suitable target for simulations on contemporary quantum hardware.

We assess the interaction between H₃S⁺ and UV light within the electric dipole approximation in response theory, and characterize dissociation paths as homo- or heterolytic by computing partial atomic charges and other properties of the excited-states. Our study contains approximations and limitations, that we endeavor to document. Notwithstanding these limitations, it illustrates that near-term quantum hardware can be effectively used to explore ground- and excited-states of molecules by means of active-space calculations. While active spaces treated in our study are still small, the underlying methodology naturally extends to larger active spaces. Furthermore, our algorithm is amenable to multiple algorithmic improvements and extensions (for example to treat dynamical correlation, solvation effects, and larger chemical systems), that draw a path towards larger and more realistic simulations.

2 Methods

Several authors have shown that the absorption cross-section of electromagnetic radiation by a molecular system can very generally be represented as a Fourier integral.^38–41

Let Ĥ be the unperturbed time-independent molecular Hamiltonian, with eigenstates Ĥ|Φ_A〉 = ε_A|Φ_A〉. If the system, initially at equilibrium at zero temperature, interacts weakly with an external electric field of frequency ω, transitions from the ground state into other quantum states |Φ_A〉 occur if the frequency of the radiation is close to Δε_A0 = ε_A − ε₀. Assuming a field with wavelength much larger than molecular dimensions, the perturbation can be written as V̂(t) = −[small mu, Greek, circumflex]·Ê(t), where [small mu, Greek, circumflex] is the dipole moment operator. According to time-dependent quantum-mechanical perturbation theory, to first order in the perturbation, the rate of transition from the ground state to any excited state is given by Fermi's golden rule^42–48 and is proportional to the dipole dynamical structure factor (DSF),

where μ_A0 = |〈Φ_A|[small mu, Greek, circumflex]|Φ₀〉|² is the transition dipole between Hamiltonian eigenstates Φ₀ and Φ_A. Combined with the excitation energy Δε_A0, it permits evaluation of the oscillator strength μ_A0Δε_A0, which in turn specifies the absorption cross section of electromagnetic radiation.^45–48 It should be noted that the scalar character of the DSF is due to the assumption of an isotropic system, for which any response is independent from the polarization vector of the incident radiation. Furthermore, while an appropriate description of photo-dissociation requires a joint quantum mechanical treatment of electrons and nuclei,^42,43 particularly near conical intersections, here, concerned with valence-electron UV/vis spectroscopy, we focus on vertical electronic transitions.

DSFs of the form shown in eqn (1) are natural targets for quantum computers, where the time-dependent dipole correlation function f(t) = 〈Φ₀|[small mu, Greek, circumflex](t)·[small mu, Greek, circumflex]|Φ₀〉 can be computed by simulating Hamiltonian time evolution.²⁶ Note that f(t) needs to be computed over a sufficiently long time interval to allow for accurate reconstruction of its Fourier transform, time evolution needs to be controllably approximated e.g. with product formulae, and computation of correlation functions involves deep quantum circuits comprising the Hadamard test or mid-circuit measurements.⁴⁹ Approximate ground- and excited states can also be computed by means of quantum diagonalization algorithms. These heuristic methods may require shallower quantum circuits, making them compatible with near-term quantum computers. Furthermore, they enable evaluation of eqn (1) directly in the frequency domain; for example, they allow to evaluate the right member of eqn (1) by summing over the set of computed excited states without computing time-dependent correlation functions. Access to excited-state wavefunctions also allows computing densities , where the operator ĉ^†_σ(x)/ĉ_σ(x) creates or destroys an electron with spin σ at position x. From these, one can extract partial atomic charges, which in turn allow characterization of the dissociation of a single SH bond as a homolytic (H₃S⁺ → H₂S⁺ + H) or heterolytic (H₃S⁺ → H₂S + H⁺) process, as illustrated in Fig. 1.

Fig. 1 Molecular fragmentation paths. In the homolytic cleavage of a single SH bond of H₃S⁺, the two electrons in the bond are divided equally between H₂S and H, leading to the formation of two radicals (top, marked by black circles). In the heterolytic cleavage, the two electrons are taken by one part of the bond, with formation of closed-shell products H₂S and H⁺.

2.1 Algorithm

In this Section, we describe the algorithms used in the present work. Additional methodological details are provided in Appendix A–C.

2.1.1 Classical preprocessing. The starting point of this study is the definition of a Hamiltonian operator Ĥ. To elucidate the electronic structure of H₃S⁺ along cleavage of a single SH bond, we performed a set of constrained geometry optimizations on a classical computer and, for each geometry, projected the electronic Hamiltonian onto an active space of 6 spatial electrons corresponding to sulfur 3p and hydrogen 1 s. Although the active-space approximation biases the electronic structure of H₃S⁺, it offers the possibility to benchmark the performance of heuristic quantum algorithms and near-term quantum hardware, while establishing a foundation for scaling to more complex quantum simulations. Additional details are provided in Appendix A.

2.1.2 Ground-state calculations. Having defined a Hamiltonian, we start the search for excited states with a preliminary ground-state calculation. To this end, we resort to the EF technique,^36,50,51 based on the idea of partitioning a register of qubits in two halves, and representing the target wavefunction as , where Û(θ) is a parameterized unitary, λ_k are a set of coefficients, and |x_k〉 are a set of computational basis states. Observables are written as linear combinations of tensor products , and their expectation values are expressed (see Appendix B for a derivation) aswhere is evaluated as

Within EF, one prepares the states |ϕ^p_kl〉 on a quantum processor, measures the matrix elements X_kl for X̂ = Â,B̂ and the expectation values in eqn (2). In this formalism, a qubit represents a spatial orbital rather than a spin-orbital, and thus the number of qubits required for a simulation is reduced by half.

The unitary Û(θ) and states |x_k〉 are Ansätze, and parameters θ are optimized variationally along with coefficients λ_k. Here, we choose x_k ∈ {|111000〉, |110100〉} to highlight entanglement across frontier active-space orbitals, and Û(θ) as a product of 6 “hop-gates” (i.e. number-conserving functionally complete 2-qubit gates). The circuits executed in this work are shown in Fig. 2a–e, and additional details are in Appendix B.

Fig. 2 Qubit layout and quantum circuits. (a) Active-space orbitals of H₃S⁺ at equilibrium geometry, obtained from a mean-field simulation of the Born–Oppenheimer Hamiltonian, and mapped onto qubits as illustrated. (b) The circuit diagram depicts an entanglement forging Ansatz, comprising an initial state preparation (light red box) followed by 3 layers of parametrized 2-qubit hop-gates in a brickwall pattern (light blue box) and a final measurement (light orange meter symbols). (c) the red box marked V in panel (b) completes the initialization of 6 qubits in a computational basis state x_k ∈ {|111000〉, |110100〉} (top, middle circuits marked x₀, x₁ respectively) or in a superposition state |ϕ^p₀₁〉 = (|x₀〉 + i^p|x₁〉)/2 (bottom circuit marked ϕ^p₀₁, where a = ⌊p/2⌋ = 0, 0, 1, 1 and b = p%2 = 0, 1, 0, 1 for p = 0, 1, 2, 3 respectively). (d) compilation of a 2-qubit hop-gate into single-qubit and cX gates. (e) depiction of the final measurement operations. (f) depiction of two 6-qubit lines on sub-grids of the ibm_kolkata device. Each circuit involves up to 19 cX gates, 42 single-qubit gates, and 8 variational parameters, and only requires gates between pairs of qubits adjacent in a linear topology.

2.1.3 Excited-state calculations. To access excited states, we extend EF to encompass the framework of quantum diagonalization algorithms, exemplified by QSE.²⁹ Within QSE, a set of excitation operators {Ê_μ}_μ are chosen, Hamiltonian and metric matrices are constructed as H_μν = 〈Ψ_θ|Ê^†_μĤÊ_ν|Ψ_θ〉 and M_μν = 〈Ψ_θ|Ê^†_μÊ_ν|Ψ_θ〉 respectively, and Hamiltonian eigenstates are approximated as , where the columns of c are solutions of the eigenvalue equation Hc_A = Mc_Aε_A. Here, we employ as excitation operators single- and double-electronic excitations, as this choice is natural for electronic systems and compatible with the representation of QSE matrices by eqn (3).

2.1.4 Evaluation of observables. The workflow outlined here allows access to ground-, excited-state, and transition matrix elements of a vast class of operators. Along with the Hamiltonian, we compute the electron number and total spin operators, respectively N̂ and Ŝ², two important constants of motion used to label excited states and classify transitions, and the density matrices (RDMs) . Transition RDMs ρ^A0 provide access to the dipole DSF eqn (1), whereas ground- and excited-state RDMs (respectively ρ⁰⁰ and ρ^AA) provide access to partial atomic charges.

2.2 Hardware experiments

Simulations were run on IBM's 27-qubit processor ibm_kolkata based on the Falcon architecture, as shown in Fig. 2f. Segments of best-performing qubits were selected monitoring average readout and cX errors, and IBM's Qiskit and runtime libraries were used to interface with quantum hardware.⁵² Along with hardware simulations, we performed noiseless and noisy simulations of quantum circuits using the statevector and qasm simulators of Qiskit, and exact diagonalization (full configuration interaction of FCI) calculations using PySCF.^53,54 To reduce decoherence effects and systematic errors occurring on quantum hardware, we resorted to a combination of error mitigation techniques, detailed below and further discussed in Appendix C.

2.2.1 Readout error mitigation. In general, measurement errors over n qubits satisfy the relation Ap_ideal = p_noisy where p_noisy and p_ideal are vectors of probabilities (the former is returned by the noisy quantum system and the latter contains probabilities in the absence of measurement errors) and A is a 2ⁿ × 2ⁿ complete assignment matrix. To mitigate readout errors in such scenario, one needs to execute 2ⁿ circuits to measure p_noisy and compute A, and solve a square system of 2ⁿ linear equations to compute p_ideal.⁵⁵ However, it is often the case that errors on multiple qubits can be well approximated using at most calibration circuits.⁵⁶ This result holds when: A can be approximated as a tensor product of n matrices of shape 2 × 2 (tensored Ansatz); A is diagonally-dominated, allowing efficient matrix-free solution of the linear system; p_noisy is well-approximated by a vector with sparsity at most , where n_s is the number of accumulated statistical samples (or shots).

2.2.2 Post-selection. Since the states |ϕ^p_kl〉 have 3 electrons, only outcomes of a computational basis measurement corresponding to binary strings with Hamming weight 3 are retained.⁵⁷ Furthermore, |ϕ^p_kl〉 is real-valued for p = 0, 2 and thus expectation values of purely-imaginary Pauli operators on such states vanish.

2.2.3 Clifford-based gate error mitigation. For particular parameter configurations, e.g.θ* = 0, the circuits in Fig. 2b–e are in the Clifford group. The expectation value of a linear combination of Pauli operators over such a circuit can be computed exactly at polynomial cost on a classical computer⁵⁸ and measured on a device, offering a pool of data from which to learn the effect of decoherence on measurement outcomes, and mitigate errors.⁵⁹ We use such data to perform an add and subtract correction, i.e. to compute

X_ideal(θ) ≃ X_hw(θ) + X_ideal(θ*) − X_hw(θ*). (4)

2.2.4 Purification. Any n-qubit density operator can be written aswhere the Bloch vector a⃑ is defined so that [small rho, Greek, circumflex] = [small rho, Greek, circumflex]^† and Tr[[small rho, Greek, circumflex]] = 1, and must be compatible with the condition [small rho, Greek, circumflex] ≥ 0. In particular, since the purity Tr[[small rho, Greek, circumflex]²] of a density operator lies⁶⁰ in the interval [2⁻ⁿ, 1], and Tr[[small rho, Greek, circumflex]²] = 2⁻ⁿ(1 + ‖a⃑‖²), then ‖a⃑‖² must lie in the interval [0, 2ⁿ − 1]. Due to decoherence and artifacts of error mitigation, we may observe Tr[[small rho, Greek, circumflex]²] ∉ [2⁻ⁿ,1] despite the target state being pure. When that happens, we scale the Bloch vector so that Tr[[small rho, Greek, circumflex]²] = 1. In the remainder of this work raw, readout-error mitigated (only technique a), and fully error-mitigated results (all of the four mitigation techniques in subsections 2.2.1 to 2.2.4) will be labeled raw, roem, and em respectively.

3 Results

3.1 Ground-state simulations

Fig. 3 shows simulation results for the ground state of H₃S⁺. Variational EF simulations (panels a and c) are on average ∼40 mE_h above FCI, with a non-parallelity error (defined as npe = max_R|ΔE(R) − 〈ΔE(R)〉| with ΔE(R) = E(R) − E_FCI(R)) of 10.2, 11 ± 3, 20 ± 3 mE_h and a binding energy (defined as E(R_max) − min_RE(R) of 163, 164 ± 4, 162 ± 4 mE_h for statevector, qasm and ibm_kolkata respectively. It should be noted that binding energies are challenging to evaluate on quantum hardware, due to a non-smooth potential energy curve for large R. The combined use of EF and QSE (panels b and d) significantly improves the agreement between variational and FCI energies, and decreases non-parallelity errors to 0.4, 2 ± 6, 6 ± 13 mE_h for statevector, qasm and ibm_kolkata respectively. Binding energies decrease slightly, respectively to 159.0, 163 ± 4, 156 ± 5 mE_h.

Fig. 3 Ground-state energies and partial atomic charges during bond cleavage. (Left) Computed total energy (a and b) and deviation between computed and FCI total energy (c and d) from EF (a and c) and QSE with single and double excitations on top of EF (b and d) using simulators (green lines, orange symbols for statevector and qasm) and quantum hardware (ibm_kolkata, red symbols). em is an abbreviations for error mitigation. (Right) Computed atomic charges on S (e and f) and the departing H (g and h) as a function of bond-length from EF (e and g) and QSE (f and h). Charges are computed with a Mulliken population analysis based on meta-Lowdin atomic orbitals. Partial atomic charges on the remaining H atoms are equal to each other, and to (1 − q_S − q_H)/2.

Fig. 3 also shows ground-state partial atomic charges. Hartree-Fock (SCF) incorrectly predicts (panels g and h) the charge q_H on the departing hydrogen to remain finite as R diverges, i.e. a heterolytic ground-state dissociation path. All other methods predict q_H to vanish as R diverges, i.e. a homolytic ground-state dissociation path. These observations are in line with the difference between experimental⁶¹ gas-phase ionization potentials of H₂S and H being (10.5–13.6) eV = −3.1 eV†.

EF inaccurately approximates the electronic structure of the H₂S⁺ moiety, leading to discrepancies between computed and exact values of q_S (panel e). QSE improves agreement of q_S with FCI, but quantitatively significant differences remain (∼0.1 a.u., panel f). The qualitative agreement between computed and exact charges is primarily due to methodological approximations, with decoherence on quantum hardware introducing additional deviations in the amount of ∼0.01 a.u. on average.

3.2 Excited-state simulations

Dipole DSFs are shown in Fig. 4. While positions and strengths of dominant peaks are in qualitative agreement between exact and simulated results, noisy simulations show uncertainties on excitation energies, which translate into broadening and overlapping of peaks. Notwithstanding such limited precision, simulations consistently indicate the absorption of ultraviolet (UV) radiation by H₃S⁺. For all geometries, pronounced peaks are present at ℏω ≃ 20 eV (0.73 E_h or 0.9 nm), in the high-energy end of UV. At R = 0.757, 1.357 Å the spectrum is supported above 13 eV (0.48 E_h or 95 nm). As R further increases, structures appear at lower energies, specifically ℏω ≃ 4–10 eV (0.15–0.37 E_h).

Fig. 4 Dipole spectral functions. Spectral function of the dipole operator from FCI (blue), statevector (green), qasm with noise model and error mitigation (orange) and hardware (ibm_kolkata, red) at the representative bond-lengths R = 0.757, 1.357, 2.057 and 3.957 Å (left to right). Spectral peaks are plotted with a broadening of 0.2 mHa (5 meV) for FCI and statevector, and a broadening reflecting the uncertainty on the excitation energy for qasm and ibm_kolkata.

To interpret the vertical electronic excitations highlighted by dipole DSFs, in Fig. 5 we show low-lying singlet and triplet excited states (respectively S₀, S₁⋯S₄ and T₁⋯T₄ in ascending order of energy at large bond-length). Simulations on classical and quantum devices predict qualitatively correct curves, with T₁, S₀ degenerate for large R and T₂, S₁ more than 100 mE_h above S₀ across dissociation, albeit with statistical uncertainties of ∼5 mE_h for noisy simulations.

Fig. 5 Excited-state energies and partial atomic charges during bond cleavage. (Panels a and b) Total energies of four low-lying excited singlet (left, panel a) and triplet (right, panel b) states from FCI (blue lines), simulators (green lines, orange symbols for statevector, qasm) and quantum hardware (ibm_kolkata, red symbols). (Panels c and d) Computed partial atomic charges on the departing H as a function of SH bond-length for the singlet (left, panel c) and triplet (right, panel d) excited states, from QSE with singles and doubles on top of the EF wavefunction, using FCI (blue lines), simulators (green lines, orange symbols for statevector, qasm), and quantum hardware (ibm_kolkata, red symbols). Dark and light colors indicate the lower- and higher-energy states in the large bond-length regime R respectively.

Vertical singlet–singlet and singlet–triplet gaps computed from Fig. 5 are listed in Table 1. Noiseless simulations tend to overestimate gaps by 10 and 20 mE_h respectively. Noisy classical and hardware simulations are in line with such trends, but feature statistical uncertainties of 5 mE_h respectively. Note also that FCI predicts S₁, S₂ and T₁, T₂ to be degenerate, but these degeneracies are lifted in QSE, even with the statevector simulation, because the excitation operators are not symmetry-adapted.

Table 1 Singlet–singlet and singlet–triplet gaps^a

Gap	FCI	s.v.	qasm	ibm_kolkata
a Vertical singlet–singlet and singlet–triplet gaps from FCI, and from QSE with singles and doubles on top of the EF wavefunction using statevector (abbreviated s.v.), qasm, and ibm_kolkata. Gaps are computed at the equilibrium geometry R = 1.357 Å and listed in Hartree units.
S1–S0	0.484	0.493	0.494 ± 5	0.492 ± 5
S2–S0	0.484	0.501	0.502 ± 6	0.499 ± 6
T1–S0	0.405	0.417	0.417 ± 5	0.416 ± 5
T2–S0	0.405	0.428	0.429 ± 5	0.426 ± 6

---

Fig. 5 also shows partial atomic charges on the departing H as a function of R. Charges exhibit discontinuous behavior around equilibrium geometry. While FCI and noiseless curves are in agreement with each other, noisy and hardware simulations significantly deviate from FCI and noiseless values, indicating the sensitivity of excited-state partial atomic charges to finite measurement error and decoherence. Nevertheless, for large R, all simulations agree that S₁, T₁ and T₂ lead to homolytic dissociation, whereas S₂ leads to heterolytic dissociation.

Comparison between Fig. 4 and 5 indicates that absorption of UV light at the equilibrium geometry (1.357 Å) causes transitions to low-lying singlet states S₁ and S₂ (indeed, the lowest-energy peaks of S(ω) are located at ℏω = ΔE_S1,S0 and ΔE_S2,S0), in turn suggesting coexistence of both homolytic and heterolytic pathways in the gas-phase dissociation of H₃S⁺, described within an active space.

4 Conclusions

In this work, we took a step towards delivering physically relevant simulations on near-term quantum devices. By integrating the EF technique for qubit reduction with quantum diagonalization algorithms exemplified by QSE, we computed ground- and excited-state wavefunctions and properties of H₃S⁺. Combining these algorithmic developments with a sequence of state-of-the-art error mitigation techniques, we experimentally realized the proposed algorithmic workflow on a superconducting quantum computer. Note that this is not a direct extension of previous work, but a careful combination and generalization of independent algorithms.

By computing dipole spectral functions and excited-state energies and partial atomic charges, we were able to elucidate the mechanism of dissociation of H₃S⁺ upon absorption of UV light.

Our study is among the earliest simulations of excited-state molecular spectra on a quantum processor.^35,62–64 Comparison against exact diagonalization indicates that the proposed methodology is capable to deliver accurate results, at least for the active space sizes currently accessible. Notwithstanding this encouraging result, a number of challenges need to be addressed to provide chemically meaningful results, particularly in connection with the photo-acid generating properties of Ph₃S⁺. This goal can be achieved by integrating additional functionalities in the algorithmic workflow considered here.

This study simulates electrons in an active space derived from a minimal basis. While simulations of this kind provide benchmarks and occasions to illustrate algorithmic workflows, useful quantum simulations require accounting for static and dynamical electron correlation in realistic basis sets; on near-term devices, this can be achieved using N-electron valence perturbation theory or otherwise approximate techniques.^65–67 Industrially relevant photodissociation processes require describing realistic functional groups such as phenyl, which on near-term devices can be done integrating the algorithms presented here with quasi-complete active space^68,69 or fragmentation techniques.^70,71 Simulations are carried out in the gas phase, whereas photo-dissociation reactions may occur in a solvent. Solvation effects can be accounted for using implicit or explicit solvation models.^72,73 Research into these algorithmic extensions and improvements is underway.

Encouragingly, the algorithmic workflow considered in this work appears useful, in conjunction with near-term quantum architectures and in combination with other algorithms, and serves to demonstrate the usefulness of hybrid quantum-classical simulation techniques in the continuing search for physically relevant simulations.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Author contributions

G. O. J. and Y.-y. O. conceived, initiated, and supervised the project. M. M. developed the quantum simulation codebase, ran experiments on emulators and quantum hardware, and wrote the initial version of the manuscript. All the authors analyzed computational and experiment results and participated in the composition of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Appendices

A Additional methodological details

For brevity, the Einstein summation convention is used in most of the equations, while vector expansions, linear combinations of tensor products, and summation over nuclei are explicitly described with the summation symbol.

A.1 Classical pre-processing. For each value of the hydrogen-sulfur bond-length studied in this work, we performed a constrained geometry optimization using density functional theory (DFT) with B3LYP-D3 functional and a def2-qzvpp basis set.

For each geometry, we carried out a restricted closed-shell Hartree-Fock (RHF or SCF) calculation with the quantum chemistry software PySCF at STO-6G level, yielding a set of orthonormal molecular orbitals (MOs) , where |χ_α〉 are atomic orbitals AOs, and a second-quantization Born–Oppenheimer Hamiltonian of the form

where indices prqs = 1⋯m label MOs, σ, τ = ↑, ↓ label spin polarizations. The nucleus–nucleus Coulomb interaction is given bywhere R_α and Z_α are the position and atomic number of nucleus α. The coefficientsspecify the one-electron part of the Hamiltonian and the electron–electron Coulomb interaction respectively. Hartree units are used throughout, the numbers of spin-up and spin-down electrons and nuclei are N_↑, N_↓, and N_nuc respectively.

The Born–Oppenheimer Hamiltonian is restricted to an active space spanned by the 6 higher-energy MOs, corresponding to linear combinations of S[3p] and H[1s] orbitals. To this end, we froze the 7 lower-energy MOs with a standard frozen-core procedure,

where indices ij, pr label core and non-core orbitals respectively.

While larger active spaces are necessary for applications of practical interest, we elected to focus on a 6-orbital active space because S is valence-isoelectronic with O, and S is not playing a hypervalent role as, for example, in SO²⁻₄. In Table 2 we show the binding energy, vertical singlet–singlet, and vertical singlet–triplet gaps of H₃S⁺ using CASCI(6e,6o) and CASCI(8e,7o). As seen, results are quantitatively but not qualitatively affected by the extension of the active space. Similarly, both CASCI calculations predict homolytic dissociation along the ground-state potential energy curve.

Table 2 Active space comparison^a

Quantity	CASCI(6e,6o)	CASCI(8e,7o)
a Binding energies (BE) and vertical singlet–singlet, singlet–triplet gaps from CASCI(6e,6o) and CASCI(8e,7o). Gaps are computed at the equilibrium geometry R = 1.357 Å and listed in Hartree units.
S1–S0	0.484	0.474
S2–S0	0.484	0.474
T1–S0	0.405	0.403
T2–S0	0.405	0.403
BE	0.164	0.157

---

Another limitation of the present study is the use of a minimal STO-6G basis set. In Table 3 we show the binding energy, vertical singlet–singlet, and vertical singlet–triplet gaps of H₃S⁺ using classical coupled-cluster with singles and doubles (CCSD) and equation-of-motion CCSD (EOM-CCSD) for ground- and excited-state calculations respectively, in a minimal STO-6G and a correlation-consistent cc-pVTZ basis. Gaps and binding energies are significantly affected by dynamical correlation, but qualitative features (ordering and degeneracy of the excited states) are preserved.

Table 3 Basis set comparison^a

Quantity	STO-6G	cc-pVTZ
a Binding energies (BE) and vertical singlet–singlet, singlet–triplet gaps from CCSD and EOM-CCSD respectively, in a STO-6G and a cc-pVTZ basis set, with 5 frozen orbitals. Gaps are computed at the equilibrium geometry R = 1.357 Å and listed in Hartree units.
S1–S0	0.475	0.359
S2–S0	0.475	0.359
T1–S0	0.403	0.302
T2–S0	0.403	0.302
BE	0.156	0.178

---

A.2 Representation of operators within EF. In this Section of the Appendix, we discuss how to represent second-quantization operators as qubit operators, and how to evaluate their expectation values within EF.
A.2.1 The Hamiltonian. In this work, we chose to partition molecular spin-orbitals into spin-up and spin-down. A Cholesky decomposition^74–77 of the electron-repulsion integral, (pr|qs) = L^γ_prL^γ_qs, and a reordering of the creation and destruction operators,are used to represent the Hamiltonian eqn (6) asand to separate operators acting on spin-up and spin-down molecular spin-orbitals asIn Jordan–Wigner representation,where  and L̂^γ are the qubit representations of Â_↑/Â_↓ and L̂^γ_↑/L̂^γ_↓ restricted to the first/last m qubits respectively.
A.2.2 Spin-summed one-body operators. A simpler formula holds for generic one-body operators,

X̂ = x_prĉ^†_pσĉ_rσ = X̂_↑ + X̂_↓. (14)

In Jordan–Wigner representation it leaves withwhere B̂ is the qubit representation of X̂_↑/X̂_↓ restricted to the first/last m qubits. eqn (15) holds for particle number x_pr = δ_pr, the spin-summed one-body density matrix, x_pr = δ_pp0δ_rr0 for element (p₀, r₀), and the charge-gauge electronic dipole operator,
A.2.3 Total spin. The total spin operator is

Ŝ2 = Ŝ₋Ŝ₊ + Ŝ_z(Ŝz + 1) (17)

where Ŝ₋ = ĉ^†_p↓ĉ_p↑, Ŝ₊ = Ŝ^†₋, and Ŝ_z = ĉ^†_p↑ĉ_p↑ − ĉ^†_p↓ĉ^↓_p. For a closed-shell wavefunction, Ŝ_z(Ŝ_z + 1) = 0, leaving

Ŝ² = _c^†_p↓ĉ_p↑ĉ^†_q↑ĉ_q↓ = N̂_↓ − ĉ^†_p↑ĉ_q↑ĉ^†_p↓ĉ_q↓ (18)

In Jordan–Wigner representation,where Ĉ is the qubit representation of N̂_↓ restricted to the last m qubits, and Ê_pq is the qubit representation of ĉ^†_p↑ĉ_q↑/ĉ^†_p↓ĉ_q↓ restricted to the first/last m qubits.
A.2.4 QSE operators. In this work, we chose the QSE excitation operators to be single- and double-electron excitations, respectively

Ê_ai,σ = ĉ^†_aσĉ_iσ,Ê_aibj,στ = ĉ^†_aσĉ^†_bτĉ_jτĉ_iσ. (20)

In Jordan–Wigner representation,for singles andfor doubles. In eqn (20) and (21), Ê_ai/Ê_aibj is the qubit representation of ĉ^†_a↑ĉ_i↑/ĉ^†_a↑ĉ^†_b↑ĉ_j↑ĉ_i↑ restricted to the first m qubits.

Given two or more operators X̂, Ŷ of the form and , i.e. compatible with eqn (2), their product can be written as

which is also compatible with eqn (2). The construction in eqn (23) is used to represent QSE operators as linear combinations of tensor products.

A.3 Determinant composition of the EF ansatz. Any wavefunction of (N_α, N_β) electrons in m spatial orbitals can be written as a linear combination of electron configurationswhich are Slater determinants, mapped onto bitstrings by conventional fermion-to-qubit mappings. Here, we elected to highlight the spin-up and spin-down parts of the configuration (respectively x and y), but other representations are possible, e.g. a partition based on groups of spatial orbitals. The EF Ansatz is formally derived from a singular value decomposition ψ_ij = ∑_μσ_μU_iμV_jμ, leading to

Thus, the EF Ansatz can reproduce any fermionic wavefunction, provided that (i) all the singular values λ_μ are retained in the representation eqn (25) and (ii) quantum circuits Û and V̂ such that |u_μ〉 = Û|x_μ〉, |v_μ〉 = V̂|y_μ〉 are available.

In practice, however, the non-zero singular values may be up to , which increases combinatorially with active space size, requiring a truncation. Furthermore, while Û and V̂ can be represented as a product of hop-gates with all-to-all connectivity,³⁶ implementing such circuits on near-term devices is technically challenging, requiring to use a heuristic Ansatz.

Determining the expressive power of the EF Ansatz in the presence of a singular-value truncation and of a heuristic Ansatz for the quantum circuits Û and V̂ is an open problem. However, an insightful limiting case can be immediately identified: when the circuits

are equal to the exponential of an anti-Hermitian one-body operator, the EF Ansatz reduces to the familiar multi-configuration self-consistent field (MCSCF) quantum chemistry method. In such limiting case, the Ansatz bitstrings correspond to a set of electronic configurations (MC), and the subsequent circuits perform an orbital optimization within the active space (SCF). In the present work, we focused on approximations to the ground-state wavefunction for which x_μ = y_μ. This choice corresponds to closed-shell Slater determinants, and cannot reproduce open-shell singlet or triplet linear combinations of Slater determinants. In this work, we prepared excited states with such character using QSE (i.e., applying suitable linear combinations of single and double excitations), but for the purpose of improving the accuracy of EF and allowing state-specific excited-state calculations, one must allow x_μ and y_μ to differ.

B Details of simulations

B.1 Ground-state EF calculations. The expectation value of the Hamiltonian is written introducing eqn (13) in eqn (2). The resulting expectation value is minimized, as a function of all free parameters (hop-gate and orbital-optimization angles and coefficients λ_k), using an in-house code⁷⁸ interfaced with the classical optimization method L-BFGS-B^79,80 and the statevector simulator of Qiskit.

The variational optimization of coefficients λ_k is performed as follows: it is observed that the energy is a second-degree polynomial in the variables λ,

where the Schmidt matrixis introduced. Therefore, for a fixed parameter configuration θ, the energy is minimized when the coefficients λ solve the following Lagrange equations,where a constraint is introduced to ensure normalization of the EF wavefunction. The solution of the Lagrange equations is simply h_ml(θ)λ_l = ελ_m, and the energy is minimized when ε is the lowest eigenvalue of h_ml(θ).

B.1.1 Ansatz design. The Ansatz in Fig. 2 was chosen to attain a balance between chemical realism and adequacy for contemporary quantum hardware. The computational basis states x_k are chosen to highlight entanglement between frontier molecular orbitals, and to ensure that preparation unitaries require linear qubit connectivity only, as per Fig. 2c. Similar considerations are made for the subsequent product of hop-gates. The hop-gate portion of the circuit consists of two blocks of gates, each acting on a sub-group of 3 qubits. Such a condition ensures that one of the reference states is preserved by the action of the EF circuit, Û(θ)|x₀〉 = |x₀〉, so that h₀₀(θ) in eqn (28) takes the formand can thus be computed classically, thereby reducing the effect of decoherence on EF simulations.³⁶
B.1.2 HOMO–LUMO orbital optimization. In this study, MOs were not used as active-space basis functions. Instead, we carried out an orbital optimization^81,82 limited to the HOMO–LUMO subspace, i.e. where R_ml(φ) acts as a SU(2) rotation on the HOMO and LUMO orbitals. The angle φ was optimized along with other variational parameters in a preliminary set of ground-state EF calculations (see Section B.1), then the unitary R was used to transform the integrals in the Born–Oppenheimer Hamiltonian with a standard transformation.
B.1.3 Evaluation of observables. Once the optimal hop-gate angles and orbital-optimization angles were computed, quantum state tomography was executed on the EF circuits in Fig. 2, for the purpose of characterizing decoherence effects, enabling purification error mitigation, and avoiding repeated measurements. More specifically, n = 6 qubits were measured in the 3ⁿ eigenbases of X, Y, and Z Pauli operators using the circuit-runner program from the runtime library of Qiskit, and the 4ⁿ entries of the Bloch vector were computed, along with their statistical uncertainties for noisy simulators, by standard post-processing. Post-selection was conducted on the probability distributions of the Pauli measurements. Readout error mitigation and Clifford error mitigation were conducted on the individual entries of the Bloch vector, and purification on the resulting Bloch vector.

Expectation values of quantities like the energy were computed from the entries of the Bloch vector through standard error propagation. For example, given an operator

and a Bloch vector describing the noiseless, noisy, or hardware simulation of |ϕ^p_kl〉 asone haswhere μ_i, v_i are the mean values and statistical uncertainties over Bloch vector entries a_i(p, k, l). Statistical uncertainties are propagated to observable properties (i.e. EF expectation values and derived quantities) with standard error propagation.

B.2 Ground- and excited-state QSE calculations. In the study of Hamiltonian eigenstates with QSE, we first measured the elements of the metric and total spin matricesas described in Sections A.2 and B.1. By solving the eigenvalue problem S_μνf_νt = M_μνf_νtσ_t, QSE matrices are projected on total spin eigenspaces,and the eigenvalue equation Hg = Mgε is solved within each subspace. In summary, QSE energies and spins were defined introducing c_μA = f_μtg_tA and writingStatistical uncertainties were assigned to ε_A, σ_A starting from the definition eqn (36) as described in the following paragraph.
B.2.1 Statistical uncertainties. Determining eigenvalues and eigenvectors of noisy matrices is a notoriously delicate procedure.^83,84 Since eigenvalues of a matrix where elements are normally distributed are not normally distributed, statistical uncertainties are difficult to estimate and are not simply associated with variances of Gaussian distributions.

In this work, we resorted to a simple numerical protocol to assign indicative statistical uncertainties to measured quantities: (i) we solve the eigenvalue equation Sf = Mfσ and Hg = Mgε using the mean values of the matrices S and M, without assigning statistical uncertainties to solutions c = fg. In this work, we obtained non-singular metric matrices, det(M) ≫ 0, so that no eigenvalue truncation was necessary. (ii) we propagate statistical uncertainties from the matrix elements H_μν, S_μν, M_μν to the numerators and denominators of eqn (36), and to the ratio between these quantities, using standard error propagation. (iii) An identical procedure was used to assign statistical uncertainties to particle numbers, RDMs, and derived quantities. Ground- and excited-state RDMs were rescaled so that their trace was statistically compatible with total particle number.

B.2.2 Evaluation of partial atomic charges. Upon computing a ground- or excited-state RDM (the former with EF, the latter with either EF or QSE), the RDM was transformed from the active-space to the MO basis with a simple unitary transformation ρ → R(φ)ρR^†(φ), then an extended RDM was generated, by padding the MO-basis RDM with contributions from frozen MOs,and the extended RDM was transformed to the AO basis, [small rho, Greek, tilde] → C[small rho, Greek, tilde]C⁻¹. Partial charges were then computed with a Mulliken population analysis based on meta-Lowdin atomic orbitals as implemented in the PySCF package. Statistical uncertainties were assigned to partial atomic charges by drawing n = 100 samples of the RDM [small rho, Greek, tilde] in the AO basis, computing partial atomic charges for each sample, and averaging results with standard statistical operations.

C Details of hardware simulations

Hardware simulations were carried out on ibm_kolkata. Jobs submitted on the hardware consisted of 150 circuits and n_s = 100, 000 shots each. We adapted the circuit-runner program from the runtime library of Qiskit. One of the unique options for the circuit-runner program is the ability to correct for measurement errors (i.e. the roem technique) automatically in the cloud.

C.1 Effect of error mitigation

In Fig. 6 and 7, we illustrate the impact of various error mitigation techniques on the result of this work. In Fig. 6a and e we focus on ground-state energies, also shown in Fig. 3 of the main text. Raw/roem noisy simulations on qasm and roem hardware simulations differ by hundreds of mE_h from noiseless results. Comparison between qasm (roem) and ibm_kolkata (roem) results indicates that noise models underestimate the impact of decoherence on observable properties, and comparison between panels (a) and (e) shows that QSE errors are less pronounced than EF errors, on average.

Fig. 6 Effect of error mitigation on ground-state properties. Deviation between statevector and noiseless qasm (light blue), noisy qasm (purple, yellow and orange for raw, readout error mitigated and fully error mitigated data), quantum hardware (ibm_kolkata, dark and light red for readout error mitigate and fully error mitigated data) values of ground-state energy (a and e), total spin (b and f), particle number (c and g), and atomic charge on the departing hydrogen (d and h), from EF (left, panels a–d) and QSE with singles and doubles on top of the EF wavefunction (right, panels e–h). Charges are computed with a Mulliken population analysis based on meta-Lowdin atomic orbitals. roem and em are abbreviations for readout and full error mitigation.

Fig. 7 Effect of error mitigation on excited-state energies. Deviation between statevector and noiseless qasm (light blue), noisy qasm (purple, yellow and orange for raw, readout error mitigated and fully error mitigated data), quantum hardware (ibm_kolkata, dark and light red for readout error mitigate and fully error mitigated data) values of S₁–S₀, T₁–S₀, S₂–S₀, T₂–S₀ (a, b, c, d) energy differences, from QSE with singles and doubles on top of the EF wavefunction (right). roem and em are abbreviations for readout and full error mitigation.

In Fig. 6b and f we focus on total spins, and observe that decoherence causes a form of singlet–triplet spin contamination, greatly reduced by error mitigation. In Fig. 6c and g we focus on particle number. Owing to post-selection, this quantity is essentially noise-free even at the level of roem data. In Fig. 6d and h we show instead partial charges on the departing H. Results indicate that partial charges. While this level of accuracy is satisfactory for qualitative applications (e.g. determining the distribution of electric charge along dissociation), higher accuracy is needed for quantitative tasks, such as determination of electrostatic properties. Fig. 7 shows the effect of error mitigation technique on excited-state energies. Decoherence tends to underestimate singlet–singlet and singlet–triplet gaps, especially for large R. Error mitigation restores agreement with noiseless results, but excited-state energies and derived gaps feature large statistical uncertainties, as documented in Table 1 of the main text. The reduction of such statistical uncertainties is an important direction of future research.

Acknowledgements

We thank Agata M. Brańczyk and Iskandar Sitdikov for generous help and guidance in understanding and using a computational package implementing the entanglement forging method. We acknowledge Riddhi Gupta for valuable interactions regarding the runtime library. We thank Joseph Latone for access to the Clifford cluster at IBM Almaden Research Center, where classical computations were carried out. We thank Andrew Eddins, Hajime Nakamura, Yukio Kawashima, Paul Nation, Zhendong Li, and Barbara Jones for helpful feedback about the manuscript.

Notes and references

1. T. Helgaker, S. Coriani, P. Jørgensen, K. Kristensen, J. Olsen and K. Ruud, Chem. Rev., 2012, 112, 543–631.
2. R. A. Friesner, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 6648–6653.
3. B. D. Gates, Q. Xu, M. Stewart, D. Ryan, C. G. Willson and G. M. Whitesides, Chem. Rev., 2005, 105, 1171–1196.
4. A. S. Gangnaik, Y. M. Georgiev and J. D. Holmes, Chem. Mat., 2017, 29, 1898–1917.
5. O. Nalamasu, M. Cheng, J. M. Kometani, S. Vaidya, E. Reichmanis and L. F. Thompson, Advances in Resist Technology and Processing VII, 1990, pp. 32–48.
6. R. Fallica and Y. Ekinci, J. Mat. Chem. C, 2018, 6, 7267–7273.
7. C. J. Martin, G. Rapenne, T. Nakashima and T. Kawai, J. Photochem. Photobiol., C, 2018, 34, 41–51.
8. K. Sambath, Z. Wan, Q. Wang, H. Chen and Y. Zhang, Org. Lett., 2020, 22, 1208–1212.
9. N. Ohmori, Y. Nakazono, M. Hata, T. Hoshino and M. Tsuda, J. Phys. Chem. B, 1998, 102, 927–930.
10. J. L. Dektar and N. P. Hacker, J. Org. Chem., 1988, 53, 1833–1835.
11. J. L. Dektar and N. P. Hacker, J. Am. Chem. Soc., 1990, 112, 6004–6015.
12. N. Klikovits, P. Knaack, D. Bomze, I. Krossing and R. Liska, Polym. Chem., 2017, 8, 4414–4421.
13. M. Jin, H. Hong, J. Xie, J.-P. Malval, A. Spangenberg, O. Soppera, D. Wan, H. Pu, D.-L. Versace and T. Leclerc, et al. , Polym. Chem., 2014, 5, 4747–4755.
14. W. Zhou, S. M. Kuebler, D. Carrig, J. W. Perry and S. R. Marder, J. Am. Chem. Soc., 2002, 124, 1897–1901.
15. J. W. Knapczyk and W. E. McEwen, J. Am. Chem. Soc., 1969, 91, 145–150.
16. J. Crivello and J. Lam, J. Polym. Sci., Polym. Chem. Ed., 1979, 17, 1059–1065.
17. S. P. Pappas, B. C. Pappas, L. R. Gatechair, J. H. Jilek and W. Schnabel, Polym. Photochem., 1984, 5, 1–22.
18. R. Davidson and J. Goodin, Eur. Polym. J., 1982, 18, 589–595.
19. J. P. F. LeBlanc, et al. , Phys. Rev. X, 2015, 5, 041041.
20. B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang and G. K.-L. Chan, Science, 2017, 358, 1155–1160.
21. M. Motta, D. M. Ceperley, G. K.-L. Chan, J. A. Gomez, E. Gull, S. Guo, C. A. Jiménez-Hoyos, T. N. Lan, J. Li, F. Ma, A. J. Millis, N. V. Prokof'ev, U. Ray, G. E. Scuseria, S. Sorella, E. M. Stoudenmire, Q. Sun, I. S. Tupitsyn, S. R. White, D. Zgid and S. Zhang, Phys. Rev. X, 2017, 7, 031059.
22. K. T. Williams, Y. Yao, J. Li, L. Chen, H. Shi, M. Motta, C. Niu, U. Ray, S. Guo and R. J. Anderson, et al. , Phys. Rev. X, 2020, 10, 011041.
23. I. M. Georgescu, S. Ashhab and F. Nori, Rev. Mod. Phys., 2014, 86, 153.
24. Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre and N. P. Sawaya, et al. , Chem. Rev., 2019, 119, 10856–10915.
25. B. Bauer, S. Bravyi, M. Motta and G. Kin-Lic Chan, Chem. Rev., 2020, 120, 12685–12717.
26. M. Motta and J. E. Rice, Wiley Interdiscip. Rev. Comput. Mol. Sci., 2021, e1580.
27. S. Lloyd, Science, 1996, 273, 1073–1078.
28. J. M. Martyn, Z. M. Rossi, A. K. Tan and I. L. Chuang, PRX Quantum, 2021, 2, 040203.
29. J. R. McClean, M. E. Kimchi-Schwartz, J. Carter and W. A. de Jong, Phys. Rev. A, 2017, 95, 042308.
30. T. Takeshita, N. C. Rubin, Z. Jiang, E. Lee, R. Babbush and J. R. McClean, Phys. Rev. X, 2020, 10, 011004.
31. J. Cohn, M. Motta and R. M. Parrish, PRX Quantum, 2021, 2, 040352.
32. N. Yoshioka, H. Hakoshima, Y. Matsuzaki, Y. Tokunaga, Y. Suzuki and S. Endo, Phys. Rev. Lett., 2022, 129, 020502.
33. E. N. Epperly, L. Lin and Y. Nakatsukasa, A theory of quantum subspace diagonalization, arXiv, 2021, preprint, arXiv:2110.07492,  DOI:10.48550/arXiv.2110.07492.
34. U. Baek, D. Hait, J. Shee, O. Leimkuhler, W. J. Huggins, T. F. Stetina, M. Head-Gordon and K. B. WhaleySay NO to optimization: a non-orthogonal quantum eigensolver, arXiv, 2022, preprint, arXiv:2205.09039,  DOI:10.48550/arXiv.2205.09039.
35. J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. E. Kimchi-Schwartz, J. R. McClean, J. Carter, W. A. de Jong and I. Siddiqi, Phys. Rev. X, 2018, 8, 011021.
36. A. Eddins, M. Motta, T. P. Gujarati, S. Bravyi, A. Mezzacapo, C. Hadfield and S. Sheldon, PRX Quantum, 2022, 3, 010309.
37. S. E. Smart and D. A. Mazziotti, Phys. Rev. Lett., 2021, 126, 070504.
38. M. Lax, J. Chem. Phys., 1952, 20, 1752–1760.
39. E. J. Heller, J. Chem. Phys., 1978, 68, 3891–3896.
40. K. C. Kulander and E. J. Heller, J. Chem. Phys., 1978, 69, 2439–2449.
41. B. R. Johnson and J. L. Kinsey, J. Chem. Phys., 1989, 91, 7638–7653.
42. R. Gordon, J. Chem. Phys., 1965, 43, 1307–1312.
43. C. Boulet and D. Robert, J. Chem. Phys., 1982, 77, 4288–4299.
44. A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt and R. J. Schoelkopf, Rev. Mod. Phys., 2010, 82, 1155.
45. V. Vitale, J. Dziedzic, S. M.-M. Dubois, H. Fangohr and C.-K. Skylaris, J. Chem. Theory Comput., 2015, 11, 3321–3332.
46. D. R. Nascimento and A. E. DePrince III, J. Chem. Theory Comput., 2016, 12, 5834–5840.
47. J. J. Goings, P. J. Lestrange and X. Li, Wiley Interdiscip. Rev. Comput. Mol. Sci., 2018, 8, e1341.
48. X. Li, N. Govind, C. Isborn, A. E. DePrince III and K. Lopata, Chem. Rev., 2020, 120, 9951–9993.
49. S. Endo, I. Kurata and Y. O. Nakagawa, Phys. Rev. Res., 2020, 2, 033281.
50. S. Bravyi, G. Smith and J. A. Smolin, Phys. Rev. X, 2016, 6, 021043.
51. P. Huembeli, G. Carleo and A. Mezzacapo, Entanglement forging with generative neural network models, arXiv, 2022, preprint, arXiv:2205.00933,  DOI:10.48550/arXiv.2205.00933.
52. G. Aleksandrowicz, T. Alexander, P. Barkoutsos, L. Bello, Y. Ben-Haim, D. Bucher, F. Cabrera-Hernández, J. Carballo-Franquis, A. Chen, C. Chen, et al., Qiskit: An open-source framework for quantum computing, 2019, https://zenodo.org/record/2562111#.XhA8qi2ZPyI.
53. Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova and S. Sharma, et al. , Wiley Interdiscip. Rev. Comput. Mol. Sci., 2018, 8, e1340.
54. Q. Sun, et al. , J. Chem. Phys., 2020, 153, 024109.
55. F. B. Maciejewski, Z. Zimborás and M. Oszmaniec, Quantum, 2020, 4, 257.
56. P. D. Nation, H. Kang, N. Sundaresan and J. M. Gambetta, PRX Quantum, 2021, 2, 040326.
57. W. J. Huggins, J. R. McClean, N. C. Rubin, Z. Jiang, N. Wiebe, K. B. Whaley and R. Babbush, Npj Quantum Inf., 2021, 7, 1–9.
58. D. Gottesman, The Heisenberg representation of quantum computers, arXiv, 1998, preprint, arXiv:quant-ph/9807006,  DOI:10.48550/arXiv.quant-ph/9807006.
59. P. Czarnik, A. Arrasmith, P. J. Coles and L. Cincio, Quantum, 2021, 5, 592.
60. M. A. Nielsen and I. L. Chuang, Quantum computation and quantum information, Cambridge University Press, 2010.
61. R. D. Johnson III, et al., NIST 101. Computational chemistry comparison and benchmark database, 1999, https://cccbdb.nist.gov.
62. Q. Gao, G. O. Jones, M. Motta, M. Sugawara, H. C. Watanabe, T. Kobayashi, E. Watanabe, Y.-y. Ohnishi, H. Nakamura and N. Yamamoto, npj Comput. Mater., 2021, 7, 1–9.
63. K. Huang, X. Cai, H. Li, Z.-Y. Ge, R. Hou, H. Li, T. Liu, Y. Shi, C. Chen and D. Zheng, et al. , J. Phys. Chem. Lett., 2022, 13, 9114–9121.
64. I. Khan, M. Tudorovskaya, J. Kirsopp, D. M. Ramo, P. Warrier, D. Papanastasiou and R. Singh, Chemically aware unitary coupled cluster with ab initio calculations on system model H1: a refrigerant chemicals application, arXiv, 2022, preprint, arXiv:2210.14834,  DOI:10.48550/arXiv.2210.14834.
65. A. Tammaro, D. E. Galli, J. E. Rice and M. Motta, J. Phys. Chem. A, 2023, 127, 817–827.
66. T. Takeshita, N. C. Rubin, Z. Jiang, E. Lee, R. Babbush and J. R. McClean, Phys. Rev. X, 2020, 10, 011004.
67. J.-N. Boyn, A. O. Lykhin, S. E. Smart, L. Gagliardi and D. A. Mazziotti, J. Chem. Phys., 2021, 155, 244106.
68. H. Nakano and K. Hirao, Chem. Phys. Lett., 2000, 317, 90–96.
69. H. Nakano, J. Nakatani and K. Hirao, J. Chem. Phys., 2001, 114, 1133–1141.
70. Y. Kawashima, E. Lloyd, M. P. Coons, Y. Nam, S. Matsuura, A. J. Garza, S. Johri, L. Huntington, V. Senicourt and A. O. Maksymov, et al. , Commun. Phys., 2021, 4, 1–9.
71. H. Ma, M. Govoni and G. Galli, npj Comput. Mater., 2020, 6, 1–8.
72. H.-P. Cheng, E. Deumens, J. K. Freericks, C. Li and B. A. Sanders, Front. Chem., 2020, 8, 587143.
73. D. Castaldo, S. Jahangiri, A. Delgado and S. Corni, Quantum simulation of molecules in solution, arXiv, 2021, preprint, arXiv:2111.13458,  DOI:10.48550/arXiv.2111.13458.
74. N. H. F. Beebe and J. Linderberg, Int. J. Quantum Chem., 1977, 12, 683–705.
75. F. Aquilante, T. B. Pedersen and R. Lindh, J. Chem. Phys., 2007, 126, 194106.
76. M. Motta and S. Zhang, Wiley Interdiscip. Rev. Comput. Mol. Sci., 2018, 8, e1364.
77. M. Motta, J. Shee, S. Zhang and G. K.-L. Chan, J. Chem. Theory Comput., 2019, 15, 3510–3521.
78. L. Bello, A. M. Brańczyk, S. Bravyi, A. Eddins, J. Gacon, T. P. Gujarati, I. Hamamura, T. Imamichi, C. Johnson, I. Liepuoniute, M. Motta, M. Rossmannek, T. L. Scholten, I. Sitdikov and S. Woerner, Entanglement forging module, https://github.com/qiskit-community/prototype-entanglement-forging, 2021.
79. C. Zhu, R. H. Byrd, P. Lu and J. Nocedal, ACM Trans. Math. Softw., 1997, 23, 550–560.
80. J. L. Morales and J. Nocedal, ACM Trans. Math. Softw., 2011, 38, 1–4.
81. W. Mizukami, K. Mitarai, Y. O. Nakagawa, T. Yamamoto, T. Yan and Y.-y. Ohnishi, Phys. Rev. Res., 2020, 2, 033421.
82. I. O. Sokolov, P. K. Barkoutsos, P. J. Ollitrault, D. Greenberg, J. Rice, M. Pistoia and I. Tavernelli, J. Chem. Phys., 2020, 152, 124107.
83. J. Lee, F. D. Malone, M. A. Morales and D. R. Reichman, J. Chem. Theory Comput., 2021, 17, 3372–3387.
84. N. S. Blunt, A. Alavi and G. H. Booth, Phys. Rev. B, 2018, 98, 085118.

---

Footnote

† Since E(H₂S⁺ + H) − E(H₂S + H⁺) = [E(H₂S⁺) − E(H₂S)] + [E(H) − E(H⁺)] = IP(H₂S) – IP(H), the inequality IP(H₂S) − IP(H) < 0 implies E(H₂S⁺ + H) < E(H₂S + H⁺).

---