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

The Heisenberg exchange coupling parameter J (H = −2JSi · Sj) 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 jS2, the wave function overlap estimation with the SWAP test, and Bayesian optimisation of the parameter j. Numerical quantum circuit simulations for H2 under a covalent bond dissociation, C, O, Si, NH, OH+, CH2, NF, O2, and triple bond dissociated N2 molecule revealed that the BxB can compute the J value within 1 kcal mol−1 of errors with less computational costs than conventional quantum phase estimation-based approaches.


Introduction
Quantum computing and quantum information processing (QC/QIP) is one of the most innovative research elds in modern science. A quantum computer uses a quantum bit (qubit) as the minimal unit of information, which can possess not only either |0i or |1i but also arbitrary superposition of the |0i and |1i states. 1 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". 2 Among the diverse topics in the eld 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 conguration 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. 3 A quantum-classical hybrid algorithm known as a variational quantum eigensolver (VQE) executable on the noisy intermediate-scale quantum (NISQ) devices 4 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][23][24][25][26][27][28][29][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 lightemitting materials, 34,35 and so on. Open shell multi-nuclear transition metal complexes are oen involved in enzymes as reactive centres. 36,37 Single molecule magnets have been extensively studied as molecular memory devices. 38 To disclose their electronic structures, sophisticated ab initio quantum chemical calculations are powerful and essential tools. However, in open shell systems carrying spin-b unpaired electrons, the wave function has a strong multi-congurational 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][40][41] to determine spin quantum numbers of arbitrary wave functions, 42 and to purify the spin contaminated wave functions, 43 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][45][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 denition 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 lowspin 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 "Bayesian exchange coupling parameter calculator with broken-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 À1 or even smaller and therefore chemical precision is necessary to evaluate the value quantitatively, and (2) determining the energy down to the ne 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 , (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 À1 of errors regardless of the magnitude and sign of the exchange coupling parameter J.
This paper is organised as follows. In Section 2, we briey 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 2 molecule under a covalent bond dissociation, singlet-triplet energy gaps of C, O, Si, NH, OH + , CH 2 , NF and O 2 are given. N 2 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.

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 DE 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 DE S-T , 54 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, a, b, and 0 specify doubly occupied, singly occupied by a spin-a electron, singly occupied by a spin-b electron and unoccupied orbitals, respectively. The multi-congurational nature of the open shell singlet wave function originates from the symmetry requirement of the spin operator S 2 . 55,56 The wave function is a simultaneous eigenfunction of the electronic Hamiltonian H and the S 2 operator, but single Slater determinant carrying spin-b unpaired electrons is not an eigenfunction of the S 2 operator. The conguration interaction (CI) methods or the multi-congurational methodologies such as complete active space self-consistent eld (CASSCF) are required to treat open shell singlet wave functions explicitly. In fact, the multireference (MR)-CI, MR-CC, and multi-reference perturbation theory (MRPT) with the CASSCF reference wave function are capable of predicting the DE S-T and the J value in a quantitative manner. [57][58][59] In this framework, the energies of high-spin and low-spin states are calculated separately and subtract to obtain the DE S-T and J value. The fact that the open shell singlet wave function has a multi-congurational character prevents us from directly treating open shell singlet states by means of single conguration 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 |J BS i. 60,61 The |J BS i 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, hS 2 i is an expectation value of the S 2 operator. For twospin systems HS corresponds to the spin-triplet state. The E HS can be easily calculated at the single conguration theory because the spin-triplet wave function with M s ¼ 1 is generally well approximated by the single Slater determinant |22/ 2aa00/0i. 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, DE S-T can be calculated from the spin-ip variants of the time-dependent DFT, CI, and equation-of-motion coupled cluster (EOM-CC) methods. [62][63][64][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 ip excitations from the high-spin reference state. In these approaches, DE S-T is calculated as an excitation energy from the high-spin state.
As described here, the calculations of DE S-T and the J value generally require two separate calculations, i.e. the energy calculations of the high-spin and low-spin states.

Quantum chemical calculations on quantum computers
Here we briey 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 denitions of the quantum gates and quantum circuits are given in the ESI. † In QPE the time evolution of a wave function |Fi is simulated conditional on the qubits for readout and the relative phase shi caused by the time evolution is extracted by means of inverse quantum Fourier transformation (QFT À1 ) and following qubit measurements. 66 The quantum state |Fi 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 |hF|J i i| 2 . Here, |J i i is the full-CI wave function of the ith electronic state. By preparing the |Fi 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 spinb unpaired electrons, the Hartree-Fock wave function |J HF i is generally a good approximation of the ground state. In open shell molecules carrying spin-b unpaired electrons, one can use a conguration state function |J CSF i that is an eigenfunction of the S 2 operator as |Fi. The |J CSF i 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 oen adopted to construct the quantum circuit corresponding to the time evolution operator U ¼ exp(ÀiHt). The time evolution operator under the rst and second order Trotter decompositions are given in eqn (11) and (12), respectively.
The quantum circuit for the time evolution operator exp(ÀiwX 1 Z 2 Z 3 X 4 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 delity 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 DE ¼ p/2 kÀ1 t for U ¼ exp(ÀiHt). If we set t ¼ 1 we need k ¼ 12 to achieve $1 kcal mol À1 of energy tolerance. In this case, we have to apply the time evolution operator up to U 2 11 ¼ exp(À2 11 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 q and feedbacks the parameter to the quantum processing unit.
In eqn (13), the |J HF i is generally used as |F 0 i. U(q) 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(q). In VQE, the energy expectation value E(q) is calculated through repeated measurements and therefore it includes statistical errors. It is known that approximately O(w max 2 /3 2 ) of measurements are required to determine the energy within a precision 3, where w max ¼ max(|w m |) in eqn (9). 26 Millions of measurements are required to calculate the energy within 1 kcal mol À1 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 À1 or even smaller. To discuss the spin state energy gap and J value within 1 kcal mol À1 of tolerance, we have to calculate the energy of individual spin states within 0.5 kcal mol À1 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 3 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.

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 specied 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 H ¼ À2JS i $S j 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 DE 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 2 . The eigenvalue of the S 2 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 shied Hamiltonian H 0 with an additional term jS 2 as in eqn (14).
H 0 and H have the simultaneous eigenfunctions, because [H, S 2 ] ¼ 0. The eigenvalue of the shied Hamiltonian H 0 is given in eqn (15).
Eqn (15) . 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 nd the j parameter in the shied Hamiltonian H 0 ¼ H + jS 2 , where the two spin states have the same eigenvalue of H 0 . The energy gap is calculated as DE LS-HS ¼ j[S HS (S HS + 1) À S LS (S LS + 1)]. Fig. 2 illustrates a schematic energy diagram of the shied Hamiltonian H 0 with J > 0. In the case of J < 0, the crossing point of the S ¼ 0 and S ¼ 1 lines appears on the le hand side of the gure. 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.
Theorem 2. If the system can be approximated by the two-site Heisenberg spin Hamiltonian with effective spinsS $ 1/2 as given in eqn (16), all spin states have the same eigenvalue of H 0 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 2 molecule (S ¼ 3/2), which involves six unpaired electrons.
From theorem 1 the high-spin and low-spin states have the same eigenenergy of the shied Hamiltonian H 0 if j ¼ J. Thus, the calculation of the J value can be accomplished by nding the j parameter so as for the broken-symmetry wave function |J BS i in eqn (18) to become an eigenfunction of H 0 .
If |J BS i is an eigenfunction of H 0 , the time evolution of |J BS i under the shied Hamiltonian H 0 merely gives rise to a phase shi without changing the structure of |J BS i. The deviation of |J BS i from the eigenfunction of H 0 can be estimated from the square overlap |hJ BS |U(j, t)|J BS i| 2 , where U(j, t) is a timeevolution operator dened in eqn (19). The square overlap |hJ BS |U(j, t)|J BS i| 2 can be efficiently evaluated on a quantum computer, by utilising a SWAP test. 73 The quantum circuit to estimate the square overlap |hJ BS -|U(j, t)|J BS i| 2 on a quantum computer is depicted in Fig. 3.
The quantum circuit in Fig. 3 consists of two parts: the time evolution of |J BS i under the shied Hamiltonian H 0 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 le conditionally interchanges the quantum states of qubits storing |J BS i and U(j, t)|J BS i if the control qubit is in the |1i state. The quantum state before measurement is given in eqn (20).
Here, DE 0 denotes an energy difference between the low-spin and high-spin states under the shied Hamiltonian H 0 .
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 lowspin 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 eld of quantum computing, Bayesian inference has been used in QPE 74,75 and Hamiltonian learning. [76][77][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 m ¼ 0 and s 2 ¼ 1 as the starting prior distribution, where m and s 2 denote the mean and variance, respectively. (2) The evolution time t of U(j, t) is set to be 1.2/s 2 for twospin systems. For triple bond dissociation of N 2 molecule we used t ¼ 0.4/s 2 . These evolution time lengths were set so that P(0|j; t) becomes a single maximum when j is taken in the range of m À s 2 to m + s 2 . (3) Draws m samples in the range of m À s 2 to m + s 2 with a constant interval, and executes the quantum circuit of Fig. 3 R 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 tted by a normal distribution and the posterior distribution P(j|0; t) was calculated by using eqn (24). (5) If m(posterior) < m(prior) À s 2 /2 or m(posterior) > m(prior) + s 2 /2, return to step (3) with the j value giving the largest P(0) as m(prior). Note that this and next steps are introduced for the purpose of algorithm stability enhance. (6) Set s 2 (posterior) ¼ s 2 (prior)/5 if s 2 (posterior) is smaller than s 2 (prior)/5. (7) Convergence check. If s 2 (posterior) is smaller than the threshold, the algorithm returns m(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 À1 ), tighter threshold should be adopted.

Results and discussion
In order to demonstrate the BxB quantum algorithm we developed a numerical quantum circuit simulation program using Python with OpenFermion 79 and Cirq 80 libraries and performed the numerical simulations for covalent bond dissociation in H 2 , triplet ground states of C, O, Si, NH, OH + , CH 2 , NF, and O 2 , and triple bond dissociation in N 2 . The simulations were performed ve times for each geometry. The computational conditions for the numerical simulations are summarised in the ESI. †

Covalent bond cleavage in H 2 molecule
To demonstrate the BxB quantum algorithm, we performed the numerical quantum circuit simulations for covalent bond dissociation in H 2 molecule in the range of 1.2Å to 3.0Å. Note that |J BS i is converged to the RHF root for R(H-H) ¼ 1.1Å and shorter. We used the STO-3G basis set and |J BS i is mapped onto four qubits. The 2J ¼ DE 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 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 À1 at all atom-atom distances. The DJ value tends to become large for the shorter H-H distance, which presumably originates from the fact that |J BS i 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 atten 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 DJ 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 AE 0.001 kcal mol À1 at R(H-H) ¼ 1.5Å, which corresponds to DJ ¼ À0.003 kcal mol À1 .
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 nal Bayesian iteration is about 250-300 atomic unit, which should be compared with t $ 8000 a.u. to calculate the DE S-T within the same precision (0.5 kcal mol À1 for the energy difference and hence 0.25 kcal mol À1 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.

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 2 , NF, and O 2 . The experimental DE S-T value is available for these systems. 81,82 Note that |J BS i is approximated by the linear combination of the lowest spin-triplet and the second lowest spin-singlet wave functions in CH 2 and O 2 , rather than the lowest spin-singlet wave function. Therefore, the J value obtained from the BxB algorithm does not correspond to the half of DE S-T in CH 2 and O 2 . To calculate the DE S-T that can be directly compared to the experimental value in CH 2 and O 2 molecules, the initial wave function should be replaced from the broken-symmetry wave function to |J 0 i as dened in eqn (25) and (26), respectively. Here, 2, a, and 0 stand for doubly occupied, singly occupied by a spin-a 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 2 all valence orbitals are included in the active space. In NF and O 2 , the valence s/s* and p/p* orbitals are considered to be in the active space. The number of qubits required for the wave function mapping is summarised in Table 1.
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 ve numerical simulations. The standard deviation of the J(BxB) is less than 0.05 kcal mol À1 . Again, the BxB algorithm reproduced the J(CAS-CI) value within 1 kcal mol À1 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 DJ(BxB À full-CI) becomes À0.20 and À0.26 kcal mol À1 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.

Triple bond dissociation in N 2 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 twosite Heisenberg spin Hamiltonian given in eqn (16) with effective spinsS $ 1/2. Here, we demonstrate the J value calculation in the triple bond dissociation in N 2 molecule with the (6e, 6o) active space consisting of valence orbitals participating in the triple bond (s/s* and p/p* orbitals). By using the STO-3G basis set, |J BS i converges to the open shell state with six unpaired electrons, those of two occupy the s-type orbitals and the rest distribute to the p-type orbitals, at the bond distance R(N-N) ¼ 2.1Å and longer. |J BS i is approximated by the linear combinations of four spin states as given in eqn (27).
By assuming the Heisenberg spin Hamiltonian with effective spinsS ¼ 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 s-type orbitals is different from that in the p-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 2 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 À1 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 †).

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 À1 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 À1 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 delity. 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 innitely 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, 83 because the evolution time length t is set to be as P(0) changes substantially in the range of m À s 2 < j < m + s 2 .
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 2 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 |J BS i 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 |J BS i. 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 AE 0.01 kcal mol À1 , which is very close to that obtained by the SWAP test with all the qubits (J ¼ 18.14 AE 0.01 kcal mol À1 ).
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 |J BS i. 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 conicts to declare.