Quantum computing quantum Monte Carlo with hybrid tensor network for electronic structure calculations

Quantum computing quantum Monte Carlo (QC-QMC) is an algorithm that can be combined with quantum algorithms such as variational quantum eigensolver (VQE) to obtain the ground state with higher accuracy than either VQE or QMC alone. Here we propose an algorithm combining QC-QMC with hybrid tensor network (HTN) to extend the applicability of QC-QMC for the system beyond the size of a single quantum device, called HTN+QMC. For HTN in a two-layer quantum-quantum tree tensor, HTN+QMC for an O (cid:0) n 2 (cid:1) -qubit trial wave function can be executed by using only a n -qubit device excluding ancilla qubits. Our algorithm is evaluated on the Heisenberg chain model, the graphite-based Hubbard model, the hydrogen plane model, and MonoArylBiImidazole using full configuration interaction QMC. We found that the algorithm can achieve energy accuracy several orders of magnitude higher than VQE or QMC, and HTN+QMC gives the same energy accuracy as QC-QMC when the system is appropriately decomposed. Moreover, we develop a pseudo-Hadamard test technique that enables efficient overlap calculations between a trial wave function and a standard basis state. In a real device experiment using the technique, we obtained almost the same accuracy as the statevector simulator, indicating the noise robustness of HTN+QMC.


I. INTRODUCTION
Computationally accurate prediction of physical properties enables us to speed up the development of functional materials such as batteries [1,2], catalysis [3], and photochemical materials [4,5].The physical properties are mainly governed by electrons in materials, and the computational cost of calculating electronic structures increases exponentially with system size in general, which often prevents classical computers from achieving the required accuracy to predict the properties.
Quantum computers are expected to solve such classically intractable problems in quantum chemistry and materials science [6].Current quantum computers are called noisy intermediate scale-quantum (NISQ) devices [7], which have limitations on the numbers of both qubits and quantum gates due to physical noise, and various approaches for the NISQ devices are proposed [8,9].One of the most popular algorithms is the variational quantum eigensolver (VQE) [10] which is used to obtain the ground state energy by minimizing energy cost function using a variational quantum circuit, where its parameters are updated by a classical computer.In contrast to the traditional quantum algorithms such as the quantum phase estimation [11], VQE requires much smaller hardware resources.However, VQE still suffers from many issues such as insufficient accuracy [12] and vanishing parameter gradients so-called barren plateau [13].* shu.kanno@quantum.keio.ac.jp In addition to the studies to avoid those issues [14][15][16][17], there are multiple variants of quantum Monte Carlo (QMC), which were proposed for further relaxing the hardware requirement [18][19][20][21][22][23][24].QMC is a computational method that uses stochastic sampling techniques to solve large quantum many-body problems such as molecular systems containing hundreds of electrons [25,26].To date, various types of QMC approaches have been proposed including variational Monte Carlo [27,28], auxiliary field quantum Monte Carlo [29,30], and full configuration interaction quantum Monte Carlo (FCIQMC) [31], where we adopt FCIQMC in this study.In FCIQMC, which is useful for quantum chemistry, a stochastic imaginary-time evolution is executed in the space of all the standard basis states (such as the Slater determinants) that can be constructed from a given spatial orbital basis.As in Refs.[18,21,22], a type of QMC called quantum computing QMC (QC-QMC) was introduced.Combined with a quantum algorithm such as VQE, QC-QMC is used to improve the accuracy of energy evaluation in ground-state calculations by mitigating the sign problem, which otherwise causes an exponential statistical error in the energy evaluation [32].Specifically, the wave function distribution, i.e., the walker distribution, is classically generated, and energy is evaluated by using a trial wave function prepared by VQE.In contrast to VQE where the accuracy depends on parameter optimization, QC-QMC has no such optimization while the accuracy relies on the quality of the trial wave function in the energy evaluation.The ground state computation without optimization is expected to lower the hardware requirements, and in fact, there was a demonstration on a 16 qubit diamond, the largest hardware experiment run nearly within the chemical accuracy [18].
While QC-QMC is expected to provide highly accurate ground-state calculations, applying this method to large-scale systems is one of the remaining challenges.Although the electronic correlations outside the active space can be obtained by classical post-processing [18], the size of the available active space is limited by that of the trial wave function.The proposals for constructing wave functions larger than the size of a quantum device include the divide-and-conquer algorithm [33,34], embedding theory [35][36][37], circuit cutting [38,39], perturbation [40], and tensor network [41,42].Hybrid tensor network (HTN) [42] is a general framework of the tensor network that can be implemented on a reduced number of qubits and gates by decomposing the wave function in the original system into smaller-sized tensors, and these tensors are processed by either quantum or classical computation, i.e., quantum/classical hybrid computation.Such a decomposition can reduce both the effective width and depth of the circuit, which is more robust to noise than the original circuit.There are several quantum algorithms or ansatze depending on the tensor network structure, which include the matrix product states [43], projected entangled pair states [44], and tree tensor networks (TN) [41,45].In the two-layer TN, which is adopted in this study, the quantum states of the subsystem in the first layer are integrated with the second layer either by a quantum or classical method.We refer to the former and latter types as quantum-quantum TN (QQTN) and quantum-classical TN, where the deep VQE [34] and entanglement forging [46,47] can be broadly classified into the quantum-quantum and quantum-classical ones, respectively.If gaining the quantum advantage (rather than noise robustness) is a priority, QQTN is preferred.
In this study, we propose an algorithm of QC-QMC in combination with HTN of a two-layered QQTN.In particular, we consider the following QQTN: ( The tensor network representation of this state is depicted in Fig. 1(a).The tensors ψ ⃗ i ( ⃗ θ U ) (φ im ( ⃗ θ Lm )) are defined using wave functions |ψ( ⃗ θ U )⟩ ( |φ im ( ⃗ θ Lm )⟩) as ψ ⃗ i ( ⃗ θ U ) = ⟨ ⃗ i|ψ( ⃗ θ U )⟩ ( ⟨ ⃗ j m |φ im ( ⃗ θ Lm )⟩), where ⃗ i = i 1 i 2 . . .i k ( ⃗ j m = j m1 j m2 . . .j mn ) being a k (n)-qubit binary string.ψ ⃗ i ( ⃗ θ U ) and φ im ( ⃗ θ Lm ) are assumed to be constructed by quantum circuits, where ⃗ θ U and ⃗ θ Lm are variational parameters of the upper tensor (blue), and the m-th lower tensor (orange), respectively.The detail  2), where we omit the parameters.See Sec.IV C for details. (2) HT N ) , (2) where the observable O is defined using tensor prod- HT N , ⃗ θ HT N ) is calculated by first performing a quantum computation on the lower tensor, classically processing the obtained results to generate the contracted operators N m , and then contracting them again with a quantum computation on the upper tensor.Here, the specific expressions of Eq. ( 2) when performing the algorithm are and where |ϕ h ⟩ is the h-th standard basis state as in a Slater determinant, and O is implicitly included in Eq. ( 4) as O = I ⊗nk .The HTN+QMC algorithm consists of two steps; we first optimize ⃗ θ HT N = ( ⃗ θ U , ⃗ θ Lm ) to prepare the trial wave function, by calculating Eq. (3).In the second step, we execute QMC by using the trial wave function ψ HT N ( ⃗ θ HT N ) , which was optimized in the first step, and Eq. ( 4).We will denote HTN+VQE when referring to the first step only and omit the parameter in |ψ HT N ⟩ hereafter.In this formalism, thanks to the contraction technique shown in Fig. 1(b), when k ∼ n, we can construct O n 2 -qubit trial wave function in QMC by using only n-qubit device excluding ancilla qubits.
We benchmarked the performance of HTN+QMC for the Heisenberg chain model, the hydrogen (H 4 ) plane model, the graphite-based Hubbard model, and MonoArylBiImidazole (MABI), which are classified into four material categories: physical or chemical (solid or molecule) and basic or applied as in Fig. 2. The first two models are commonly used as the benchmark [18,49], Graphite is a two-dimensional layered material and is used as an anode material for lithium-ion batteries [50], and MABI is a model system of the photochromic radical dimer PentaArylBiImidazole (PABI) [51].We also execute HTN+QMC on a real device ibmq kolkata by the newly developed algorithm to calculate the overlap between the trial wave function and the standard basis state required in HTN+QMC.We call this algorithm a pseudo-Hadamard test because the Hadamard test type circuit is used in the algorithm.We mention that the proposed algorithm can be applied to other types of HTN structures and QMC, and the overlap technique can be extended to calculate a transition amplitude between two wave functions.

A. Benchmarking models
The performance of HTN+QMC is benchmarked with the Heisenberg chain model, the graphite-based Hubbard model, the hydrogen plane model, and MABI.The Heisenberg chain model, as shown in Fig. 2(a), is defined as a chain of k clusters consisting of four sites where The intra-and inter-cluster interactions are 1 and J inter , respectively.In the benchmark, we consider k = 2 and 3, i.e., 8-and 12-qubit models and J inter = 0.2, 0.4, . . ., 2.0 Hartree. Figure 2(b) shows our graphite-based Hubbard model (the graphite model hereafter), where two layers of graphene sheets are modeled with periodically aligned unit cells, in which two carbon atoms reside in each layer.Using two qubits to represent the up and down spin orbitals (p z orbitals) in each carbon atom, we define the Hamiltonian with 8 qubits as where q is the spin-orbital index for the p z orbital in the carbon, and q = 1, 2, 3 and 4 (5, 6, 7 and 8) corresponds to the first (second) layer.a † q (a q ) is the creation (annihilation) operators on the q-th site, and n q is the number operator n q = a † q a q .t 1 and t 2 are the hopping energy between the first and second nearest neighbor sites corresponding to the intra-and inter-layer interaction energy, respectively, and U is the on-site Coulomb energy.The prefactors for t 1 and t 2 are due to periodic boundary conditions, e.g., the prefactor for t 2 is 2 because two inter-layer interactions exist per carbon (one inside the unit cell and the other outside the unit cell).The reason for t 2 being only on two indices, q = 1 and 2, is that graphite is AB stacking.We determine the value of t 1 , t 2 , and U by the electronic structure calculation.For details see Appendix D.
The Hamiltonians for the hydrogen plane model (8qubit model) in Fig. 2(c) and MABI (12-qubit model) in Fig. 2(d) were constructed using their restricted Hartree-Fock orbitals of their ground states as molecular orbital bases.In the hydrogen plane model, the two hydrogen molecules are arranged vertically as in the index (ID) 0 of Fig. 2(c).Then the intermolecular and intramolecular distances are shortened and lengthened, respectively, until the shape of the hydrogen plane becomes a square as in ID 3. The two structures are interpolated between ID 0 and 3, thus benchmarking a total of four structures.
Next, we present the specific conditions for constructing the benchmarking models.The decomposition settings used for individual models are as follows; for the Heisenberg chain model, the cluster and even-odd settings refer to decomposing the model into subsystems by the cluster and even and odd indices, respectively, as illustrated in Fig. 3(a); for the graphite model, the horizontal and vertical settings refer to decomposing the model horizontally and vertically with respect to the sheet (in the ab-axis plane in Fig. 2(b)), respectively, as in Fig. 3(b); for the chemical models (hydrogen plane model and MABI), the HOMO-LUMO, alpha-beta, and occ-unocc settings refer to decomposing the model using the pair of HOMO−x and LUMO+x as units (x = 0 and 1 in the hydrogen plane model and x = 0, 1, and 2 in MABI), based on alpha and beta spin-orbitals, and based on occupied and unoccupied orbitals, respectively, as in Fig. 3(c), where the orbitals of the chemical models are shown in Appendix D; the no-decomposition setting refers to the calculation without decomposing the original system, i.e., not using HTN.The cluster setting for the Heisenberg model, the horizontal setting for the graphite model, and the HOMO-LUMO setting for chemical models are adopted as a default.
Figure 4 shows the quantum circuit called real amplitude ansatz, a popular hardware efficient ansatz employed for devices with linear connectivity layout.Figure (a) shows an original circuit, which is used for the statevector, whereas, in order to run on real devices, a variant with ancilla qubit was considered as in figure (b), where details on the operation will be given in Sec.II B 3. In both ansatze, the rotation angle of the RY gate (i.e., the parameter) is updated at the VQE and HTN+VQE iterations while in the second step of QC-QMC and HTN+QMC, it is fixed.Because HTN+VQE demands quite a high computational cost due to its iterative process, the real device (ibmq kolkata) was used only for the overlap calculation in the second step, i.e., not for HTN+VQE computations.
The depth represents the value for the real amplitude ansatze, which are k-qubit and n-qubit circuits for upper and lower tensors, respectively, in HTN, whereas nkqubit circuits in the no-decomposition setting (in which ancilla qubit is not counted).Here, k and nk represent the number of subsystems and the number of system In all the QMC methods (QMC, QC-QMC, and HTN+QMC), we evaluate the mean and standard deviation of the (mixed) energy over 5,000 to 10,000 iterations.Hereafter (absolute) deviation from the exact ground state is referred to by energy difference.See Appendix E for the details of the conditions of VQE and QMC.

B. Numerical results
First, we demonstrate the performance improvements in HTN+QMC for HTN+VQE for the Heisenberg chain model, which is then followed by the analysis of all the benchmarking models.Then for the Heisenberg chain model, we discuss the HTN+QMC performance focusing on the dependencies on the decomposition settings.Finally, we describe the operational principle of the pseudo-Hadamard test technique and present the results of the real device experiments for the hydrogen plane model and MABI using that technique.In the HTN+QMC result, the energy difference from the exact value is 5.3 × 10 −3 Hartree, which is an improvement by two orders of magnitude from the result for HTN+VQE.In addition, compared to the QMC result, the HTN+QMC results show higher accuracy and FIG. 4. Circuits of real amplitude ansatz in this study.dH , dN , and dH are the depth in HTN, the no-decomposition setting, and HTN (real device), i.e., the block surrounded by the dotted line is repeated dH or dN , and dH times, respectively.(a) The circuit in statevector procedure.All the lines represent system qubits.(b) The circuit in real device procedure.The topmost line represents an ancilla qubit and the other lines represent system qubits.less energy fluctuation.Specifically, the energy differences (standard deviations) for QMC and HTN+QMC are 5.5 × 10 −1 and 5.3 × 10 −3 Hartree (7.7 × 10 −1 and 4.8 × 10 −2 Hartree), respectively.Therefore, the energy accuracy would be improved in HTN+QMC compared to either HTN+VQE or QMC alone, and we will discuss the details of this improvement in the next subsection.

Performance improvements
Here we briefly comment on the other models, for which the computational results are described in detail in Appendix F. Table I shows the results summarized on HTN of the Heisenberg chain models with k = 2 and 3, the graphite model, the hydrogen plane model with ID 3, and MABI.In all the models, the difference for HTN+QMC is one to several orders of magnitude smaller than that for either HTN+VQE or QMC alone.Especially for the hydrogen plane model, the difference for HTN+QMC is one order of magnitude smaller than that for QMC, even though the fidelity for HTN+VQE is 0.44, almost the same as that for the leading single standard basis state in the exact ground state (single reference state hereafter), 0.45.The above result indicates that the performance of QMC may be improved even if the trial wave function does not have high fidelity, and it is important to evaluate the quality of the trial wave function through the QMC computation.

Dependences on the decomposition settings
We consider the Heisenberg chain model with k = 2, n = 4, and J inter = 1.0, which was computed in the cluster, even-odd, and no-decomposition settings.Prior to the results, we introduce an average interaction strength between the subsystems G mr as a simple measure of decomposition defined as where δ amr = 1 if Pauli X, Y , or Z operator is included in r P amr of the Hamiltonian, which is defined in Sec.IV D, in more than two subsystems, and δ amr = 0 otherwise.G mr in the cluster setting and even-odd setting are 1.5 and 10.5 Hartree, respectively, i.e., the cluster setting seems to be suitable for the trial wave function preparation.In Appendix G, we also show the values of G mr for the graphite model, the hydrogen plane model with ID 3, and MABI, and we have confirmed that the results, especially in the physical models, are almost consistent with intuition.For example, in the graphite model, G mr of the horizontal setting has one order of magnitude smaller than that of the vertical setting.The chemical models especially in the hydrogen plane model, however, give results that differ from intuition, which indicates that more sophisticated measures might be needed in complicated systems.The mutual information, which represents the degree of dependence between random variables, could be used as the alternative measure [52,53].We also show the values using the exact ground state in the appendix, which tends to be more consistent with intuition than G mr .Although the exact ground states cannot be obtained in realistic situations, in a previous study [52], the classical construction of approximate value using the density matrix renormalization group (DMRG) is proposed.Figure 6(a), (b), and (c) show the decomposition dependencies of the fidelity, the energy difference, and the standard deviation for the Heisenberg chain model.In cases of the decomposition settings, the fidelity (Fig. 6(a)) in the cluster setting 0.92 is four times higher than in the even-odd setting 0.22, as expected from the values of G mr .We also calculated the bipartite entanglement entropy of the exact ground state when decompos-ing the cluster and even-odd, which is 0.66 and 3.46.It is evident that an appropriate choice of the decomposition setting is crucial in the performance of the trial wave function.
In addition, the fidelity for the cluster setting is equal to or slightly higher than that for the no-decomposition setting with d N = 10 whereas the fidelity for the evenodd setting is lower than for the no-decomposition setting with d N = 2. Here, the numbers of the parameters for the cluster setting with d H = 4 and the no-decomposition setting with d N = 10 are nk(d H +1)+k(d H +1) = 50 and nk(d N +1) = 88, respectively.That is, the cluster setting with fewer parameters shows comparable performance to the no-decomposition setting.The increase in the fidelity is obviously related to the decrease in the energy difference in HTN+VQE and HTN+QMC as in Fig. 6(b), and that in the standard deviation in HTN+QMC as in Fig. 6(c).A similar tendency for the fidelity, energy difference, and standard deviation can be found in the graphite model, the hydrogen plane model, and MABI, see Appendix G for the details.From the results, we can expect that if the system is appropriately decomposed, i.e., if the interaction between the subsystems is small, HTN can prepare a trial wave function that performs as well as or better than the wave function generated by the quantum circuit of the original system size.Note that the fidelities for the hydrogen plane model with ID 3 are almost the same or lower than that of the single reference state depending on the decomposition settings as shown in Appendix G.One of the ways to increase the fidelity is an improvement of the initial wave function used in the VQE or HTN+VQE.As shown in the appendix, the fidelity can increase in some cases by using the initial state which is close to the Hartree-Fock state.In addition, while the parameters of all the tensors are sequentially optimized in this study, separate optimization of one tensor by another could be an alternative [54].
The high performance in the cluster setting seems to arise from that the wave function represented by the QQTN in Fig. 1(a) has a problem-inspired structure while the ansatz actually used for each tensor is hardware efficient one.In order to analyze the wave function in detail, we first compare the distribution of the absolute coefficient of the standard basis state, which was calculated by using one of the 10 random seeds in Fig. 6 (a), (b), and (c).As shown in Fig. 6(d), the cluster setting exhibits a distribution very close to that of the exact ground state, in comparison to the even-odd setting.Note that for the even-odd setting, the distribution was calculated with the basis state encoding, where the qubit index was reordered from 8, 7, 6, 5/4, 3, 2, 1 to 8, 6, 4, 2/7, 5, 3, 1 as shown in Even-Odd of Fig. 3(a).This qubit index was then changed back to that of Cluster in the plotting, allowing for direct comparison with the even-odd setting.
We then examine the wave function, which is represented as a linear combination of tensor products of the subsystems.In the case of the cluster setting for the Heisenberg chain model, we can approximately describe the exact ground state |ψ g−C ⟩ by four dominant terms with respective coefficients of c 1 = −0.37,c 2 = 0.24, c 3 = −0.23,and c 4 = 0.18, as The state is energetically favored because of the (sub)antiferromagnetic spin configurations, where there is at most one spin pair with the same parity (e.g., 00 and 11), considering that all terms in the Hamiltonian of Eqs. ( 5), (6), and ( 7) are positive.Each term consists of a primitive basis state (in bold) and the spin and spatial inverted ones, where the ground state redefined in this primitive basis {|0101⟩ , |1010⟩} ⊗ {|0101⟩ , |0110⟩} is regarded as a separable state (with the same number of up/down spins in each subsystem), which can be accurately prepared by the lower tensors with one leg connecting to the upper tensor.In case of the even-odd setting, where the qubit index is reordered and the number of up (or down) spins in the subsystem takes 0 to 4, the state in such a primitive basis, e.g., , is no longer separable, resulting in a lower fidelity than in the cluster setting.Through the above mechanism, the decomposition setting affects the fidelity, where matching between the tensor network state and the target state is a crucial factor, being in analogy with the discussion in the area low with tensor networks.Note that the number of legs limits the dimension of the basis for state preparation.
The HTN+VQE results with the cluster setting in Fig. 6(d) show that the coefficients for the above 12 basis states are of the same magnitude as c 1 , c 2 , . . ., c 4 .In contrast, the even-odd setting gives a much less accurate wave function; half of the basis states in Eq. (10) was negligibly small in magnitude, and the coefficient of |0101⟩ ⊗ |0101⟩ is 0.78, which is so large leading to quite different distribution from the exact.
Next, Fig. 6(e) shows the distribution of the nodecomposition setting, where the distributions of the ground state and cluster settings are reproduced from the figure (a) for comparison.The depth is set to d N = 6 such that the numbers of parameters in the nodecomposition and cluster settings become comparable, that is, 50 and 48, respectively.The distribution for the no-decomposition setting deviates from that for the exact wave function, and half of the basis states in Eq. (10) were negligibly small in magnitude.In general, the nodecomposition setting would achieve higher performance than HTN.However, in the cases where the number of ansatz parameters is restricted, HTN may perform better if we find an ansatz and decomposition suited to the structure of the system.In fact, the coefficient of the basis state |00101101⟩, which has the 10-th largest magnitude in the no-decomposition setting, is 0.14, whereas it is only -0.057 in the ground state and is -0.00018 in the cluster setting.Thus, the cluster setting can efficiently prepare a trial wave function that incorporates the correlations of the system with fewer parameters by eliminating the basis states that have a small contribution to the ground state of the system.Note that these observations are verified by the calculation over 10 random seeds.

Real device experiments
We consider performing the proposed HTN+QMC algorithm on a real device.For this purpose, we need to reduce the computational cost on the overlap between the trial wave function and the standard basis state (as in Eq. ( 4)).Conventionally, the overlap ⟨ψ|ϕ h ⟩ between a wave function |ψ⟩ = U ψ |0⟩ ⊗ns and the standard basis state |ϕ h ⟩ = U ϕ h |0⟩ ⊗ns , can be calculated using the Hadamard test as in Fig. 7(a), where U ψ typically consists of a collection of one-and two-qubit gates depending on the ansatz and n s is the number of system qubits.The controlled-U ψ contains a much larger number of multiqubit gates compared to U ψ because all one-and twoqubit gates in U ψ are modified to the controlled gates.Several techniques to calculate the overlaps in fermionic systems have been proposed, that can avoid deep circuits by using the particle preserving ansatz [21,55].However, the techniques are inadequate for HTN+QMC because even for fermionic systems, the number of electrons is not preserved in the subsystem, and thus particle nonpreserving ansatze (e.g., hardware efficient ansatz) must be used.Note that U ϕ h requires at most n s two-qubit gates because it applies the CNOT gates from the ancilla to target qubits, the latter of which must be set to |1⟩ in |ϕ h ⟩.
We develop a gate-efficient technique to calculate the overlap executable on an arbitrary ansatz, called the pseudo-Hadamard test.See Appendix H for the extension of the technique to HTN.
In fact, the gate Ũψ can be determined through the VQE, which is formulated as where the index of the ancilla qubit is set to 0. Fig- shows the quantum circuit used to run such VQE, where O iq (i q = 0, 1, . . ., n s ) is an observable of the i q -th qubit.The constraints can be relaxed through the appropriate selection of ansatz, such as the real amplitude ansatz without rotation gates on the ancilla qubit, as shown in Fig. 4(b).In this case, the VQE can be simplified as We adopt this simplified formalism hereafter.Note that Ũψ |0⟩ |0⟩ ⊗ns takes either |0⟩ |0⟩ ⊗ns or − |0⟩ |0⟩ ⊗ns in the constraint, but the signs cancel out in a numerator and denominator of E mix in both QC-QMC and HTN+QMC, allowing one to ignore the sign.In general, the sign might be able to be determined based on the overlap between a trial wave function and a specific standard basis state, e.g., the Hartree-Fock state [21].We also mention the extensions of the technique; a transition amplitude of ψ O ϕ h can be calculated by the measurements of each system qubit in the basis corresponding to O iq for the pseudo-Hadamard test circuit; We applied the pseudo-Hadamard test technique to the HTN+QMC calculations on a real device ibmq kolkata for the hydrogen plane model with ID 3 and MABI.The best trial wave function was selected from 100 seeds of HTN+VQE runs under the constraint ( dH = 4 and statevector execution), where at most five qubits including an ancilla qubit are used in the overlap calculation.See Appendix H 1 for details on the extension to HTN. Figure 8(a) shows the result for the hydrogen plane model.A smaller fluctuation around the exact ground state energy was found in HTN+QMC than in QMC.The standard deviations obtained by the QMC, HTN+QMC (statevector) and HTN+QMC (ibmq kolkata) methods are 6.7×10 −3 , 1.5×10 −3 , and 2.0×10 −3 Hartree, respectively.Importantly, there is only a slight difference in the fluctuation between the real device and statevector, indicating that there would be some robustness of HTN+QMC to noises.The previous study [18] showed that QC-QMC is robust to the depolarization noise, and we found by the analysis shown in Appendix H 2 that such robustness also holds in the current settings of HTN+QMC.The MABI results in Fig. 8(b) also showed high accuracy in the real device experiment, where the standard deviations are 2.6 × 10 −3 , 1.6 × 10 −3 , and 1.8 × 10 −3 Hartree in the above three methods, respectively.These numerical and theoretical results suggest that HTN+QMC can be a promising candidate for calculating large systems beyond the scale of quantum devices while maintaining accuracy.Note that although we utilized the readout error mitigation and the dynamic decoupling in the overlap calculation, we confirmed that a comparable level of accuracy was achieved even without those options in HTN+QMC for the hydrogen plane model, where the standard deviation without the options is 1.5 × 10 −3 Hartee.

III. DISCUSSION
We proposed an algorithm HTN+QMC that combines QC-QMC with HTN for accurately calculating problems in quantum chemistry beyond the size of a quantum device: (1) QC-QMC can perform electronic structure calculation with higher accuracy than either VQE or QMC alone by using a trial wave function prepared on quantum devices, (2) HTN enables us to construct a large wave function beyond the size of a quantum device by decomposing the wave function into smaller-size tensors, which are computed by hybrid quantum/classical methods.Specifically, the two-layered QQTN, which is one of the HTN structures, allows QC-QMC to prepare O n 2qubit trial wave function by using only O(n) qubit devices.
The trial wave function prepared by HTN+VQE is used to execute HTN+QMC, where FCIQMC is adopted as a QMC method.As demonstrated on the benchmarking models, we found that HTN+QMC exhibits energy accuracy that is several orders of magnitude higher than HTN+VQE or QMC.In addition, HTN+QMC will outperform QC-QMC (i.e., without HTN), if the interaction between the decomposed subsystems becomes small, that is, if the target state is close to a separable or low entangled product state in a suitably chosen decomposition setting.We developed the gate-efficient pseudo-Hadamard test technique involving an ancilla qubit, for the purpose of conducting the HTN+QMC experiments on real devices; the at most five qubits experiment showed an accuracy comparable to the statevector simulation for the hydrogen plane models (8 qubit model) and MABI (12 qubit model).
While this study assumed that the size of the target system is larger than that of the quantum device, there may be cases in which the proposed algorithm should be used even if the size of the target system is the same as that of the quantum device.For example, a quantum computer with thousands of qubits will appear in the near future [58], but due to noise in the quantum device, it is possible that an accurate solution may not be obtained in the calculation of a chemical model of that size.In such a case, decomposing the system into subsystems of about hundreds of qubits would lead to a more accurate result than that obtained by directly executing on a thousand qubits.
The research on QC-QMC has not yet been extensive, and the application of techniques developed for NISQ devices to QC-QMC is a possible future research item, as in HTN in this study.Another interesting direction is the application of QC-QMC to the fields outside the electronic structure calculations such as machine learning for designing advanced materials.

A. Variational quantum eigensolver
VQE takes a variational approach to obtain the ground state of a given Hamiltonian, where a trial wave function ψ( ⃗ θ) is generated by a quantum circuit with variational parameters ⃗ θ, which is repeatedly updated using the classical computer until a termination condition is satisfied, e.g., the expectation value of the Hamiltonian ψ( ⃗ θ) H ψ( ⃗ θ) takes a minimum.The expressibility of ψ( ⃗ θ) depends on the quantum circuit, what is called ansatz.There are several problem-inspired ansatze such as the unitary coupled cluster ansatz [59] and Hamiltonian variational ansatz [60], and despite the high performance, the ansatze require deeper quantum gates.In real device experiments, compromised strategy with hardware-efficient ansatze [61] is often taken.
The Hamiltonian H in the electronic structure problems can be represented as where c a is the a-th coefficient of H, P ab ∈ {X, Y, Z, I} is the a-th Pauli or identity operator on the b-th site.In this study, c a ∈ R is assumed for all the models (c a ∈ C in general).The number of terms is estimated to be O N 4 so with N so being the number of spin-orbitals [62].We choose the Jordan-Wigner mapping [63] for the fermionqubit translation.

B. Quantum computing quantum Monte Carlo
In QC-QMC Ref. [18], a stochastic algorithm such as the imaginary-time evolution is used to iteratively update the discretized coefficients of the wave function i.e., walkers.There are two approaches to integrate the quantum computations: one is to use it only for the energy evaluation [21], and the other for the energy evaluation and walker control [18,21,22].The latter approach, which is more effective in mitigating the sign problem, is computationally expensive [21].We adopt the former approach which is easier to execute on a real device.The following mixed energy E mix is used as a common energy estimator: where |ξ⟩ is a trial wave function.|ψ QM C ⟩ denotes a wave function generated by QMC, which is defined by where w h is the h-th coefficient and is expressed using discretized units (walkers where H |ψ g ⟩ = E g |ψ g ⟩ is used.When the sign problem occurs, the statistical error for evaluating E mix becomes exponentially large [32].For example in FCIQMC, for basis states with small w h , the sign of w h fluctuates over the iterations by competing positive and negative walker counts [64,65].However, by using a trial wave function |ξ⟩ that incorporates the electron correlations correctly, we can mitigate the sign problem.A trivial case is that, if the trial wave function coincides with the ground state, i.e., |ξ⟩ = |ψ g ⟩, then we can obtain E g for any |ψ QM C ⟩ with zero variance [66].See Appendix B for a detail of the derivation.In contrast, the trial wave functions that are conveniently available in classical computations, such as the Hartree-Fock state, the linear combination of mean-field states, and the Jastrow-type states [25] may easily bring about the sign problem.Therefore, if we can only prepare a suboptimal trial wave function in the exponentially large Hilbert space, we can considerably minimize the errors in the energy estimation. In QC-QMC, the trial wave function |ξ⟩ is prepared by a quantum algorithm.Specifically, by decomposing |ψ QM C ⟩, Eq. ( 15) is rewritten as The matrix elements ⟨ϕ h | H |ϕ h ′ ⟩ can be obtained by trivial classical calculations.In contrast, quantum computation is required to calculate the overlap ⟨ξ|ϕ h ⟩, for which we can apply the Hadamard test or the classical shadow [67,68].We discuss in Appendix B that the variance of the mixed energy will decrease with the fidelity of the trial wave function in a simplified case.In this study, fidelity is defined as the square of the overlap between the target and ground states assuming the absence of noises.

C. Hybrid tensor network
In the two-layer QQTN, the original nk-qubit system is decomposed into k subsystems of n qubit states φ im , which are combined together and integrated with the k-qubit state |ψ⟩.These states can be associated with tensors, which are defined as φ im ⃗ jm = ⃗ j m φ im and ψ ⃗ i = ⃗ i ψ with the vector indices ⃗ j m = j m1 j m2 . . .j mn and ⃗ i = i 1 i 2 . . .i k in n-qubit and k-qubit binary strings, respectively.We call the former a lower tensor and the latter an upper tensor, respectively.The wave function |ψ HT N ⟩ is then decomposed into tensor products of φ im which is the same as Eq. ( 1) although omitting the parameters, where the number of vector coefficients {ψ ⃗ i } is 2 Lk with L denoting the number of bits (legs) used to construct each i m .L = 1 is adopted in this study.If L = n, |ψ HT N ⟩ becomes the general nk-qubits wave function.
When L ≪ n, |ψ HT N ⟩ lives in a subspace much smaller than the entire 2 nk dimensional Hilbert space, although it can be larger than the subspace consisting only of the classical tensor due to the exponentially large ranks of ψ ⃗ i and φ im ⃗ jm , and the performance of the QQTN depends on the decomposition setting of the system, which determines the entanglement between the subsystems.
There is a variety of tensor implementations in a quantum circuit.For example, we can define a lower tensor as ⊗n , etc., where the first two formulations give the index in initial states and the third a unitary matrix by using unitary matrices U Lm and U im Lm .In the present study, we choose ⊗n−1 and |ψ⟩ = U U |0⟩ ⊗k for lower and upper tensors, respectively.
In our QQTN formulation, an observable is defined as a tensor product O = m,r O mr , where O mr is observable on the r-th qubit of the m-th subsystem (lower tensor).Transition amplitude is then defined as which is the same as Eq. ( 2) although omitting the parameters, where ψ are the two different of states for |ψ HT N ⟩.As will be explained in the next section, T is used both in VQE and QMC to evaluate an observable and overlap, respectively.We first calculate r O mr φ im (2)  for the lower tensors by quantum computations, construct N m = N 00 N 01 N 10 N 11 , which is a trivial classical process, and then integrate the results as T = ψ (1)  m N m ψ (2) on the upper tensor by quantum computations.Appendix C shows the details of this procedure, which is based on the circuit of the Hadamard test as in Ref. [69].Our algorithm requires only max(n, k) qubits except for ancilla qubits.In this study, only 8k + 2 terms are measured in the T calculation, i.e., the overhead for calculating the expectation value is a linear scale in the system size nk.
Note that the number of measurements is 2 × 4 L k + 2 in general.

D. Proposed algorithm: HTN+QMC
The procedure of HTN+QMC consists of two steps.2. Perform QC-QMC by using the obtained trial wave function |ψ HT N ⟩; the quantum computer is used to compute E mix to accurately estimate the ground state energy.
Here, we will denote HTN+VQE when only the first step is mentioned.HTN+QMC can be performed for nkqubit system by using only O(max(n, k)) qubits except for ancilla qubits.Specifically, if k ∼ n, an O n 2 -qubit trial wave function is prepared only by an O(n)-qubit circuit.In both steps, the Hamiltonian and mixed energy can be evaluated through the calculations of T in Eq. ( 20).
In the first step, by splitting the index b in Eq. ( 14) into the two indices m and r, we can rewrite H and its expectation value as In the second step, the overlap ⟨ξ|ϕ h ⟩ = ⟨ψ HT N |ϕ h ⟩ in Eq. ( 18) is calculated by substituting ψ (1) 20), where the circuit parameters are fixed to the values obtained in the first step.More specifically, we can prepare an arbitrary basis state ⊗nk by setting U (2) Lm = m X jmr(h) , where j mr (h) is a function of h, j mr (h) takes a value on {0, 1}, and X is a Pauli X operator.Having the overlaps ⟨ψ HT N |ϕ h ⟩ of all |ϕ h ⟩ corresponding to the walkers appearing during the QMC execution at hand, we can perform QMC through iterative evaluations of the mixed energy (Eqs.( 15) and ( 16)).In FCIQMC, the wave function |ψ QM C ⟩ is defined as where w h is the real coefficient of h-th standard basis state (such as the Slater determinant) |ϕ h ⟩.Those coefficients will be updated iteratively through the imaginary time evolution according to which can be written in a discrete form as is the Kronecker delta, τ is an imaginary time, ∆τ is a time increment, and S is an energy shift.The first and second terms are contributions from the off-diagonal and diagonal elements of the Hamiltonian, respectively.
The procedure of FCIQMC is illustrated in Fig. 9.The walker with the positive (negative) sign is depicted as red (blue) rectangular.The imaginary time in the i Qth iteration is denoted by τ i Q .For the walkers in the i Q -th iteration, called the parent walkers, the following operations are performed to obtain the i Q + 1-th walkers.
• Spawning step [the off-diagonal terms in Eq. (A3)]: For each parent walker with index h, a new walker with index h ′ is generated with a probability of • Death/cloning step [the diagonal term in Eq. (A3)]: For each parent walker with an index h, remove (copy) walkers of the same index with probability copy) the walkers with probability 1.
• Annihilation step: Merge the walkers obtained in the spawning and death/cloning steps.The walkers with opposite signs are canceled out.
In the death/cloning step, the initial value of S is a constant (e.g., the value of the Hartree-Fock energy), and when the number of walkers exceeds the set value N shif t , S is updated every A iteration as S( ) so that the number of the walkers remains constant (called variable shift mode), where S(τ i Q ) is S on the time τ i Q , ζ is a damping parameter, and N W (τ i Q ) is the number of walkers on τ i Q .We note that although the computational cost, i.e., the number of walkers required by this algorithm, may increase exponentially with the size of the system in general, there are several proposals to make the algorithms executable for practical systems [1][2][3][4].For example, the number of walkers required for N 2 molecule can be reduced by three orders of magnitude (10 8 to 10 5 ) by restricting spawning only on the basis state for which the number of walkers exceeds a certain threshold [2].
Finally, the other calculation conditions that are not described in Sec.II A are the following.The initial value of the energy shift S was set to the energy of the leading single reference state in the standard basis state of the ground state (the ground state when constrained to six electrons in MABI).The imaginary time increment ∆τ was 0.001 Hartree −1 for the Heisenberg chain and graphite models, and 0.1 Hartree −1 for the hydrogen plane model and MABI.The parameters for the variable shift were N shif t = 1000, A = 5, and ζ = 0.1.

Appendix B: Statistical analysis of mixed energy in simple cases
Here, in the simple cases, we discuss the statistical error of the mixed energy, which arises from preparing the trial wave function by a quantum algorithm.Let ϵ and (1 − F ) denote the deviations of the wave function generated by QMC |ψ QM C ⟩ and the trial wave function |ξ⟩ from the ground state |ψ g ⟩, the wave functions can be represented as and respectively.F 2 is a fidelity of |ξ(F )⟩, ψ ⊥ s and ψ ⊥ s ′ are the orthogonal states with the coefficients |ψ g ⟩, u s and v s ′ , and ϵ, F ∈ [0, 1].The mixed energy E mix (ϵ; F ) can be then described as with the ground state energy E g and eigenvalue of ψ ⊥ s denoted by E s .Below, we assume E g = 0 without loss of generality.Then the mean and variance of E mix (ϵ; F ), E(E mix (ϵ; F )) and V ar(E mix (ϵ; F )), are respectively represented as and where the integral variable of the expected value is the QMC time step, on which ϵ and u s depend.Now consider E(E mix (ϵ; F )) and V ar(E mix (ϵ; F )) in the two cases.In the first case of F = 1, both E(E mix (ϵ; F )) and V ar(E mix (ϵ; F )) are trivially zero, that is, if the trial wave function is the ground state, E mix is equal to the ground state energy regardless of |ψ QM C (ϵ)⟩ (excluding ϵ = 1).
In the second case, where E mix (ϵ; F ) fluctuates around the ground state energy in a small amount, assuming E(E mix (ϵ; F )) = 0 and the expected value of the higher order terms of ϵ be negligible, the variance becomes 2 is monotonically decreasing, and thus when the fidelity of the trial wave function in the quantum algorithm is larger than that in the classical algorithm, we can expect a decrease in the variance of the mixed energy.The calculation for the upper tensor.
Appendix C: Calculation of the transition amplitude by using HTN We explain the procedure for calculating T in Eq. ( 20) based on the Hadamard test [5], which is reproduced as The first part includes the calculations of the matrices N i ′ m im for the lower tensors.Figure 10(a) shows the quantum circuit for this purpose when the state in the lower tensor is defined by U init is one of the circuits in the lower panel, which is selected for each of the matrix elements N 00 , N 01 , N 10 , or N 11 , and creates the superposed state Lm and U Lm operations, and measured in the measurement basis corresponding to O mr .We can obtain the real and imaginary parts of )) by specifying respective X and Y for the measurement basis on the ancilla qubit.Calculation of N i ′ m im for each lower tensor, therefore, requires 4 × 2 = 8 different measurements, where factor 4 comes from N 00 , N 01 , N 10 , and N 11 , and the factor 2 from Re(N i ′ m im ) and Im(N i ′ m im ).The second part is the calculation of the upper tensor.In terms of the upper tensor state, T can be rewritten as where the 2×2 non-Hermitian matrix N m is reformulated by the singular value decomposition (SVD), where and V m are unitary matrices, and D m is a diagonal matrix.Since N m is a 2 × 2 matrix, the SVD can be executed classically.The quantum circuit for calculating Eq. (C3) is shown in Figure 10(b), where the superposed state gives the diagonal element of D m , i.e., (D m ) 00 or (D m ) 11 , respectively [5].The real and imaginary parts of T can be obtained by specifying X and Y for the measurement basis on the ancilla qubit, respectively.The calculation for the upper tensor thus requires measurements on two different bases.In total, the calculations of the lower and upper tensor require 8k +2 sets of measurements, i.e., the overhead for calculating the expectation value is a linear scale for the system size nk.Note that when the Hamiltonian and wave function are real as in this study, the calculations of the real parts in the lower and upper tensors are unnecessary, reducing the overhead to 4k+1.We additionally note that the QQTN can be expressed as a nk-qubit quantum circuit ( m U (l) ⊗n ) under the current assumption of φ im(l) and ψ (l) , where Lm operates qubits corresponding to m-th subsystem (i.e., m-th |0⟩ ⊗n ), and U U operates the qubits corresponding the first qubit of each subsystem [6].As a final remark, the robustness of the barren plateau in this type of circuit has been reported [7].
The MABI structure was determined by S 0 state geometry optimization at the three-state-averaged CASSCF(14e,14o)/6-31G level using Molpro 2015 [23,24], where MABI geometry in Cartesian coordinates is shown in Table II.The PySCF package [25] was adopted for constructing the Hamiltonian for the hydrogen plane models and MABI with STO-3G and 6-31G being as the respective basis functions.The orbitals of the hydrogen plane model with ID 2 and MABI are shown in Fig. 11.The CASCI(6e,6o) problem was considered for preparing the 12-qubit Hamiltonian in MABI.HTN+QMC codes were implemented by the Python, especially the NumPy [27] and SciPy [28] packages, and Qiskit [29] and Qulacs [30] were adopted as quantum circuit simulators.The sequential least squares programming optimizer (SLSQP) in the SciPy package was used in the parameter optimizations.In MABI, although there are six electrons in the ground state in the CASCI(6e,6o) problem, the number of electrons in the ground state of the corresponding qubit Hamiltonian is not six.Thus VQE/AC [31] algorithm was used to perform energy optimization with the constraints imposed on the number of electrons.The ground state of MABI refers to the eigenstate with the smallest eigenvalue acquired by restricting the electron number to six.
In all the QMC methods (QMC, QC-QMC, and HTN+QMC), the initial number of walkers was one, and the single reference state was selected as the initial configuration.In the hydrogen plane model (except for ID 3) and MABI, we checked out that the Hartree-Fock state was chosen as the single reference state.The trial wave function was the single reference state, and the states obtained by VQE and HTN+VQE were used for QC-QMC and HTN+QMC, respectively.The maximum number of iterations in all the types of QMC was 10,000.

Appendix F: Benchmark results
We benchmark the models with circuit depth d H = 1, 2, . . ., 6.As an example, we first explain the benchmark results of the Heisenberg chain model in the cluster setting with k = 2, especially J inter = 0.2, 1.0, and 2.0 Hartree. Figure 12(a), (b), and (c) show the energy difference, standard deviation, and fidelity versus circuit d H , respectively.For HTN+VQE in Fig. 12(a), on the one hand, the energy difference is more than 10 −1 Hartree when J inter is large, e.g., for J inter = 1.0 Hartree, the difference is 4.3 × 10 −1 Hartree even with d H = 6.On the other hand, all the values of the difference for HTN+QMC are less than 10 −1 Hartree with d H ≥ 2. In addition, the difference for HTN+QMC with d H ≥ 3 is more than one order of magnitude less than that of QMC.In Fig. 12(b), the decrease in the standard deviation from QMC is also confirmed in HTN+QMC.From the results, we find that HTN+QMC can evaluate the energy more accurately than the previous classical and quantum algorithms.worsens as the fidelity decreases, for example of d H = 4, the fidelity with J inter = 0.2, 1.0, and 2.0 Hartree is 1.00, 0.92, and 0.49, respectively, and the energy difference with J inter = 2.0 for HTN+QMC in Fig. 12(a) is one order of magnitude larger than that of J inter = 0.2 and 1.0.Nevertheless, the energy difference with d H = 4 and J inter = 2.0 for HTN+QMC is one order of magnitude less than that for QMC, i.e., the accuracy of HTN+QMC is still higher than that of QMC even when the fidelity is low.Therefore, even if the trial wave function prepared in HTN+VQE is not highly accurate, HTN+QMC can be expected to provide better energy accuracy than the classical calculation.Note that we also calculate the bipartite entropy for the models in Fig. 12(d), which tends to increase with J inter .
Finally, we show the benchmark results for the Heisenberg chain model with k = 3, the graphite model, the hydrogen plane model, and MABI as in Figs. 13, 14, 15, and 16, respectively.The results of bipartite entanglement entropy are shown for the models with k = 2. From these data, we determined that d H = 4 would provide sufficient accuracy and adopted this value as the default.
Appendix G: Remaining results of the decomposition settings We discuss the results of the Heisenberg chain model with k = 2 and J inter = 1.0, the graphite model, the hydrogen plane model with ID 3, and MABI.Table III show the values of G mr defined in Eq. ( 9) and the mutual information.The mutual information for two and three subsystems I(A; B) and I(A; B; C) are defined as H(B, C, A) − H(B, C), and H(A) is von Neumann entropy calculated from the reduced density matrix for the subsystem A. We transformed the wave function, which was originally constructed in an NumPy array to the reduced density matrix using PennyLane [32].Figures 17,18,and 19 show the results of the analysis for the graphite model, the hydrogen plane model with ID 3, and MABI, respectively.We also executed the VQE and HTN+VQE using the initial state close to the Hartree-Fock state; the unitary matrix U corresponding to the real amplitude ansatz with all parameters set to zero satisfies U |0⟩ ⊗nk = |0⟩ ⊗nk , and thus an arbitrary basis state (including Hartree-Fock state) |ϕ h ⟩ can be constructed in the state preparation for VQE by applying Pauli X gates on the appropriate qubits after U .In the following, we show that the same state construction is also possible in HTN+VQE.In case that the basis state corresponds to a computational basis state |ϕ h ⟩ = m,r |j mr (h)⟩ = ( m,r X jmr(h) ) |0⟩ nk , the state can be constucted by substituting U Lm for r X jmr(h) U Lm in HTN, where j mr (h) ∈ {0, 1}: assuming |ψ HT N ⟩ in Eq. ( 19) becomes We performed the VQE and HTN+VQE from the initial states close to the Hartree-Fock state by using the initial parameters which are chosen randomly from [0, 0.01), and the results for the hydrogen plane model with ID 3 and MABI are shown in Figs.18 and 19, respectively.where H3) ŨLm and ŨU are unitary gates that are specified in detail in the next paragraph.
The expectation value of the observable can be rewritten in terms of the upper and lower tensors (both involv-ing the ancilla) as and from the corresponding measurement results denoted by M 0 , M 1 , M + , and M y+ , respectively, we obtain the matrix elements Ñ 00 = M 0 , Ñ ⊗k .Alternatively, Ñm is Hermitian and can be measured using the eigenvalue decomposition of Ñm classically.Now we move on to the overlap T h between the quantum state and standard basis state, which is represented similarly to the observable as The circuit for calculating Ñ i hm in the lower tensor is shown in Fig. 21(a).Since there is no tensor index in |ϕ hm ⟩, only two different circuits are needed for Ũinit (the lower panel of the figure).U ϕ hm is composed of CNOT gates from ancilla to target qubits, which should be set to |1⟩ in |ϕ hm ⟩.Next, we explain the calculation of the upper tensor.Ñ i hm is a tensor with one leg, i.e., can be represented using a two-element vector as Ñ i hm =    qubits are measured only X basis in this study since the trial wave function is real.Finally, the rest of the calculation conditions, which are not described in Sec.II A, are listed below.We checked on the four-qubit Heisenberg chain model to confirm that the amplitudes close to those of the exact ground state were obtained by using the pseudo-Hadamard test when a highly accurate wave function is obtained in HTN+VQE.We used at most five qubits in real device execution, where the initial layout in the device was [12,15,18,21,23], and the number of shots was 4000.We combined the constraints in Eq. (H2) as m ⟨n Lm ⟩ + ⟨n U ⟩ = 0 in the implementation since the left side of each constraint in Eq. (H2) is non-negative for any ŨLm and ŨU because they consists of the expectation values for the particle number.The overlaps are calculated only for the Slater determinants of four and six electrons in the hydrogen plane model and MABI, respectively, with the rest of the amplitudes set to zeros.In the calculation of MABI, in addition to the constraints on the overlap calculation, an additional constraint enforcing the number of electrons to be six was imposed.

Analysis of noise resilience
We discuss the reason for the robustness of HTN+QMC to the noise.In the circuit of Fig.

(H10)
Since T h can be represented by substituting Ũψ = ( m ŨLm ) ŨU , U ϕ h = m U ϕ hm , and n s = nk in Fig. 7(b), we can assume the above discussion holds true, where ŨLm operates the ancilla qubit and qubits corresponding to m-th subsystem, and ŨU operates on the ancilla qubit and qubits corresponding the first qubit of each subsystem.Then the ratio of T h for different standard basis states is Since E mix is an estimator that includes the normalization of the wave function in HTN, HTN+QMC is unaffected by depolarizing noise.Figure 22(a) and (b) show the results of T h in the hydrogen plane model without and with the normalization, respectively.The fidelity is restored from 0.09 to 0.49 through normalization, whereas the fidelity is 0.46 in the statevector.As a result, the noise robustness of HTN+QMC, which was experimentally observed, was supported theoretically and numerically.

FIG. 1 .
FIG. 1. Representation and calculation of hybrid tensor network in the two-layer QQTN.(a) The hybrid tensor network representation of the two-layered QQTN state.(b) The calculation of T in Eq. (2), where we omit the parameters.See Sec.IV C for details.

FIG. 2 .
FIG. 2. Models for benchmarking HTN+QMC.The structure of graphite is drawn by VESTA [48].The qubit indices are labeled in (a) and (b), see Sec.II A for details.(a) The Heisenberg chain model.(b) The graphite-based Hubbard model.(c) The hydrogen plane model.(d) MonoArylBilmidazole (MABI).
ucts O = m,r O mr , and O mr is an observable attached on the r-th qubit of the m-th subsystem (r = 1, 2, . . ., n).As shown in Fig. 1(b), the transition amplitude T ( ⃗ θ

FIG. 3 .
FIG. 3. Decomposition settings.n = 4 and k = 2 are assumed for all the figures.(a) The settings for the Heisenberg chain model.(b) The settings for the graphite-based Hubbard model.(c) The settings for chemical models in the case of the eightqubit model.

Figure 5 (
Figure 5(a) shows the result of executing HTN+VQE for the Heisenberg chain model in the cluster setting with d H = 4, k = 2, and J inter = 1.0 Hartree.As in the inset, the energy difference from the exact ground state energy is 4.4 × 10 −1 Hartree at the end of the optimization.Figure 5(b) shows the result for HTN+QMC.In the HTN+QMC result, the energy difference from the exact value is 5.3 × 10 −3 Hartree, which is an improvement by two orders of magnitude from the result for HTN+VQE.In addition, compared to the QMC result, the HTN+QMC results show higher accuracy and

FIG. 5 .
FIG. 5. Results of the energy in the Heisenberg chain models of the cluster setting with n = 4, k = 2, Jinter = 1.0, and dH = 4.The black dashed line represents the energy of the exact ground state.(a) The result for HTN+VQE (blue line).The inset is an enlarged view along the y-axis.(b) The results for QMC (blue line) and HTN+QMC (orange line).

FIG. 6 .
FIG. 6. Results of the analysis for the Heisenberg chain model with Jinter = 1.0.(a), (b), and (c) are the decomposition dependencies.The average and error bar is obtained over 10 different random seeds used for the initial parameters in VQE or HTN+VQE.The fidelity of the single reference state is plotted by the black bar in (a), and the energy difference and standard deviation by the black cross markers in (b) and (c), respectively.The green (blue) circle and cross markers denote the results for the decomposition (no-decomposition) setting in HTN+VQE and HTN+QMC (VQE and QC-QMC), respectively, with dH = 4 (dN = 2, 4, . . ., 10).In HTN+VQE and HTN+QMC, n = 4 and k = 2 for all the settings.(d) and (e) are the distributions of the wave function for the Heisenberg chain model with Jinter = 1.0.The basis index is the value when the standard basis state is expressed in the decimal number, e.g., |10000010⟩ corresponds to 130.The black, red, green, and blue bars represent the coefficients for the exact ground state, HTN+VQE result with the cluster setting, HTN+VQE result with the even-odd setting, and VQE result with the no-decomposition setting, respectively.(a) The fidelity.(b) The energy difference.(c) The standard deviation.(d) The comparison of the exact ground state, cluster setting, and even-odd setting.(e) The comparison of the exact ground state, cluster setting, and no-decomposition setting (dN = 6).
FIG. 7. Quantum circuits for calculating the overlap.The topmost line represents the ancilla qubit and the other lines represent the system qubits.(a) The Hadamard test.(b) The pseudo-Hadamard test.(c) The VQE circuit with constraint to determine Ũψ used for the pseudo-Hadamard test.The ancilla qubit is set to |0⟩ when calculating the constraints.

FIG. 8 .
FIG. 8. Results of the energy on the real device execution for the hydrogen plane model and MABI.HOMO-LUMO setting and dH = 4 are adopted for the models.The blue line represents the QMC result, the light green and red lines respectively represent the HTN+QMC results of statevector and real device procedures, and the black dashed line represents the exact ground state energy.(a) The result of the hydrogen plane model with n = 4 and k = 2. (b) The result of MABI with n = 4 and k = 3.The inset in each figure is an enlarged view along the y-axis.

=
and ⟨ψ HT N | H |ψ HT N ⟩ = a c a ⟨ψ HT N | m,r P amr |ψ HT N ⟩ , (22) respectively.The expectation value is evaluated through Eq. (20) by setting ψ |ψ HT N ⟩ and replacing O mr with P amr , where the coefficient index a is regarded as implicitly included in the expression of O mr .The wave function |ψ HT N ⟩ is prepared by parameterized quantum circuits equivalent to the unitary matrices U Lm and U U .
Appendix A: Procedure of full configuration interaction quantum Monte Carlo

FIG. 9 .
FIG.9.Procedure of FCIQMC.The walker of |ϕ h ⟩ with a plus (minus) sign is represented as a red (blue) rectangle.

FIG. 10 .
FIG. 10.Quantum circuits for calculating T .The topmost line represents an ancilla qubit and the other lines represent system qubits.(a) The calculation for the lower tensor.(b) The calculation for the upper tensor.
and k m=1 V m , and measured in the computational basis.The measurement result, |0⟩ or |1⟩

FIG. 11 .
FIG. 11.Orbitals for the hydrogen plane model and MABI, where the structures are drawn by Jmol [26].(a) The orbitals of the hydrogen plane model with ID 2. (b) The orbitals of MABI.

Figure 12 (
FIG. 12. Results of the Heisenberg chain model of the cluster setting with n = 4, k = 2, and dH = 1, 2, . . ., 6.The values of Jinter are shown in color bars.The circle and cross marks denote the results of HTN+VQE and HTN+QMC, respectively.The dashed lines in (a) and (b) [(c)] denote the results of QMC [a single reference state].(a) The energy difference.(b) The standard deviation.(c) The fidelity.(d) The bipartite entanglement entropy.
FIG. 13. Results of the Heisenberg chain model of the cluster setting with n = 4, k = 3, and dH = 1, 2, . . ., 6.The values of Jinter are shown in color bars.The circle and cross marks denote the results of HTN+VQE and HTN+QMC, respectively.The dashed lines in (a) and (b) [(c)] denote the results of QMC [a single reference state].(a) The energy difference.(b) The standard deviation.(c) The fidelity.

FIG. 15 .
FIG. 15. Results of the hydrogen plane model of the HOMO-LUMO setting with n = 4, k = 2, and dH = 1, 2, . . ., 6.The values of ID are shown in color bars.The circle and cross marks denote the results of HTN+VQE and HTN+QMC, respectively.The dashed lines in (a) and (b) [(c)] denote the results of QMC [a single reference state].(a) The energy difference.(b) The standard deviation.(c) The fidelity.(d) The bipartite entanglement entropy.

FIG. 17 .
FIG. 17. Results of the analysis for the graphite model.The average values and error bars were calculated over 10 different random seeds used for the initial parameters in VQE or HTN+VQE execution.For reference, values obtained from the single reference state and those from QMC are shown by the black bar in (a) and cross in (b) and (c), respectively.The green (blue) circle and cross denote the results from the decomposition (nodecomposition) setting in HTN+VQE and HTN+QMC (VQE and QC-QMC), respectively.In HTN+VQE and HTN+QMC, n = 4 and k = 2 for all the settings.(a) The fidelity.(b) The energy difference.(c) The standard deviation.

FIG. 19 .
FIG. 19.Results of the analysis for MABI.The average values and error bars were calculated over 10 different random seeds used for the initial parameters in VQE or HTN+VQE execution.For reference, values obtained from the single reference state and those from QMC are shown by the black bar in (a) and cross in (b) and (c), respectively.The green circle, green cross, blue circle, and blue cross respectively denote the results from the decomposition, decomposition, nodecomposition, and no-decomposition settings in HTN+VQE, HTN+QMC, VQE, and QC-QMC, and the yellow circle, yellow cross, purple circle, and purple cross respectively denote those when the initial state is close to the Hartree-Fock state in HTN+VQE, HTN+VQE, VQE, and VQE.In HTN+VQE and HTN+QMC, n = 4, k = 3 for the HOMO-LUMO setting, and n = 6, k = 2 for the other two settings.(a) The fidelity.(b) The energy difference.(c) The standard deviation.

FIG. 20 .
FIG. 20.Quantum circuits for calculating HTN+VQE to execute the pseudo-Hadamard test.The topmost line represents an ancilla qubit and the other lines represent system qubits.The circuit with the ancilla qubit set to |0⟩ is used when calculating constraints.(a) The calculation for the lower tensor.(b) The calculation for the upper tensor.

iÑhm
FIG. 21.Quantum circuits for calculating T h by using pseudo-Hadamard test.The topmost line represents an ancilla qubit and the other lines represent system qubits.(a) The calculation for the lower tensor.(b) The calculation for the upper tensor.

FIG. 22 .
FIG. 22. Results of T h in the hydrogen plane model without and with normalization in the real device procedure.The light green and red bars represent the values in the statevector and real device procedure, respectively.(a) The result without normalization.(b) The result with normalization.
). |ϕ h ⟩ is the h-th standard basis state as in a Slater determinant and can be represented by a binary string i.e., a computational basis.The procedure of generating the wave function |ψ QM C ⟩ depends on the method of QMC.In this work, we chose FCIQMC.See Appendix A for the details of its procedure.E mix is not exactly equivalent to the expectation value of the Hamiltonian ⟨ψ QM C |H|ψ QM C ⟩ / ⟨ψ QM C |ψ QM C ⟩ (that is, not the pure estimator), but we can obtain the exact ground state energy E g when |ψ QM C ⟩ approaches to the ground state |ψ g ⟩, i.e., |ψ QM C ⟩ ∼ |ψ g ⟩ as

TABLE II .
MABI geometry.The unit is Angstrom.