A quantum algorithm for spin chemistry: a Bayesian exchange coupling parameter calculator with broken-symmetry wave functions

Abstract

The Heisenberg exchange coupling parameter J (H = −2JS_i · S_j) characterises the isotropic magnetic interaction between unpaired electrons, and it is one of the most important spin Hamiltonian parameters of multi-spin open shell systems. The J value is related to the energy difference between high-spin and low-spin states, and thus computing the energies of individual spin states are necessary to obtain the J values from quantum chemical calculations. Here, we propose a quantum algorithm, B̲ayesian ex̲change coupling parameter calculator with b̲roken-symmetry wave functions (BxB), which is capable of computing the J value directly, without calculating the energies of individual spin states. The BxB algorithm is composed of the quantum simulations of the time evolution of a broken-symmetry wave function under the Hamiltonian with an additional term jS², the wave function overlap estimation with the SWAP test, and Bayesian optimisation of the parameter j. Numerical quantum circuit simulations for H₂ under a covalent bond dissociation, C, O, Si, NH, OH⁺, CH₂, NF, O₂, and triple bond dissociated N₂ molecule revealed that the BxB can compute the J value within 1 kcal mol⁻¹ of errors with less computational costs than conventional quantum phase estimation-based approaches.

---

1. Introduction

Quantum computing and quantum information processing (QC/QIP) is one of the most innovative research fields in modern science. A quantum computer uses a quantum bit (qubit) as the minimal unit of information, which can possess not only either |0〉 or |1〉 but also arbitrary superposition of the |0〉 and |1〉 states.¹ By utilising the quantum superposition and entangled quantum states, quantum computers can afford to solve certain problems exponentially faster than the classical counterparts. In 2019, a research team of Google Inc. claimed that they achieved “quantum supremacy”.²

Among the diverse topics in the field of QC/QIP, quantum chemical calculations of atoms and molecules are one of the most anticipated applications in the near future. A quantum algorithm based on a quantum phase estimation (QPE) to solve the full configuration interaction (full-CI) in polynomial time, which can give the variationally best wave function within the basis set being used but the computational cost to solve on classical computers increases exponentially against the system size, was reported in 2005.³ A quantum–classical hybrid algorithm known as a variational quantum eigensolver (VQE) executable on the noisy intermediate-scale quantum (NISQ) devices⁴ was proposed in 2014.^5,6 Since then, many reports on the reduction of computational cost with a speedup by improving the quantum algorithms^7–21 have appeared and relevant experimental demonstrations using various quantum devices^22–30 have been documented.

Despite of the rapid progress of the theory for quantum chemical calculations on quantum computers (QCC-on-QCs), a method to efficiently treat open shell electronic structures is still in the stage of infancy. Open shell systems are ubiquitous in chemistry. For example, organic biradicals can be used as prototypes of molecular spin quantum computer,^31,32 polarising agents in dynamic nuclear polarisation (DNP),^32,33 organic light-emitting materials,^34,35 and so on. Open shell multi-nuclear transition metal complexes are often involved in enzymes as reactive centres.^36,37 Single molecule magnets have been extensively studied as molecular memory devices.³⁸ To disclose their electronic structures, sophisticated ab initio quantum chemical calculations are powerful and essential tools. However, in open shell systems carrying spin-β unpaired electrons, the wave function has a strong multi-configurational character and the Hartree–Fock wave function is not a good approximation of the ground state. We have developed theoretical methods to construct spin symmetry-adapted wave functions,^39–41 to determine spin quantum numbers of arbitrary wave functions,⁴² and to purify the spin contaminated wave functions,⁴³ on quantum computers.

In molecules carrying two or more unpaired electrons, the determination of the ground state spin multiplicity and estimation of the spin state energy gap are important tasks. Experimentally, the spin state energy gap is evaluated e.g., by simulating the temperature dependence of magnetic susceptibility or ESR forbidden transition intensity,^44–46 by assuming a Heisenberg spin Hamiltonian^47–51 given in eqn (1).

Here, i and j run over unpaired electrons, J_ij is an exchange coupling parameter, and S_i is an electron spin operator acting on the i-th unpaired electron. Under this definition two spins prefer being ferromagnetic when J_ij is positive. To obtain the exchange coupling parameter J from quantum chemical calculations, one must calculate the energies of high-spin and low-spin states and subtract to obtain the energy difference. The situation has been the same even if we use a quantum computer. However, as we propose in this paper, a quantum algorithm “B̲ayesian ex̲change coupling parameter calculator with b̲roken-symmetry wave functions (BxB)” is capable of calculating the J value directly, without evaluating the energies of individual spin states. The direct calculation of spin state energy gaps and an associated exchange coupling parameter J is very important in QCC-on-QCs, because (1) the J value is usually in the order of kcal mol⁻¹ or even smaller and therefore chemical precision is necessary to evaluate the value quantitatively, and (2) determining the energy down to the fine order of magnitude on quantum computer is extremely cost demanding. Note that quantum algorithms for the direct estimation of the energy gap by combining Ramsey-type measurement and Rabi oscillation experiments or quantum annealing were proposed recently.^52,53

The key technique used in the BxB quantum algorithm consists in (1) quantum simulations of the time evolution of a broken-symmetry wave function under the Hamiltonian with an additional term jS², (2) a SWAP test for the estimation of the wave function overlap, and (3) Bayesian optimisation of the j parameter. Importantly, the proposed approach is easier to be implemented in the process of quantum computing, compared with conventional QPE-based full-CI, and it can afford to calculate the J value within 1 kcal mol⁻¹ of errors regardless of the magnitude and sign of the exchange coupling parameter J.

This paper is organised as follows. In Section 2, we briefly review the theoretical methods for the calculation of exchange coupling parameter J on classical computers, and quantum chemical calculations on quantum computers. Then, we introduce a new quantum algorithm for the J value calculations. In Section 3, numerical simulation results of the J value calculations for H₂ molecule under a covalent bond dissociation, singlet–triplet energy gaps of C, O, Si, NH, OH⁺, CH₂, NF and O₂ are given. N₂ molecule under triple bond dissociation is also discussed as a representative example of systems with more than two unpaired electrons. Conclusions and future perspectives are given in Section 4.

2. Theory

2.1 Quantum chemical calculations of the exchange coupling parameter J on classical computers

Let us assume a biradical molecule as a chemical entity which has two unpaired electrons. The electronic ground state is either a spin-singlet (S = 0) or spin-triplet (S = 1), depending on the nature of magnetic interaction between unpaired electrons. Here, S is a spin quantum number. The strength of the magnetic interaction between unpaired electrons is characterised by an exchange coupling parameter J, by assuming the Heisenberg spin Hamiltonian given in eqn (1). The eigenvalue of the Heisenberg spin Hamiltonian is 3J/2 and −J/2 for the spin-singlet and triplet states, respectively.

The quantum chemical calculation is a powerful tool to disclose the electronic structures of open shell molecules and to understand their magnetic interactions. In quantum chemical calculations, the J value in two-spin systems can be computed from the singlet–triplet energy gap ΔE_S–T, by using the relationship given in eqn (2).

It should be emphasised that the spin-singlet state considered here has two unpaired electrons and thus open shell. If the lowest singlet state is closed shell, the coupled cluster method like CCSD(T) can give reliable ΔE_S–T,⁵⁴ but in the case of the open shell singlet, the Hartree–Fock wave function is not a good approximation and the conventional coupled cluster method is unsuitable. The two-spin open shell singlet wave function is mainly described by the linear combination of two Slater determinants:

Here, 2, α, β, and 0 specify doubly occupied, singly occupied by a spin-α electron, singly occupied by a spin-β electron and unoccupied orbitals, respectively. The multi-configurational nature of the open shell singlet wave function originates from the symmetry requirement of the spin operator S².^55,56 The wave function is a simultaneous eigenfunction of the electronic Hamiltonian H and the S² operator, but single Slater determinant carrying spin-β unpaired electrons is not an eigenfunction of the S² operator. The configuration interaction (CI) methods or the multi-configurational methodologies such as complete active space self-consistent field (CASSCF) are required to treat open shell singlet wave functions explicitly. In fact, the multi-reference (MR)-CI, MR-CC, and multi-reference perturbation theory (MRPT) with the CASSCF reference wave function are capable of predicting the ΔE_S–T and the J value in a quantitative manner.^57–59 In this framework, the energies of high-spin and low-spin states are calculated separately and subtract to obtain the ΔE_S–T and J value.

The fact that the open shell singlet wave function has a multi-configurational character prevents us from directly treating open shell singlet states by means of single configuration theory such as Hartree–Fock (HF) and density functional theory (DFT). Instead, the J value can be calculated by using Yamaguchi's equation given in eqn (4) in conjunction with the broken-symmetry (BS) wave function |Ψ_BS〉.^60,61 The |Ψ_BS〉 is a spin-mixed wave function, and in two-spin systems it is a linear combination of the spin-triplet wave function with M_s = 0 and open shell singlet wave function as given in eqn (5).

Here, 〈S²〉 is an expectation value of the S² operator. For two-spin systems HS corresponds to the spin-triplet state. The E_HS can be easily calculated at the single configuration theory because the spin-triplet wave function with M_s = 1 is generally well approximated by the single Slater determinant |22⋯2αα00⋯0〉. Thus, the energy calculations of the high-spin and BS states are required to obtain the J value from Yamaguchi's equation.

Apart from the methods above, ΔE_S–T can be calculated from the spin-flip variants of the time-dependent DFT, CI, and equation-of-motion coupled cluster (EOM-CC) methods.^62–65 In these methods the high-spin state is calculated by DFT, HF, and CC, respectively, and the low-spin state is described by the spin flip excitations from the high-spin reference state. In these approaches, ΔE_S–T is calculated as an excitation energy from the high-spin state.

As described here, the calculations of ΔE_S–T and the J value generally require two separate calculations, i.e. the energy calculations of the high-spin and low-spin states.

2.2 Quantum chemical calculations on quantum computers

Here we briefly review the two major theoretical methods for quantum chemical calculations on quantum computers, namely QPE-based full-CI and VQE. We also discuss the conventional approach to calculate the J value on a quantum computer.

A quantum circuit for QPE-based full-CI is given in Fig. 1. In Fig. 1, the horizontal lines specify a qubit or n-qubit quantum register, and squares, circles and vertical lines represent quantum gates. Detailed definitions of the quantum gates and quantum circuits are given in the ESI.†

Fig. 1 A quantum circuit of the quantum phase estimation-based full-CI calculations.

In QPE the time evolution of a wave function |Φ〉 is simulated conditional on the qubits for readout and the relative phase shift caused by the time evolution is extracted by means of inverse quantum Fourier transformation (QFT⁻¹) and following qubit measurements.⁶⁶ The quantum state |Φ〉 is projected onto the eigenstate of Hamiltonian by the measurements and thus the full-CI energy can be obtained. Which electronic state is obtained in QPE is probabilistic, depending on the square overlap |〈Φ|Ψ_i〉|². Here, |Ψ_i〉 is the full-CI wave function of the i-th electronic state. By preparing the |Φ〉 to have sufficiently large overlap with the targeted electronic state, the QPE is capable of giving the full-CI energy of a desired electronic state with high probability. In closed shell singlet molecules around their equilibrium geometry or high-spin molecules having no spin-β unpaired electrons, the Hartree–Fock wave function |Ψ_HF〉 is generally a good approximation of the ground state. In open shell molecules carrying spin-β unpaired electrons, one can use a configuration state function |Ψ_CSF〉 that is an eigenfunction of the S² operator as |Φ〉. The |Ψ_CSF〉 can have a large overlap with the full-CI wave function, and it can be easily prepared on quantum computers.^39,40

In QCC-on-QCs the second-quantised Hamiltonian in eqn (6) is transformed into a qubit Hamiltonian by using Jordan–Wigner transformation (JWT),^3,67 Bravyi–Kitaev transformation (BKT),^68,69 or other methods,^3,68,70 and the wave function is mapped onto a quantum register. In this paper, we used the JWT for wave function mapping. In the JWT, the fermionic creation and annihilation operators a^†_p and a_q are transformed to the tensor products of Pauli operators (Pauli strings) as given in eqn (7) and (8), respectively, and the qubit Hamiltonian is expressed by the linear combination of Pauli strings as in eqn (9) and (10).

Here, h_pq and h_pqrs in eqn (6) denote one- and two-electron integrals, respectively. X_i, Y_i, and Z_i are Pauli operators acting on the i-th qubit. Under the JWT, each qubit represents the occupation number of a particular spin orbital, and the number of qubits required for the wave function mapping is equivalent to the number of spin orbitals included in the calculations.

The qubit Hamiltonian contains many noncommutative Pauli strings and Trotter decomposition is often adopted to construct the quantum circuit corresponding to the time evolution operator U = exp(−iHt). The time evolution operator under the first and second order Trotter decompositions are given in eqn (11) and (12), respectively.

The quantum circuit for the time evolution operator exp(−iwX₁Z₂Z₃X₄t) is illustrated in Fig. S1 in ESI† as an example.

The QPE-based full-CI method is a powerful tool in QCC-on-QCs, but there are several challenging issues for the implementation on real quantum devices. First, the quantum circuit is so deep that only small molecules with minimal basis sets can be currently calculated within decoherence time of qubits. Second, the quantum circuit contains many controlled-R_z gates and they should be executed with high fidelity to obtain the reliable energy. Third, the evolution time t becomes very long when we try to determine the energy to a small order of magnitudes. The k-qubit QPE described in Fig. 1 can compute the full-CI energy with an energy tolerance ΔE = π/2^k−1t for U = exp(−iHt). If we set t = 1 we need k = 12 to achieve ∼1 kcal mol⁻¹ of energy tolerance. In this case, we have to apply the time evolution operator up to U²¹¹ = exp(−2¹¹iH). From these reasons it is believed that quantum error correction code is necessary for the implementation of QPE.

In contrast to QPE, VQE is computationally less demanding, and it is expected to be executable on NISQ devices available in the near future. In VQE, the quantum processing unit repeatedly performs the wave function preparation using the parametrised quantum circuit and following measurements to calculate the energy expectation value in eqn (13). The classical processing unit performs a variational optimisation of the parameter θ and feedbacks the parameter to the quantum processing unit.

E(θ) = 〈Φ₀|U^†(θ)HU(θ)|Φ₀〉 (13)

In eqn (13), the |Ψ_HF〉 is generally used as |Φ₀〉. U(θ) determines the quality of wave function. So far, unitary coupled cluster (UCC)^5,6,71 and heuristic ansatzes^25,26,72 have been mainly used to construct U(θ). In VQE, the energy expectation value E(θ) is calculated through repeated measurements and therefore it includes statistical errors. It is known that approximately O(w_max²/ε²) of measurements are required to determine the energy within a precision ε, where w_max = max(|w_m|) in eqn (9).²⁶ Millions of measurements are required to calculate the energy within 1 kcal mol⁻¹ of tolerance in VQE for molecules containing second-row atoms or heavier.

Noteworthily, to calculate the Heisenberg exchange coupling parameter J by using QPE or VQE, two separate energy calculations for the high-spin and low-spin states are necessary, just as on classical computers. Here, we emphasise that the calculation of the J value precisely is much more difficult on quantum computers than on classical computers, because the J value is generally in the order of kcal mol⁻¹ or even smaller. To discuss the spin state energy gap and J value within 1 kcal mol⁻¹ of tolerance, we have to calculate the energy of individual spin states within 0.5 kcal mol⁻¹ of precision. In both QPE and VQE, the computational cost steeply depends on the requested energy precisions; the length of the time evolution must be doubled in QPE, and the measurement number should be quadrupled in VQE, to make the energy precision ε to be half.

We emphasise that most of chemistry problems including the J value calculation focus on the energy difference or relative energy, rather than the total energy itself. From the viewpoint of computational cost reduction and accuracy improvement, it is very useful if we can calculate the energy difference between two electronic states directly in one calculation, without inspecting the total energy of individual electronic states. In the next section, we show that the direct calculation of the J value is possible by using the BxB quantum algorithm, if the system can be approximated by the two-site Heisenberg spin Hamiltonian.

2.3 A BxB quantum algorithm for the calculation of Heisenberg exchange coupling parameter J

As discussed above, the calculation of the spin state energy gap and Heisenberg exchange coupling parameter J generally requires two separate calculations to evaluate the total energies of high-spin and low-spin states, regardless of whether calculations are performed on classical or quantum computers. However, because quantum computers utilise quantum superposition states as the computational resource, we can calculate the spin state energy gap and J value directly on quantum computers.

Unless otherwise specified we focus on two-spin systems. The theory can be easily extended to the systems carrying more than two unpaired electrons if the system is approximated by the two-site Heisenberg model with effective spins S̃ ≥ 1/2.

According to the discussions in Section 2.1, the ground state of two-spin systems is either spin-singlet (J < 0) or spin-triplet (J > 0) state. The energy difference between the lowest spin-singlet and triplet states ΔE_S–T equals to 2J, as given in eqn (2). Importantly, both the spin-singlet and triplet wave functions are simultaneous eigenfunctions of the Hamiltonian and the spin operator S². The eigenvalue of the S² operator is S(S + 1), where S is a spin quantum number. By using this feature, we can reformulate the problem to calculate the J value to another form.

Let us focus on the shifted Hamiltonian H′ with an additional term jS² as in eqn (14).

H′ = H + jS² (14)

H′ and H have the simultaneous eigenfunctions, because [H, S²] = 0. The eigenvalue of the shifted Hamiltonian H′ is given in eqn (15).

E′ = E + jS(S + 1) (15)

Eqn (15) stands for and . From this relationship we can derive the following important theorem:

Theorem 1. The problem to calculate the energy gap between two states belonging to different spin quantum numbers is equivalent to the problem to find the j parameter in the shifted Hamiltonian H′ = H + jS², where the two spin states have the same eigenvalue of H′. The energy gap is calculated as ΔE_LS–HS = j[S_HS(S_HS + 1) − S_LS(S_LS + 1)].

Fig. 2 illustrates a schematic energy diagram of the shifted Hamiltonian H′ with J > 0. In the case of J < 0, the crossing point of the S = 0 and S = 1 lines appears on the left hand side of the figure. Importantly, theorem 1 holds not only for two-spin systems but also systems carrying more than two unpaired electrons. Furthermore, the following theorem can also be derived.

Fig. 2 An energy diagram of the shifted Hamiltonian H′ = H + jS² with j as a variable. J > 0 assumed.

Theorem 2. If the system can be approximated by the two-site Heisenberg spin Hamiltonian with effective spins S̃ ≥ 1/2 as given in eqn (16), all spin states have the same eigenvalue of H′ if j = J.

Theorem 2 is derived from the fact that the spin state energy gap in Heisenberg spin Hamiltonian of eqn (16) is in eqn (17). In Section 3, we treat the triple bond dissociation of a N₂ molecule (S̃ = 3/2), which involves six unpaired electrons.

ΔE_p–q = E_S=p − E_S=q = J[q(q + 1) − p(p + 1)] (17)

From theorem 1 the high-spin and low-spin states have the same eigenenergy of the shifted Hamiltonian H′ if j = J. Thus, the calculation of the J value can be accomplished by finding the j parameter so as for the broken-symmetry wave function |Ψ_BS〉 in eqn (18) to become an eigenfunction of H′.

|Ψ_BS〉 = c_S=0|Ψ_S=0〉 + c_S=1|Ψ_S=1〉 (18)

If |Ψ_BS〉 is an eigenfunction of H′, the time evolution of |Ψ_BS〉 under the shifted Hamiltonian H′ merely gives rise to a phase shift without changing the structure of |Ψ_BS〉. The deviation of |Ψ_BS〉 from the eigenfunction of H′ can be estimated from the square overlap |〈Ψ_BS|U(j, t)|Ψ_BS〉|², where U(j, t) is a time-evolution operator defined in eqn (19). The square overlap |〈Ψ_BS|U(j, t)|Ψ_BS〉|² can be efficiently evaluated on a quantum computer, by utilising a SWAP test.⁷³

U(j, t) = exp(−iH′t) = exp{−i(H + jS²)t} (19)

The quantum circuit to estimate the square overlap |〈Ψ_BS|U(j, t)|Ψ_BS〉|² on a quantum computer is depicted in Fig. 3.

Fig. 3 A quantum circuit to estimate the square overlap |〈Ψ_BS|U(j, t)|Ψ_BS〉|² by utilising a SWAP test.

The quantum circuit in Fig. 3 consists of two parts: the time evolution of |Ψ_BS〉 under the shifted Hamiltonian H′ and a SWAP test. Here, n is equal to the number of active spin orbitals in our implementation. The controlled-SWAP gate appearing as the third quantum gate from left conditionally interchanges the quantum states of qubits storing |Ψ_BS〉 and U(j, t)|Ψ_BS〉 if the control qubit is in the |1〉 state. The quantum state before measurement is given in eqn (20).

Assuming eqn (5) ( in eqn (18)), the probability to obtain the |0〉 state P(0) is calculated as in eqn (21).

Here, ΔE′ denotes an energy difference between the low-spin and high-spin states under the shifted Hamiltonian H′.

ΔE′ = E_S=0 − (E_S=1 + 2j) (22)

For more general situations in which the broken-symmetry wave function is described by a linear combination of many electronic states, the probability P(0) is given by eqn (23), where c_i is the expansion coefficient of the i-th electronic state in the broken-symmetry wave function.

Eqn (21)–(23) reveal that P(0) oscillates depending on j and t, and it gives the maximum value when j = J. Thus, we can calculate the J value by seeking the j parameter giving the maximum P(0), which can be efficiently done by Bayesian inference, as described below. Note that even though the dynamical electron correlation effect is prominent and contributions from electronic states other than high-spin and low-spin states under study to the broken-symmetry wave function are not negligible, the j value giving the maximum P(0) does not change much, as long as the corresponding expansion coefficient c_i is small.

In the field of quantum computing, Bayesian inference has been used in QPE^74,75 and Hamiltonian learning.^76–78 The Bayesian phase estimation (BPE) is robust to experimental imperfections, noise, and decoherence, and it outperforms iterative QPE algorithms. Quantum Hamiltonian learning is a quantum algorithm to efficiently infer the Hamiltonian of an untrusted quantum simulator using a trusted one.

In the BxB quantum algorithm, we infer the j parameter of U(j, t) giving the maximum P(0) in the quantum circuit depicted in Fig. 3. A posterior distribution P(j|0; t) is calculated as eqn (24) from Bayes' theorem.

Here, P(0|j; t) is a likelihood function given in eqn (23). P(j) denotes a prior distribution of j. We can calculate the likelihood function P(0|j; t) by repeatedly executing the quantum circuit in Fig. 3 with various j values.

The BxB algorithm is executed in the following steps. (1) Gives the prior distribution P(j). In this work, we used a normal distribution with μ = 0 and σ² = 1 as the starting prior distribution, where μ and σ² denote the mean and variance, respectively. (2) The evolution time t of U(j, t) is set to be 1.2/σ² for two-spin systems. For triple bond dissociation of N₂ molecule we used t = 0.4/σ². These evolution time lengths were set so that P(0|j; t) becomes a single maximum when j is taken in the range of μ − σ² to μ + σ². (3) Draws m samples in the range of μ − σ² to μ + σ² with a constant interval, and executes the quantum circuit of Fig. 3R times with a given t and j to calculate P(0|j; t). In this work, we used m = 21 and R = 1000. (4) The likelihood function P(0|j; t) is fitted by a normal distribution and the posterior distribution P(j|0; t) was calculated by using eqn (24). (5) If μ(posterior) < μ(prior) − σ²/2 or μ(posterior) > μ(prior) + σ²/2, return to step (3) with the j value giving the largest P(0) as μ(prior). Note that this and next steps are introduced for the purpose of algorithm stability enhance. (6) Set σ²(posterior) = σ²(prior)/5 if σ²(posterior) is smaller than σ²(prior)/5. (7) Convergence check. If σ²(posterior) is smaller than the threshold, the algorithm returns μ(posterior) as the estimate of the j value, and otherwise returns to step (2) with the calculated posterior distribution as the prior distribution of the next iteration. In this study, the convergence threshold is set to be 0.001 Hartree. It should be noted that the convergence threshold controls the precision of the calculated J value. In case the J value of the molecule under study is small (e.g., in the order of cal mol⁻¹), tighter threshold should be adopted.

3. Results and discussion

In order to demonstrate the BxB quantum algorithm we developed a numerical quantum circuit simulation program using Python with OpenFermion⁷⁹ and Cirq⁸⁰ libraries and performed the numerical simulations for covalent bond dissociation in H₂, triplet ground states of C, O, Si, NH, OH⁺, CH₂, NF, and O₂, and triple bond dissociation in N₂. The simulations were performed five times for each geometry. The computational conditions for the numerical simulations are summarised in the ESI.†

3.1 Covalent bond cleavage in H₂ molecule

To demonstrate the BxB quantum algorithm, we performed the numerical quantum circuit simulations for covalent bond dissociation in H₂ molecule in the range of 1.2 Å to 3.0 Å. Note that |Ψ_BS〉 is converged to the RHF root for R(H–H) = 1.1 Å and shorter. We used the STO-3G basis set and |Ψ_BS〉 is mapped onto four qubits. The 2J = ΔE_S–T value calculated from the full-CI/STO-3G level and that obtained from the BxB quantum algorithm simulations, and difference of the J values from the simulations and full-CI are plotted in Fig. 4.

Fig. 4 The 2J = ΔE_S–T values of H₂ molecule under the covalent bond cleavage calculated from the full-CI/STO-3G level (solid line) and from the quantum circuit simulations of the BxB algorithm (blue circles). Inset: the difference of the J values obtained from the quantum circuit simulations and full-CI/STO-3G.

Fig. 4 clearly shows that the BxB algorithm enables us to very accurately calculate the J value, and noticeably the deviations from the J(full-CI) is less than 0.5 kcal mol⁻¹ at all atom–atom distances. The ΔJ value tends to become large for the shorter H–H distance, which presumably originates from the fact that |Ψ_BS〉 is not well described by the linear combination of the lowest singlet and triplet states, and any contributions from other electronic states become non-negligible. As a result, the likelihood function becomes flatten at shorter bond lengths. Also, for longer bond distances the high-spin and low-spin wave functions are expanded by the same sets of Slater determinants. In this case, Trotter decomposition errors on the high-spin and low-spin states are quite similar and they can be cancelled out when we calculate the energy difference. For shorter bond distances such Trotter error cancellations are not sufficient. Thus, the ΔJ value can be improved by adopting tighter threshold for convergence check and the larger number of Trotter decomposition slices. In fact, by changing the convergence threshold to 0.0001 Hartree and setting the time for single Trotter step to 0.2 atomic unit, we obtained J(BxB) = −33.752 ± 0.001 kcal mol⁻¹ at R(H–H) = 1.5 Å, which corresponds to ΔJ = −0.003 kcal mol⁻¹.

It should be also emphasised that Bayesian optimisation converged very fast: only 5–8 iterations are required to achieve the convergence (see Fig. S2 in the ESI†). The evolution time t in the final Bayesian iteration is about 250–300 atomic unit, which should be compared with t ∼ 8000 a.u. to calculate the ΔE_S–T within the same precision (0.5 kcal mol⁻¹ for the energy difference and hence 0.25 kcal mol⁻¹ for the total energy) by using the conventional QPE. Although the BxB algorithm requires the repeated execution of the quantum circuit to evaluate the likelihood function, the computational cost reduction owing to the shortened evolution time is remarkable.

3.2 Atoms and small molecules with spin-triplet ground states

Next, we studied the ground state spin-triplet atoms and small molecules. The target systems are C, O, Si, NH, OH⁺, CH₂, NF, and O₂. The experimental ΔE_S–T value is available for these systems.^81,82 Note that |Ψ_BS〉 is approximated by the linear combination of the lowest spin-triplet and the second lowest spin-singlet wave functions in CH₂ and O₂, rather than the lowest spin-singlet wave function. Therefore, the J value obtained from the BxB algorithm does not correspond to the half of ΔE_S–T in CH₂ and O₂. To calculate the ΔE_S–T that can be directly compared to the experimental value in CH₂ and O₂ molecules, the initial wave function should be replaced from the broken-symmetry wave function to |Ψ₀〉 as defined in eqn (25) and (26), respectively. Here, 2, α, and 0 stand for doubly occupied, singly occupied by a spin-α electron, and unoccupied orbitals, respectively.

The numerical quantum circuit simulations and quantum chemical calculations were carried out using an active space approximation. For neutral atoms (C, O, and Si) and NH, OH⁺, and CH₂ all valence orbitals are included in the active space. In NF and O₂, the valence σ/σ* and π/π* orbitals are considered to be in the active space. The number of qubits required for the wave function mapping is summarised in Table 1.

Table 1 Theoretical and experimental J values of triplet ground state atoms and small molecules

Systems	Number of qubits^a	J(STO-3G)/kcal mol^−1		J(6-311++G**)/kcal mol^−1		J(Exptl.)/kcal mol^−1
		CAS-CI	BxB	CAS-CI	BxB
a Number of qubits required for wave function mapping, and it corresponds to the number of spin orbitals in the active space. b Ref. 81. c Ref. 82.
C	8	22.76	22.59(1)	18.28	18.14(1)	14.57^b
O	8	29.77	29.60(2)	26.07	25.92(2)	22.68^b
Si	8	15.57	15.43(1)	12.50	12.41(1)	9.00^b
NH	10	27.96	27.54(2)	24.80	24.16(2)	19.5^c
OH^+	10	31.68	31.17(1)	31.27	30.49(2)	25.25^c
CH2	12	32.37	31.96(4)	24.97	24.64(3)
NF	12	22.69	22.57(3)	20.64	20.74(2)	17.15^c
O2	12	9.39	9.12(5)	9.69	9.59(4)

---

Because the experimental J values are available for these systems, we have examined basis set dependence on the J values by using STO-3G and 6-311++G** basis sets. The theoretical and experimental J values of the triplet ground state species under study are summarised in Table 1. In Table 1, the J value obtained from the BxB algorithm is the average value of five numerical simulations. The standard deviation of the J(BxB) is less than 0.05 kcal mol⁻¹. Again, the BxB algorithm reproduced the J(CAS-CI) value within 1 kcal mol⁻¹ of errors. The deviation of the J(BxB) value from the J(CAS-CI) is somewhat large in NH and OH⁺ with 6-311++G** basis set. The Trotter decomposition is majorly responsible for this deviation. In fact, by doubling the Trotter slice N in eqn (12), the ΔJ(BxB − full-CI) becomes −0.20 and −0.26 kcal mol⁻¹ for NH and OH⁺, respectively.

As expected, the J value calculated by using 6-311++G** basis set is closer to the experimental J, compared with the minimal basis set STO-3G.

3.3 Triple bond dissociation in N₂ molecule

Applicability of the BxB algorithm is not limited to the two-spin systems. As stated in theorem 2, the BxB algorithm enables us to infer the J value if the system can be approximated by the two-site Heisenberg spin Hamiltonian given in eqn (16) with effective spins S̃ ≥ 1/2. Here, we demonstrate the J value calculation in the triple bond dissociation in N₂ molecule with the (6e, 6o) active space consisting of valence orbitals participating in the triple bond (σ/σ* and π/π* orbitals).

By using the STO-3G basis set, |Ψ_BS〉 converges to the open shell state with six unpaired electrons, those of two occupy the σ-type orbitals and the rest distribute to the π-type orbitals, at the bond distance R(N–N) = 2.1 Å and longer. |Ψ_BS〉 is approximated by the linear combinations of four spin states as given in eqn (27).

By assuming the Heisenberg spin Hamiltonian with effective spins S̃ = 3/2, the spin-triplet (S = 1), spin-quintet (S = 2), and spin-septet (S = 3) states locate at 2J, 6J and 12J, respectively, above the spin-singlet ground state. In reality, the spin state energy gap is not exactly following this relationship because the exchange interaction between electrons in the σ-type orbitals is different from that in the π-type orbitals. The BxB algorithm calculates an effective J, rather than individual J values.

The energy diagrams of the low-lying electronic states of N₂ molecule at the bond length R(N–N) range from 2.1 to 3.0 Å are depicted in Fig. 5. The BxB algorithm predicted the trend of the spin state energy gap in a quantitative manner. Inset of Fig. 5 illustrates the J values obtained from the BxB algorithm and the CAS-CI/STO-3G calculation. The J value obtained from the BxB algorithm agree with the CAS-CI results within 0.2 kcal mol⁻¹ of tolerance. The deviation from the J(CAS-CI) is smaller for the longer N–N distance, where the Heisenberg spin Hamiltonian with effective spins becomes a good approximation. Even for the systems with six unpaired electrons, Bayesian optimisation converged in about 10 iterations and the evolution time t is shorter than 300 atomic unit (see Fig. S3 in the ESI†).

Fig. 5 The energy diagram of the low-lying electronic states in the triple bond dissociation of N₂. Solid lines specify the CAS-CI/STO-3G values, and circles denote the energy levels calculated from the J(BxB) obtained from the numerical quantum circuit simulations. Inset: the J value obtained from the CAS-CI/STO-3G calculations and from the BxB quantum algorithm simulations. The J(CAS-CI) is an average value of J computed from the relationship ΔE_S–T = 2J, ΔE_S–Q = 6J and ΔE_S−7 = 12J.

4. Conclusions

In this work, we proposed a new quantum algorithm BxB for the calculation of the exchange coupling parameter J in Heisenberg spin Hamiltonian on quantum computers. The BxB algorithm enables us to calculate the J value within 1 kcal mol⁻¹ of energy tolerance, regardless of the magnitude and absolute sign of the J value. The BxB algorithm has several advantages against the conventional J value calculations using QPE: the BxB algorithm directly calculates the spin state energy gap without inspecting the total energy of the individual spin states. The evolution time t to obtain the J value within 1 kcal mol⁻¹ of tolerance is about 300 atomic unit in BxB, which is noticeably shorter than that in QPE (t ∼ 4000 a.u. or even longer). The conventional QPE requires good initial guess wave functions having sufficiently large overlap with the full-CI wave function because the energy of which electronic state is obtained is probabilistic. The BxB algorithm uses a broken-symmetry wave function consisting of one Slater determinant, and it is very easy to prepare on quantum computers. The QPE-based full-CI requires many controlled-R_z operations to calculate the energy, which should be executed on quantum computers with high fidelity. By contrast, no controlled-R_z operation is needed in the quantum simulation of the time evolution of wave functions in the BxB algorithm.

In this study, the numerical quantum circuit simulations were performed under the noise-free model, in which all quantum gates are executed perfectly and qubits possess infinitely long coherence time. In real quantum computing devices, noises and errors arising from inaccurate quantum gate operations and decoherences are inevitable. However, we expect that the BxB algorithm is more robust against these errors compared with the conventional QPE, because in the BxB algorithm the wave functions of high-spin and low-spin states are manipulated simultaneously by utilising the quantum superposition. In the BxB algorithm, randomly occurring errors and noises act similarly on the high-spin and low-spin states, and such errors can be cancelled out to some extent in the calculation of energy gaps. Also, the BxB algorithm is robust against measurement errors, because the BxB algorithm calculates a posterior distribution, rather than a point estimate. It should be also emphasised that the number of iterations required for achieving the convergence does not increase against the system size and the number of basis functions, because the BxB algorithm optimises a single parameter j. In addition, the BxB algorithm is free from the barren plateaus problem that is present in variational quantum eigensolver,⁸³ because the evolution time length t is set to be as P(0) changes substantially in the range of μ − σ² < j < μ + σ².

The BxB algorithm needs (2n + 1) qubits to implement the computation for the system having the number of spin orbitals, n. This requirement hampers the BxB algorithm compared with the iterative QPE, which uses (n + 1) qubits. However, the (n + 1) qubits out of (2n + 1) can be initialised just before the SWAP test, and the number of qubits which are necessary to retain the quantum superposition during the time evolution simulation is the same for the iterative QPE and the BxB algorithm. Also, the S² operator only acts on the singly occupied molecular orbitals (SOMOs) and therefore it is possible to perform the SWAP test using qubits which store the occupation number of SOMOs in |Ψ_BS〉 only. By adopting this approximation the number of qubits can be reduced to (n + 2 × n_spin + 1), where n_spin is the number of unpaired electrons in |Ψ_BS〉. Using this strategy for the calculation of the J value in carbon atom with 6-311++G** basis set, we obtained J(BxB) = 18.17 ± 0.01 kcal mol⁻¹, which is very close to that obtained by the SWAP test with all the qubits (J = 18.14 ± 0.01 kcal mol⁻¹).

The BxB algorithm in the current form is applicable only for the systems whose electronic structure can be well approximated by the two-site Heisenberg spin Hamiltonian. However, this limitation originates from |Ψ_BS〉. Importantly, theorem 1 is general and it always holds. This fact insists that we can calculate the energy gap between two states belonging to different spin multiplicities using the BxB-analogue quantum algorithms, if we can prepare the quantum state in superposition of high-spin and low-spin states. Combinations of the BxB algorithm and adiabatic state preparation (ASP) technique^23,84 can potentially open the door to such applications. Another important extension of theory is applications to the molecules with three or more spin centres like triradicals and tetraradicals. The possibilities of these methodologies will be discussed in the forthcoming papers.

Author contributions

K. Sugisaki, K. Sato, and T. Takui planned and conducted the project. K. Sugisaki proposed the theory and performed numerical simulations. All the other authors discussed the results. K. Sugisaki and T. Takui wrote the paper.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by AOARD Scientific Project on “Molecular Spins for Quantum Technologies” (Grant FA2386-17-1-4040, 4041), USA, and by KAKENHI Scientific Research B (Grant No. 17H03012) and Scientific Research C (Grant No. 18K03465) from JSPS, Japan, and JST PRESTO “Quantum Software” project (Grant No. JPMJPR1914).

Notes and references

1. M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information.Cambridge University Press, Cambridge, 2000.
2. F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven and J. M. Martinis, Quantum supremacy using a programmable superconducting processor, Nature, 2019, 574, 505–510.
3. A. Aspuru-Guzik, A. D. Dutoi, P. J. Love and M. Head-Gordon, Simulated quantum computation of molecular energies, Science, 2005, 309, 1704–1707.
4. J. Preskill, Quantum computing in the NISQ era and beyond, Quantum, 2018, 2, 79.
5. M.-H. Yung, J. Casanova, A. Mezzacapo, J. McClean, L. Lamata, A. Aspuru-Guzik and E. Solano, From transistor to trapped-ion computers for quantum chemistry, Sci. Rep., 2014, 4, 3589.
6. A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik and J. L. O'Brien, A variational eigenvalue solver on a photonic quantum processor, Nat. Commun., 2014, 5, 4213.
7. Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis and A. Aspuru-Guzik, Quantum chemistry in the age of quantum computing, Chem. Rev., 2019, 119, 10856–10915.
8. S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin and X. Yuan, Quantum computational chemistry, Rev. Mod. Phys., 2020, 92, 015003.
9. B. Bauer, S. Bravyi, M. Motta and G. K.-L. Chan, Quantum algorithms for quantum chemistry and quantum materials science, Chem. Rev., 2020, 120, 12685–12717.
10. I. Kassal, S. P. Jordan, P. J. Love, M. Mohseni and A. Aspuru-Guzik, Polynomial-time quantum algorithm for the simulation of chemical dynamics, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 18681–18686.
11. J. D. Whitfield, J. Biamonte and A. Aspuru-Guzik, Simulation of electronic structure Hamiltonians using quantum computers, Mol. Phys., 2011, 109, 735–750.
12. M. B. Hastings, D. Wecker, B. Bauer and M. Troyer, Improving quantum algorithms for quantum chemistry, Quantum Inf. Comput., 2015, 15, 1–21.
13. R. Babbush, N. Wiebe, J. McClean, J. McClain, H. Neven and G. K.-L. Chan, Low-depth quantum simulation of materials, Phys. Rev. X, 2018, 8, 011044.
14. I. D. Kivlichan, J. McClean, N. Wiebe, C. Gidney, A. Aspuru-Guzik, G. K.-L. Chan and R. Babbush, Quantum simulation of electronic structure with linear depth and connectivity, Phys. Rev. Lett., 2018, 120, 110501.
15. M. Ganzhorn, D. J. Egger, P. Barkoutsos, P. Ollitrault, G. Salis, N. Moll, M. Roth, A. Fuhrer, P. Mueller, S. Woerner, I. Tavernelli and S. Filipp, Gate-efficient simulation of molecular eigenstates on a quantum computer, Phys. Rev. Appl., 2019, 11, 044092.
16. R. Babbush, D. W. Berry, J. R. McClean and H. Neven, Quantum simulation of chemistry with sublinear scaling in basis size, npj Quantum Inf., 2019, 5, 92.
17. I. G. Ryabinkin, T.-C. Yen, S. N. Genin and A. F. Izmaylov, Qubit coupled cluster method: a systematic approach to quantum chemistry on a quantum computer, J. Chem. Theory Comput., 2018, 14, 6317–6326.
18. H. R. Grimsley, S. E. Economou, E. Barnes and N. J. Mayhall, An adaptive variational algorithm for exact molecular simulations on a quantum computer, Nat. Commun., 2019, 10, 3007.
19. R. Xia, T. Bian and S. Kais, Electronic structure calculations and the Ising Hamiltonian, J. Phys. Chem. B, 2018, 122, 3384–3395.
20. S. McArdle, T. Jones, S. Endo, Y. Li, S. C. Benjamin and X. Yuan, Variational ansatz-based quantum simulation of imaginary time evolution, npj Quantum Inf., 2019, 5, 75.
21. J. Argüello-Luengo, A. González-Tudela, T. Shi, P. Zoller and J. I. Cirac, Analogue quantum chemistry simulation, Nature, 2019, 574, 215–218.
22. B. P. Lanyon, J. D. Whitfield, G. G. Gillett, M. E. Goggin, M. P. Almeida, I. Kassal, J. D. Biamonte, M. Mohseni, B. J. Powell, M. Barbieri, A. Aspuru-Guzik and A. G. White, Towards quantum chemistry on a quantum computer, Nat. Chem., 2010, 2, 106–111.
23. J. Du, N. Xu, X. Peng, P. Wang, S. Wu and D. Lu, NMR implementation of a molecular hydrogen quantum simulation with adiabatic state preparation, Phys. Rev. Lett., 2010, 104, 030502.
24. P. J. J. O'Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik and J. M. Martinis, Scalable quantum simulation of molecular energies, Phys. Rev. X, 2016, 6, 031007.
25. A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow and J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature, 2017, 549, 242–246.
26. G. Zhu, Y. Subaşı, J. D. Whitfield and M. Hafezi, Hardware-efficient fermionic simulation with a cavity-QED system, npj Quantum Inf., 2018, 4, 16.
27. J. E. Rice, T. P. Gujarati, T. Y. Takeshita, J. Latone, M. Motta, A. Hintennach and J. M. Garcia, Quantum chemistry simulations of dominant products in lithium-sulfur batteries, 2020, arXiv, preprint, arXiv:2001.01120, https://arxiv.org/abs/2001.01120.
28. F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, S. Boixo, M. Broughton, B. B. Buckley, D. A. Buell, B. Burkett, N. Bushnell, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, S. Demura, A. Dunsworth, E. Farhi, A. Fowler, B. Foxen, C. Gidney, M. Giustina, R. Graff, S. Habegger, M. P. Harrigan, A. Ho, S. Hong, T. Huang, W. J. Huggins, L. Ioffe, S. V. Isakov, E. Jeffrey, Z. Jiang, C. Jones, D. Kafri, K. Kechedzhi, J. Kelly, S. Kim, P. V. Klimov, A. Korotkov, F. Kostritsa, D. Landhuis, P. Laptev, M. Lindmark, E. Lucero, O. Martin, J. M. Martinis, J. R. McClean, M. McEwen, A. Megrant, X. Mi, M. Mohseni, W. Mruczkiewicz, J. Mutus, O. Naaman, M. Neeley, C. Neill, H. Neven, M. Y. Niu, T. E. O'Brien, E. Ostby, A. Petukhov, H. Putterman, C. Quintana, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, D. Strain, K. J. Sung, M. Szalay, T. Y. Takeshita, A. Vainsencher, T. White, N. Wiebe, Z. J. Yao, P. Yeh and A. Zalcman, Hartree-Fock on a superconducting qubit quantum computer, Science, 2020, 369, 1084–1089.
29. Y. Wang, F. Dolde, J. Biamonte, R. Babbush, V. Bergholm, S. Yang, I. Jakobi, P. Neumann, A. Aspuru-Guzik, J. D. Whitfield and J. Wrachtrup, Quantum simulation of helium hydride cation in a solid-state spin register, ACS Nano, 2015, 9, 7769–7774.
30. C. Hempel, C. Maier, J. Romero, J. McClean, T. Monz, H. Shen, P. Jurcevic, B. P. Lanyon, P. Love, R. Babbush, A. Aspuru-Guzik, R. Blatt and C. F. Roos, Quantum chemistry calculations on a trapped-ion quantum simulator, Phys. Rev. X, 2018, 8, 031022.
31. S. Nakazawa, S. Nishida, T. Ise, T. Yoshino, N. Mori, R. D. Rahimi, K. Sato, Y. Morita, K. Toyota, D. Shiomi, M. Kitagawa, H. Hara, P. Carl, P. Höfer and T. Takui, A synthetic two-spin quantum bit: g-engineered exchange-coupled biradical designed for controlled-NOT gate operations, Angew. Chem., Int. Ed., 2012, 124, 9998–10002.
32. K. Sato, R. Hirao, I. Timofeev, O. Krumkacheva, E. Zaytseva, O. Rogozhnikova, V. M. Tormyshev, D. Trukhin, E. Bagryanskaya, T. Gutmann, V. Klimavicius, G. Buntkowsky, K. Sugisaki, S. Nakazawa, H. Matsuoka, K. Toyota, D. Shiomi and T. Takui, Trityl-aryl-nitroxide-based genuinely g-engineered biradicals, as studied by dynamic nuclear polarization, multifrequency ESR/ENDOR, arbitrary wave generator pulse microwave waveform spectroscopy, and quantum chemical calculations, J. Phys. Chem. A, 2019, 123, 7507–7517.
33. K.-N. Hu, H. Yu, T. M. Swager and R. G. Griffin, Dynamic nuclear polarization with biradicals, J. Am. Chem. Soc., 2004, 126, 10844–10845.
34. H. Namai, H. Ikeda, Y. Hoshi, N. Kato, Y. Morishita and K. Mizuno, Thermoluminescence and a new organic light-emitting diode (OLED) based on triplet–triplet fluorescence of the trimethylenemethane (TMM) biradical, J. Am. Chem. Soc., 2007, 129, 9032–9036.
35. H. Namai, H. Ikeda, Y. Hoshi and K. Mizuno, Thermoluminescence originating from the singlet excited state of 1,4-diarylcyclohexane-1,4-diyls: a potentially general strategy for the observation of short-lived biradicals, Angew. Chem., Int. Ed., 2007, 46, 7396–7398.
36. J. Yano and V. Yachandra, Mn₄Ca cluster in photosynthesis. Where and how water is oxidized to dioxygen, Chem. Rev., 2014, 114, 4175–4205.
37. Z. Li, J. Li, N. S. Dattani, C. J. Umrigar and G. K.-L. Chan, The electronic complexity of the ground-state of the FeMo cofactor of nitrogenase as relevant to quantum simulations, J. Chem. Phys., 2019, 150, 024302.
38. L. Bogani and W. Wernsdorfer, Molecular spintronics using single-molecule magnets, Nat. Mater., 2008, 7, 179–186.
39. K. Sugisaki, S. Yamamoto, S. Nakazawa, K. Toyota, K. Sato, D. Shiomi and T. Takui, Quantum chemistry on quantum computers: a polynomial-time quantum algorithm for constructing the wave functions of open-shell molecules, J. Phys. Chem. A, 2016, 120, 6459–6466.
40. K. Sugisaki, S. Yamamoto, S. Nakazawa, K. Toyota, K. Sato, D. Shiomi and T. Takui, Open shell electronic state calculations on quantum computers: a quantum circuit for the preparation of configuration state functions based on Serber construction, Chem. Phys. Lett.: X, 2019, 1, 100002.
41. K. Sugisaki, S. Nakazawa, K. Toyota, K. Sato, D. Shiomi and T. Takui, Quantum chemistry on quantum computers: a method for preparation of multiconfigurational wave functions on quantum computers without performing post-Hartree–Fock calculations, ACS Cent. Sci., 2019, 5, 167–175.
42. K. Sugisaki, S. Nakazawa, K. Toyota, K. Sato, D. Shiomi and T. Takui, Quantum chemistry on quantum computers: quantum simulations of the time evolution of wave functions under the S² operator and determination of the spin quantum number S, Phys. Chem. Chem. Phys., 2019, 21, 15356–15361.
43. K. Sugisaki, K. Toyota, K. Sato, D. Shiomi and T. Takui, A probabilistic spin annihilation method for quantum chemical calculations on quantum computers, Phys. Chem. Chem. Phys., 2020, 22, 20990–20994.
44. B. Bleaney and K. D. Bowers, Anomalous paramagnetism of copper acetate, Proc. R. Soc. London, Ser. A, 1952, 214, 451–465.
45. A. Rajca, Organic diradicals and polyradicals: from spin coupling to magnetism?, Chem. Rev., 1994, 94, 871–893.
46. H. Iwamura, N. Koga and K. Matsuda, Studies of the electronic structure of polyradicals by means of their magnetic properties, Pure Appl. Chem., 1998, 70, 1953–1960.
47. W. Heisenberg, Mehrkörperproblem und Resonanz in der Quantenmechanik, Z. Phys., 1926, 38, 411–426.
48. P. A. M. Dirac, On the theory of quantum mechanics, Proc. R. Soc. London, Ser. A, 1926, 112, 661–677.
49. J. H. Van Vleck, A survey of the theory of ferromagnetism, Rev. Mod. Phys., 1945, 17, 27–47.
50. P. A. M. Dirac, Quantum mechanics of many-electron systems, Proc. R. Soc. London, Ser. A, 1929, 123, 714–733.
51. T. Arai, Exchange interaction and Heisenberg's spin Hamiltonian, Phys. Rev., 1962, 126, 471–488.
52. A. E. Russo, K. M. Rudinger, B. C. A. Morrison and A. D. Baczewski, Evaluating energy differences on a quantum computer with robust phase estimation, 2020, arXiv, preprint, arXiv:2007.08697, https://arxiv.org/abs/2007.08697.
53. Y. Matsuzaki, H. Hakoshima, K. Sugisaki, Y. Seki and S. Kawabata, Direct estimation of the energy gap between the ground state and excited state with quantum annealing, 2020, arXiv, preprint, arXiv:2007.10561, https://arxiv.org/abs/2007.10561.
54. B. Hajgató, M. Huzak and M. S. Deleuze, Focal point analysis of the singlet–triplet energy gap of octacene and larger acenes, J. Phys. Chem. A, 2011, 115, 9282–9293.
55. G. W. Pratt Jr, Eigen functions of S² by a spin operator method, Phys. Rev., 1953, 92, 278–288.
56. R. Pauncz, The Construction of Spin Eigenfunctions. An Exercise Book.Kluwer Academic/Plenum Publishers, New York, 2000.
57. G. Chung and D. Lee, Theoretical study of some bis-verdazyl diradicals: singlet–triplet energy gap, Chem. Phys. Lett., 2001, 350, 339–344.
58. M. Winkler, Singlet–triplet energy splitting and excited states of phenylnitrene, J. Phys. Chem. A, 2008, 112, 8649–8653.
59. M. G. Delcey, K. Pierloot, Q. M. Phung, S. Vancoillie, R. Lindh and U. Ryde, Accurate calculations of geometries and singlet–triplet energy differences for active-site models of [NiFe] hydrogenase, Phys. Chem. Chem. Phys., 2014, 16, 7927–7938.
60. K. Yamaguchi, The electronic structures of biradicals in the unrestricted Hartree–Fock approximation, Chem. Phys. Lett., 1975, 33, 330–335.
61. K. Yamaguchi, Y. Takahara and T. Fueno, Ab-initio molecular orbital studies of structure and reactivity of transition metal-oxo compounds, in Applied Quantum Chemistry, ed. V. H. Smith, H. F. Schaefer and K. Morokuma, Springer, Dordrecht, 1986, pp. 155–184.
62. L. V. Slipchenko and A. I. Krylov, Singlet–triplet gaps in diradicals by the spin-flip approach: a benchmark study, J. Chem. Phys., 2002, 117, 4694–4708.
63. Z. Rinkevicius and H. Ågren, Spin-flip time dependent density functional theory for singlet–triplet splittings in σ,σ-biradicals, Chem. Phys. Lett., 2010, 491, 132–135.
64. C. U. Ibeji and D. Ghosh, Singlet–triplet gaps in polyacenes: a delicate balance between dynamic and static correlations investigated by spin-flip methods, Phys. Chem. Chem. Phys., 2015, 17, 9849–9856.
65. N. Orms and A. I. Krylov, Singlet–triplet energy gaps and the degree of diradical character in binuclear copper molecular magnets characterized by spin-flip density functional theory, Phys. Chem. Chem. Phys., 2018, 20, 13127–13144.
66. D. S. Abrams and S. Lloyd, Quantum algorithm providing exponential speed increase for finding eigenvalues and eigenvectors, Phys. Rev. Lett., 1999, 83, 5162–5165.
67. P. Jordan and E. Wigner, Über das Paulische Äquivalenzverbot, Z. Phys., 1928, 47, 631–651.
68. J. T. Seeley, M. J. Richard and P. J. Love, The Bravyi–Kitaev transformation for quantum computation of electronic structure, J. Chem. Phys., 2012, 137, 224109.
69. A. Tranter, S. Sofia, J. Seeley, M. Kaicher, J. McClean, R. Babbush, P. V. Coveney, F. Mintert, F. Wilhelm and P. J. Love, The Bravyi–Kitaev transformation: properties and applications, Int. J. Quantum Chem., 2015, 115, 1431–1441.
70. K. Setia and J. D. Whitfield, Bravyi–Kitaev superfast simulation of electronic structure on a quantum computer, J. Chem. Phys., 2018, 148, 164104.
71. A. G. Taube and R. J. Bartlett, New perspectives on unitary coupled-cluster theory, Int. J. Quantum Chem., 2006, 106, 3393–3401.
72. Y. Li, J. Hu, X.-M. Zhang, Z. Song and M.-H. Yung, Variational quantum simulation for quantum chemistry, Adv. Theory Simul., 2019, 1800182.
73. H. Buhrman, R. Cleve, J. Watrous and R. de Wolf, Quantum fingerprinting, Phys. Rev. Lett., 2001, 87, 167902.
74. N. Wiebe and C. Granade, Efficient Bayesian phase estimation, Phys. Rev. Lett., 2016, 117, 010503.
75. S. Paesani, A. A. Gentile, R. Santagati, J. Wang, N. Wiebe, D. P. Tew, J. L. O'Brien and M. G. Thompson, Experimental Bayesian quantum phase estimation on a silicon photonic chip, Phys. Rev. Lett., 2017, 118, 100503.
76. N. Wiebe, C. Granade, C. Ferrie and D. Cory, Quantum Hamiltonian learning using imperfect quantum resources, Phys. Rev. A: At., Mol., Opt. Phys., 2014, 89, 042314.
77. N. Wiebe, C. Granade, C. Ferrie and D. G. Cory, Hamiltonian learning and certification using quantum resources, Phys. Rev. Lett., 2014, 112, 190501.
78. J. Wang, S. Paesani, R. Santagati, S. Knauer, A. A. Gentile, N. Wiebe, M. Petruzzella, J. L. O'Brien, J. G. Rarity, A. Laing and M. G. Thompson, Experimental quantum Hamiltonian learning, Nat. Phys., 2017, 13, 551–555.
79. J. R. McClean, N. C. Rubin, K. J. Sung, I. D. Kivlichan, X. Bonet-Monroig, Y. Cao, C. Dai, E. S. Fried, C. Gidney, B. Gimby, P. Gokhale, T. Häner, T. Hardikar, V. Havlíček, O. Higgott, C. Huang, J. Izaac, Z. Jiang, X. Liu, S. McArdle, M. Neeley, T. O'Brien, B. O'Gorman, I. Ozfidan, M. D. Radin, J. Romero, N. P. D. Sawaya, B. Senjean, K. Setia, S. Sim, D. S. Steiger, M. Steudtner, Q. Sun, W. Sun, D. Wang, F. Zhang and R. Babbush, OpenFermion: the electronic structure package for quantum computers, Quantum Sci. Technol., 2020, 5, 034014.
80. Cirq, https://github.com/quantumlib/Cirq.
81. A. R. Striganov and N. S. Sventitskii, Tables of Spectral Lines of Neutral and Ionized Atoms.Springer, New York, 1968.
82. K. Huber, Molecular Spectra and Molecular Structure IV. Constants of Diatomic Molecules. Springer, New York, 1979.
83. J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush and H. Neven, Barren plateaus in quantum neural network training landscapes, Nat. Commun., 2018, 9, 4812.
84. E. Farhi, J. Goldstone, S. Gutmann and M. Sipser, Quantum computation by adiabatic evolution, arXiv, 2000, preprint, arXiv:quant-ph/0001106, https://arxiv.org/abs/quant-ph/0001106.

---

Footnote

† Electronic supplementary information (ESI) available: Definitions of quantum gates, construction of the quantum circuit of the time evolution operator, computational conditions for numerical simulations of the BxB algorithm, convergence behaviour of the BxB algorithm. See DOI: 10.1039/d0sc04847j

---