Variational quantum eigensolver with embedded entanglement using a tensor-network ansatz

In this paper, we introduce a tensor network (TN) scheme into the entanglement augmentation process of the synergistic optimization framework by Rudolph et al. [arXiv:2208.13673] to build its process systematically for inhomogeneous systems. Our synergistic approach first embeds the variational optimal solution of the TN state with the entropic area law, which can be perfectly optimized in conventional (classical) computers, in a quantum variational circuit ansatz inspired by the TN state with the entropic volume law. Next, the framework performs a variational quantum eigensolver (VQE) process with embedded states as the initial state. We applied the synergistic to the ground-state analysis of the all-to-all coupled random transverse-field Ising, XYZ, Heisenberg model, employing the binary multiscale entanglement renormalization ansatz (MERA) state and branching MERA states as TN states with entropic area law and volume law, respectively. We then show that the synergistic accelerates VQE calculations in the three models without an initial parameter guess of the branching-MERA-inspired ansatz and can avoid a local solution trapped by a standard VQE with the ansatz in the Ising model. The improvement of optimizers for MERA in all-to-all coupled inhomogeneous systems, enhancement, and potential synergistic applications are also discussed.


I. INTRODUCTION
The ground state of quantum many-body systems is essential to understand quantum phenomena in many-body physics.However, solving quantum many-body problems is generally difficult with conventional (classical) computers because they require the treatment of Hilbert spaces that grow exponentially with an increasing system size of N .
Using quantum computers has been attracting attention to address the problems' limitations.Actually, quantum algorithms such as the quantum phase estimation [1], the qubitization [2], and the quantum singular value transformation [3] with quantum acceleration, unable us to access the ground state for large quantum many-body systems that intractable for classical computers.
However, current quantum computers are noisy intermediate-scale quantum (NISQ) devices [4] that suffer from errors and are limited in size.Quantum algorithms that can be executed using shallow quantum circuits are required to make them more practical.In this context, one of the NISQ-aware methods developed is the Variational Quantum Algorithm (VQA) [5].This method involves constructing variational models through parameterized quantum circuits (also called ansatz) on a quantum computer and optimizing the cost function using classical computers.VQAs, specifically the variational quantum eigensolver (VQE), have been applied to search for the ground state of quantum systems [6].
To obtain meaningful results using VQAs on NISQ, the ansatz utilized to represent the target state with the required precision should have few internal parameters and be as shallow as possible concerning the circuit depth.This condition is crucial for mitigating the noise effects inherent to NISQ devices and preventing the occurrence of the barren plateau (BP) problem [7], which is associated with an overly excessive variational space.Although the BP problem can be circumvented by starting the calculation from appropriate initial parameters [8], it is generally nontrivial to set such parameters.Therefore, a systematic procedure for setting up a high-quality ansatz is desired.
As a strategy to satisfy this requirement, the procedures of putting more classical computer resources into a series of processes of VQE are currently attracting attention [9,10].The VQE with a synergistic framework [10] is a quantum-classical hybrid algorithm that has received much attention to avoid BP problems and works efficiently with less use of NISQ.The synergistic framework has two key concepts: sophisticated variational methods on classical computers ideally provide the initial state for the VQE, and quantum gates added during the VQE calculation and the quantum circuit in which the classically optimized solution is efficiently embedded have different topologies.The second concept is referred to as the entanglement augmentation process.
The use of the tensor network (TN) method [11][12][13][14] for the first concept of the synergy framework is a natural idea because the TN method for quantum states repre-sents the wave function as a contraction of many tensors with small degrees of freedom (TN representation) and is a way to realize the variational method of large quantum many-body systems on a classical computer; there is an ongoing effort to convert TN states into quantum circuit representations today [15][16][17].
Actually, in Ref. [10], the authors used the matrix product state (MPS) [18][19][20], which has a onedimensional (1D) topology, as a variational method on classical computers, and employed all-to-all topologies of fully parameterized SU(4) gates for the augmentation process.However, as mentioned in Ref [10], constructing an all-to-all topology is a nontrivial task and is likely not scalable on NISQ devices.Therefore, to achieve a synergistic framework irrespective of the task that can be implemented on near-to-mid-term quantum computers, it is necessary to propose a systematic entanglement augmentation process that can qualitatively and dramatically change the entanglement structure of the ansatz by only adding very sparsely coupled quantum gates.
In this paper, we introduce the concept of a TN structure into this entanglement augmentation process and propose a procedure to efficiently and systematically augment the entanglement representation capability of a quantum circuit ansatz.Specifically, we focus on embedding a TN state with the entropic area law [21] after variational calculation on a classical computer in a portion of the entanglement-augmented quantum circuit ansatz, which is equivalent to a TN state with the entropic volume law.This entangled embedding VQE (EEVQE) approach can accelerate the calculations of VQE with a TN-inspired ansatz with an entropic-volume law.Discussions of embedding focused on TN with the entropic volume law and its connection to VQE, which is an attempt not seen in pioneering works [17,22,23] related to our study.Investigations of the TN state with the entropic volume law are essential for analyzing quantum states in complex problems, such as quantum chemical systems, long-range interacting systems, random spin systems, and nonlinear problems [24] which involve non-uniform interactions among all qubits, and simulations of real-time evolution for quantum systems because the quantum states may not satisfy the entropic area law in these systems.For specific numerical verification of EEVQE, we employ a 1D binary multi-scale entanglement renormalization ansatz (MERA) [25,26] for variational methods on a classical computer and use branching MERA networks [27,28] as the quantum circuit ansatz after entanglement augmentation.The target Hamiltonians are all-to-all coupled random transversefield Ising/XYZ/Heisenberg models whose ground states can violate the entropic area law.As we expect, through the EEVQE process, we confirm for three models that TN states after the variational optimization on classical computers can be further improved in terms of variational energy while benefiting from the initialization-problemfree and accelerating VQE calculations.Furthermore, we also confirm that EEVQE can avoid traps in locally opti-mal solutions (or broadly BP) in the random transverse field Ising model.
As part of our efforts to advance the synergistic framework, we compared optimizers to update MERA for allto-all coupled random systems.This mainly includes the Evenbly-Vidal method [26], the standard in MERA optimization, its improvements developed in this work (refer to Sec.III.2), and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method, which is the standard in parameterized quantum circuit ansatz optimization.Our comparison results showed that the BFGS method was the most effective in reducing the variational energy for all the studied models.This result indicates the importance of introducing the knowledge of optimizers used in VQE to accelerate the computation of variational methods with TN on the classical computer in the application of the synergistic framework to investigate low-energy states of quantum chemical calculations, which is the most desired application of VQE for all-coupled inhomogeneous systems.
The remainder of this paper is organized as follows: in Sec.II, we review the basics of the components of EEVQE, namely the TN structures of MERA and branching MERA, Evenbly-Vidal sequential optimization method, VQE, and the synergistic framework; in Sec.III, we describe our propositions, i.e., the procedure of embedding MERA states into a quantum circuit ansatz inspired by branching MERA states and the procedure of a modified Evenbly-Vidal method; in Sec.IV, after investigating the performance of optimizers for MERA states in a specific set of all-to-all coupled models, we present the results of a performance study of EEVQE for the models.Finally, in Sec.V, we summarize the study and discuss future studies on EEVQE.

II. REVIEWS FOR COMPONENTS OF EEVQE
In this section, we review the fundamentals of the components of EEVQE: the TN structures of 1D binary MERA and branching MERA, optimization methods for MERA, the calculation procedure of VQE, and the synergistic framework.One familiar with these components may skip to the next section, which discusses how to embed the optimized MERA state in branching-MERAinspired quantum circuits, and the scrutiny of optimization methods for MERA in the all-to-all coupled model.Here, let us consider the 1D binary MERA, whose schematic diagram is depicted in Fig. 1(a).The binary MERA network for finite-size systems consists of fourth-, third-, and second-order tensors, called disentanglers u, isometries w, and top tensors t (see Fig. 1).In the following, to simplify the discussion, the degrees of freedom (bond dimension) of each leg of all tensors are assumed to be the same positive integer χ (see the original paper in MERA [26] to extend the discussion to the case where the bond dimension is tensor-dependent).Each u, w, and t on the network always satisfies the following orthonormality or isometric conditions where δ αα ′ with α ∈ {a, b, c, d} is the Kronecker delta.These conditions are shown schematically in Fig. 1(b).
Owing to these conditions, the quantum state |Ψ⟩ represented by MERA is always normalized, that is, ⟨Ψ|Ψ⟩ = 1.This study considers a position-dependent TN, where each tensor has different tensor elements depending on its placement in the network.
To achieve a more expansive variational space, the 1D branching binary MERA in Fig. 2 introduces a branching structure into isometry while keeping layer structures of the MERA.In the case of the branching binary MERA, the orthonormality of the isometry in which the bifurcation is introduced is identical to that of the disentangler u.Since the branching structure can be arbitrarily introduced in each isometry layer in the binary MERA as shown in Fig. 2, variations of 2 n−1 -pattern branching structures are possible in a 2 n -site system.This paper fo- cuses only on a full-branching MERA, which introduces branching in all isometry layer , for example, branch 2 in Fig. 2 because only the full-branching MERA satisfies the entropic volume law [24].Of course, the more branching we adopt, the larger Hilbert space the branching MERA can search so we may consider branching patterns depending on the target system.

II.2. Optimization scheme for MERA
This subsection presents the Evenbly-Vidal algorithm, which is used as the standard optimization of MERA.Given the Hamiltonian H of the targeted system, the variational energy E can be evaluated as ⟨Ψ|H|Ψ⟩, where |Ψ⟩ is the variational state in terms of the binary MERA.We now consider that the Hamiltonian consists of the sum of the two-body interactions h ij between i-and j-th sites, namely H = i<j h ij .In this case, the expectation value of each interaction e ij = ⟨Ψ|h ij |Ψ⟩can be summed to obtain the variational energy E. For practical calculations, we should perform the minimum contraction, which depends on the (i, j)-pairs by utilizing the isometric condition of u, w, t, or causal cone as shown in Fig 3, and the typical computational cost of evaluating e ij is known to be O(χ 9 log N ) if the system size is sufficiently large [26].For example, given a MERA, suppose we want to optimize an isometry w while keeping the rest of the tensors fixed.The energy E depends quadratically on w and w † , namely where ⟨i < j⟩ refers to the set of site pairs specifying the two-body pair interaction h ij that contributes to the optimization of w, and E ij is the environment tensor obtained by hollowing out the diagrams corresponding to w and w † from the causal cone necessary to evaluate e ij [26], and C = ⟨i<j⟩ e ij with i<j = ⟨i<j⟩ + ⟨i<j⟩ is a constant term with respect to the update of w.However, no known algorithm exists to solve a quadratic problem while maintaining isometric constraints.
Evenbly and Vidal developed a method for optimizing a single tensor of MERA based on a linearizing approximation [26].In this approach, we temporally regard w and w † as independent tensors and optimize w while keeping w † fixed, namely where the tensor Υ w is called the environment tensor of the w, for example, as shown in Fig 4.
To achieve a unique global minimization of Ẽ(w), a singular value decomposition (SVD) is applied to the environment tensor.The tensor w new that minimizes Ẽ(w new ) is always given by as like being refer to Fig. 5.For this update to always work well, the Hamiltonian should be redefined as H γ = H−γI, where γ is sufficiently large so that H γ is negative definite.However, the optimization step size is scaled by γ −1 [29], so we should choose γ as small as possible practically.In the Evenbly-Vidal method in non-uniform MERA, this linearized update is performed sequentially for each tensor that constructs the MERA.Similarly, we can obtain e ij and optimize TN with the branching binary MERA; the typical numerical cost of evaluating e ij is reported to be O(χ 12 N ) if the system size is sufficiently large [27].Compared with the computational cost of MERA and branching MERA, the latter is exceedingly high, especially in higher-dimensional systems [28].This has been the primary factor preventing us from analyzing the performance of branching MERA using classical computers.However, in the case of quantum computers, there is no significant difference in the computational cost with respect to the number of gates.Specifically, the numbers of gates in the MERA and branching MERA are O(N ) and O(N log N ), respectively.

II.3. Variational Quantum Eigensolver
The variational quantum eigensolver (VQE) [6] is a quantum-classical hybrid variational method for finding the ground state energy of a target Hamiltonian H on an N -qubit system.In the VQE, the variational wave function for the system is expressed as a product of parameterized quantum unitary gates Diagram of the environment tensor Υw for tensor w.This is the summation of the tensor contraction, except w from eij, and (i, j) is involved with w in a causal cone.where θ = (θ 1 , θ 2 , • • • , θ K ) and ûk (θ k ) means k-th quantum gate acting on single or multiple qubits specified by the user and has variational parameters is the number of internal degrees of freedom for ûk .In this study, each ûk is an SU(4) gate, where θ k becomes a 15-dimensional real vector, acting on the i k -th and j k -th qubits with 1 ≤ i k < j k ≤ N .The variational wave function |Ψ⟩ is always normalized due to the unitarity of the gates.
A goal of VQE is to minimize the variational energy where with respect to θ.Since numerically exact values of the partial derivatives (2.12) with the orthonormal unit vector e k,l for shifting only the parameter θ k,l are easily obtained by the parameter shift rule [30], we can employ various gradient-based algorithms as optimizers of the VQE.In practical NISQ applications, systematic perturbations linked with gate operations and statistical errors arising from the limited number of measurements influence energy evaluation.However, we could take positive views that numerous methodologies for error mitigation have been developed to obtain the derivatives given by Eq. (2.12).In VQE, the calculation of expectation values {E(θ), E(θ ± π 2 e k,l )} is performed on a quantum computer because the calculation is a bottleneck in simulations using classical computers.Subsequently, the classical computer updates the parameters using the gradient ∇E(θ).Of course, based on the parameter shift rules, we may choose an optimizer that does not use the gradient or the Hessian.Finally, VQE achieves a ground-state search by iterative calculations between the procedures of classical and quantum computers.
The hyperparameters of the VQE are the initial values of the variational parameters, including the design of the quantum gate set, choice of the optimizer, and hyperparameters of the optimizer (see the review [31] for more information).These hyperparameters should be properly controlled according to the quantum computer's ability to perform high-precision calculations in VQE.In particular, when using NISQ devices, it is essential to attempt to decrease the total number of quantum gates within an acceptable range of numerical accuracy to reduce the effect of noise.

II.4. The synergistic framework
The problem with VQAs in NISQ devices is that statistical and systematic errors are associated with the measurement and noise, respectively.A naive approach to solving this problem is to reduce the dependence on NISQ devices in existing quantum-classical hybrid algorithms.For example, Okada et al. proposed an efficient protocol for the topological phase analysis of transverse field clusters and toric code models on quantum circuits [9].In this protocol, all optimization of the quantum circuit, including the calculation of expectation values, is performed on a classical computer using causal cone structures, as discussed in Fig. 3 for each local observable.In addition, the expectation value evaluation of nonlocal observables, which are difficult to handle on a classical computer, is performed using the classically optimized quantum circuit.
However, this protocol has constraints on the target system and ansatz.For example, the measured observables used in the optimization are local, and the ansatz comprises local gates with constant depth.To avoid this difficulty, first, the synergistic framework described in Ref. [10] employs the MPS after the variational calculation with TN methods on the classical computer and performs quantum circuit encoding of the MPS by sequential two-qubit-gate decomposition with 1D topology [32,33] as the initial ansatz of the VQE.Second, the framework augments the quantum circuit encoded from the MPS with an additional parameterized quantum circuit.In this step, the internal parameters of the additional quantum circuit are set to zero, so that they correspond to the identity operators.Finally, the framework performs VQE with the augmented circuit as the initial ansatz.The crucial properties of the framework can be summarized in the following two points: the VQE is trapped in the BP when all parameters of the augmented ansatz are given randomly but can avoid the BP by using the encoded quantum circuit; the all-to-all gate topology, which is different from the topology of the MPS, of the quantum circuit added after the quantum circuit encoding effectively contributes to better VQE performance.These importance values are also reported in Ref. [16,23].
It is worth noting that in Ref. [10] of the synergistic framework, the SU(4) gate product is created through a sequential process during the quantum-circuit encoding of the MPS.This process makes it difficult to fully utilize classical computing power as a divide-and-conquer approach cannot be naively introduced.Additionally, although all-to-all topology gates are most commonly used for entanglement augmentation, it is not a systematic procedure.To achieve a quantum advantage in nonuniform systems, such as quantum chemical systems, by the synergistic framework, TN states can be easily converted to quantum circuits.The entanglement structure of the variational quantum state must be changed by adding a few gates during the entanglement augmentation process.

III. ENTANGLED EMBEDDING VQE
To this end, we propose an entangled embedding VQE (EEVQE) protocol, using binary MERA for classical variational optimization and branching MERA augmented from binary MERA on a quantum circuit.The protocol can be summarized in four steps: 1. Perform the variational optimization of the MERA state with the entropic area law on the classical computer for the target Hamiltonian.
2. Obtain the quantum circuit representation of the optimized MERA state using the technique of quantum circuit encoding.
3. Embed the encoded circuit to an entanglementaugmented quantum circuit inspired by the branching MERA state with entropic volume law.
4. Perform VQE with the initial state given by the augmented circuit.
The following subsection will explain how to incorporate MERA states into a quantum circuit ansatz with the same degrees of freedom as the branching MERA state.Additionally, we will introduce a modified version of the Evenbly-Vidal technique for updating MERA during classical variational optimization.It will be used to compare optimizers for updating MERA in all-to-all coupled random systems, as discussed in Sec.IV.2.

III.1. Quantum circuit representation of MERA and embedding process
The disentangler u, isometry w, and top tensor t with the bond dimension χ = 2 n of MERA equivalent to the SU(2 2n ) rotation operator (2n-qubit gate) have different condition of input qubits |0⟩ ⊗n states as shown in Fig. 6(a).We should take a large n in the case of complex problems and want more accuracy; for example, in [26] they use up to n = 4 roughly to benchmark for the translational symmetry 1D critical system, but the larger n we take, the higher the computational cost for all EEVQE processes.The quantum circuit representations of the binary MERA and branching MERA according to the above rule are depicted in Fig. 6 processes of embedding the isometry and top tensors elements, which form the MERA state after variational optimization on a classical computer, into a unitary matrix of SU(2 2n ) involves using the modified Gram-Schmidt method (See Appendix.A).
In the case of n = 1, by use of Cartan decomposition, any SU(4) operator ûij for i-and j-th qubits can be decomposed into a kind of time evolution operator D ij with XYZ-type interaction and SU(2) rotation operators R i for one qubit [34], that means where i , σ y i , σ z i ) are the Pauli operators for i-th qubit.To execute 2n-qubit gates equivalent to the unitary matrix on a quantum computer, decomposing the 2nqubit gates into a product of two-qubit SU(4) gates is essential because once we have the product, we can break down each SU(4) gate into CNOT gates and single-qubit rotational gates [35], which are commonly used in NISQ today.Computer-assisted search and numerical optimization for the circuit decomposition have been widely studied [36][37][38][39][40][41][42][43][44][45][46][47].We can encode the MERA state in a quantum circuit by using decomposition techniques on each SU(2 2n ) gate individually and simultaneously in classical parallel computing.
To explain the embedding process in detail, we refer to Fig. 6(b), where the MERA state is embedded through quantum gates within the blue gate.The branching MERA is represented by a quantum circuit with an additional orange quantum gate alongside the blue gate.The embedding process is accomplished by setting the internal parameters such that all orange circuits function as identity operators.In the case of the VQE calculation, with the circuit as the initial condition, getting stuck in a local solution, we introduce weak noise to the orange (and blue) quantum gates to avoid such a situation [10,16].
Studying the effectiveness of branching MERA through a quantum computer is beneficial, as it entails lower computational costs than MERA simulations on a classical computer.MERA requires an order of SU(2 2n ) gates of O(N ) while branching MERA requires O(N log N ).These costs refer to the computational requirements for producing a quantum state in a quantum computer by operating quantum gates on individual qubits.However, suppose that the computer can compute mutually commutative quantum gates that are spatially separated.Both computational costs can be reduced to O(log N ), which corresponds to the number of layers of the isometry and the disentangler.Consequently, the relative computational cost of branching MERA is similar to that of MERA.

III.2. Improved Procedures for Sequential Tensor
Optimization in MERA Since this study mainly focuses on applications to allto-all coupled random systems, the classical computational part of the optimization of the inhomogeneous MERA becomes important and nontrivial.We modify the optimization method (Evenbly-Vidal method) employed in the original paper [26] to handle this situation.This procedure updates the tensor to be optimized using the singular vector obtained by singular value decomposition of the environment tensor constructed by all tensor contractions except the tensor to be optimized.However, we must perform a considerable number of sweep optimizations and may be trapped in a local minimum with this method, even for MERA networks with sufficient bond dimensions χ to represent the ground state of the target system, reflecting the effect of linearizing the tensor optimization.Therefore, it is unfavorable for EEVQE.
An approach to reduce the number would be to perform diagonalization of the effective Hamiltonian without linearizing the cost function with respect to only the top tensor of the MERA network.The effective Hamiltonian is defined as follows: Suppose that a quantum circuit C providing the MERA state |Ψ⟩ = C |0⟩ ⊗N is shown in Fig. 6(b) with n = 1, and we consider updating the top tensor u 13 labeled as No. 13.Then, we introduce 13 and Ψ ′ mn = ⟨m, n| u 13 |0⟩ ⊗2 with m, n ∈ {0, 1} and rewrite |Ψ⟩ and the expectation value E = ⟨Ψ| H |Ψ⟩ as follows where ) is the effective Hamiltonian referred to when updating the internal degrees of freedom of u 13 .
Here, we assume that the total Hamiltonian H = i<j h ij is given by the sum of the two-body Hamiltonian h ij between i-and j-th sites.Then, we can evaluate the effective Hamiltonian H ′ required to optimize circuit 13 as the sum of the local two-body effective Hamiltonian ) and its diagrammatic representation is shown in Fig. 7.After the contractions, we perform the diagonalization of the effective Hamiltonian H ′ and the eigenvectors {Φ ′ mn } corresponding to the ground states of H ′ are computed.Finally, the update of circuit 13 is completed by embedding (u 13 ) 0,0 mn := Φ ′ mn (3.6) and performing the unitarization as noted in Sec.III.1.
It should be mentioned here that such embedding of eigenvectors obtained by the diagonalization of the effective Hamiltonian into the tensor to be optimized is difficult to perform for tensors other than the top tensor because the MERA network, unlike the tree tensor network (TTN) [48][49][50], contains a loop structure.On the other hand, diagonalization applies to optimizing parameterized quantum gates that act only on the zero kets of the initial state, regardless of the details of the variational quantum circuit ansatz.

IV.1. Hamiltonians and numerical setup
In order to investigate the properties of our procedure, we consider the following all-to-all coupled random transverse field Ising model XYZ model and Heisenberg model in a random field with 1 = (1, 1, 1), where the coupling constants J ij , J ij = (J x ij , J y ij , J z ij ) and the magnetic fields h i and h i = (h x i , h y i , h z i ) are given by uniform random numbers in the range [-1,1).In the actual calculation, the Hamiltonians are pre-processed to be negative definite by shifting the constant terms γ ij , which is defined as The quantity evaluated hereafter is the relative error from the exact ground state energy and its random-averaged values ∆, where we prepare 10 2 realizations for the random coupling constants and magnetic fields.Then, the initial parameters of the MERA state for each realization are common random values.
For the variational update of the MERA state on the classical computer, we adopt the modified Evenbly-Vedal method as discussed in Sec.III.2 to benefit from the TN method on the classical computer (See Sec.IV.2).The sequential update schedule for the tensors was set from the bottom to the top layer, as shown in Fig. 1(a), adopting a left-to-right order within each layer, and the number of iterations for the sequential update is up to 10 3 .In this study, we use the Julia version of the ITensor library [51] for TN calculations and only consider χ = 2 to focus on the exact circuit encoding of the MERA state.
Then, we use the quantum-simulation software Cirq [52] for the quantum circuit encoding of the SU(4) unitary matrix using Cartan's decomposition.Note that in this study, the χ = 2 tensors are analytically encoded exactly into 2-qubit gates so that we can parallelize the procedure regardless of the target structure.The VQE procedure employs the BFGS method, which is a quasi-Newtonian method.The hyperparameter is only the number of BFGS iterations in this step, and the number is taken up to 10 4 in our study.
In this study, we focus on how well the VQE calculation of branching MERA with the optimized MERA as the initial wave function is performed when the effect of noise is not considered, so we utilize a quantum circuit simulator qulacs [53] and the numerical analysis software library scipy [54] for the VQE calculations.

IV.2. Benchmarks of MERA optimization procedure
Since classical optimization of inhomogeneous MERA is important, as mentioned before, we first compare which optimization method performs better in three optimization methods -Evenbly-Vidal, modified Evenbly-Vidal, and BFGS in the random Ising/XYZ model with all-to-all couplings.The results are shown in Fig. 8, where the colored lines and the light-shaded region represent ∆ and its variance, respectively.
The results of our study demonstrate that the modified Evenbly-Vidal method exhibits a marginally superior convergence speed compared to the original method in both models.In addition, our findings indicate that the modified method enhances the convergence error in the transverse Ising model.Moreover, we compared the modified Evenbly-Vidal method and the BFGS for the variational optimization of the MERA state.The results show that BFGS is more effective in reducing the variational energy during the initial optimization phase than the modified Evenbly-Vidal method for all models analyzed.This agrees with the expectation that the sequential optimization based on the TN method may be relatively easy to trap in local minima.
Also, we briefly confirm that the performance of the modified Evenbly-Vidal method improves as the weight of the degree of freedom of the top tensors in TN states increases in Appendix.B.

IV.3. EEVQE calculations
In this subsection, we conduct benchmark EEVQE calculations by augmenting gates to form the branching MERA and executing the VQE for the ground-state search of Hamiltonians in Eqs. from the variational optimization of the MERA state to the VQE calculation with the branching-MERA ansatz.
In all three models, the optimization of the MERA state is almost completed up to 10 3 iterations; in fact, we confirm that this claim is valid by increasing the number of iterations up to 5 × 10 3 .After the quantum circuit encoding of the MERA state, we perform the VQE calculation, and the log 10 ∆ decreases by 0.15, 0.23, and 0.28 in the all-to-all coupled random Ising model, XYZ model, and Heisenberg model reflecting the entanglement augmentations.In particular, since the XYZ and Heisenberg models have more quantum fluctuation, the intensity of off-diagonal components may be comparable with its diagonal components, and the accuracy improvements by the entanglement augmentation are more significant than for the Ising model.
The right panel of Fig. 9 focuses on the VQE calculation and compares the behaviors of ∆ with two initial conditions: optimized MERA states and random branching MERA states.In the case of the all-to-all coupled random transverse-field Ising model, there is a moder- ately high initial-state dependence on the converged ∆, and the VQE with a random branching MERA as the initial state is trapped by local minima.It is a specific example where the EEVQE results show an advantage over the standard VQE.On the other hand, in the case of the all-to-all coupled random XYZ and Heisenberg models, there is no significant difference in the converged ∆ for both initial conditions.However, up to approximately 10 2 iterations of VQE, EEVQE can reduce ∆ compared to VQE, where the random branching MERA states are the initial states.This result suggests that the short-time VQE computation using the NISQ device is superior in our procedure, which is consistent with the current trend in quantum algorithm development.

V. CONCLUSION AND DISCUSSION
We have developed an entangled embedded VQE (EEVQE) method that uses a branching binary MERA with a one-dimensional entropic volume law as the ansatz for VQE calculations.In the EEVQE method, the binary MERA is optimized using a modified Evenbly-Vidal method on a classical computer and serves as the initial state for VQE.Unlike the original synergistic frameworks [10] that use the TN structure only for the initial state, EEVQE incorporates an entanglement augmentation topology based on the TN structure during gate addition.We investigated the performance of our method on the all-to-all coupled random transverse field Ising, XYZ, and Heisenberg models, and evaluated the random average ∆ of the relative error between the variational energy and the exact ground state energy under computational conditions (χ, N ) = (2, 16) to obtain the exact quantum circuit representation of the branching MERA embedded with the optimized MERA and the exact ground state of each model.
In the numerical simulation, we first examined the benchmark of the MERA optimization procedure on the classical computer using three methods: the Evenbly-Vidal algorithm, its modification, and the BFGS method as a fundamental knowledge of MERA for non-uniform systems.The results show that the BFGS method is superior to the other methods for all the models studied.This trend is considered a phenomenon specific to nonuniform systems since similar analysis in uniform systems [16] shows that sequential local updating by the Evenbly-Vidal method is superior to the gradient-based methods.This property becomes crucial when using the synergistic framework to study the ground states of quantum chemical calculations, which is one of the goals of VQE for all-to-all coupled inhomogeneous systems.
Second, we verified that the EEVQE method reduces log 10 ∆ by 0.15, 0.23, 0.28 in the all-to-all coupled random Ising, XYZ, and Heisenberg models, reflecting the entanglement augmentations while being free from the initialization problem of VQE calculations.Furthermore, we confirm that EEVQE can prevent becoming stuck in a local minimum in the all-to-all coupled random Ising model.At the same time, the VQE with a randomly initialized branching-MERA ansatz may become trapped in the minimum.Also, a comparison of VQE with random initialized branching MERA and branching MERA with embedded an optimized binary MERA as the initial state shows that the latter is superior to the former up to about 10 2 iterations of the update in VQE.In Appendix C, we also experimented with the EEVQE for the state that yields the entropic volume law to reinforce the validity of the previous results.These properties meet the requirements of the quantum variational algorithm using the NISQ device, which is to improve the accuracy of the best solution on a classical computer by a small amount of use of a quantum computer.
In the present study, since we applied one-dimensional branching MERA to a non-uniform system with all-toall coupled interactions, obtaining a highly accurate ∆ was not easy with χ = 2, even for the TN that satisfies the entropic volume law.It is not because of the barren plateau but because one-dimensional MERA and branching MERA with χ = 2 cannot represent the target ground state.Even branching MERA has a log-scale depth of the circuit, so it is considered not to be trapped in the barren plateau.Therefore, there remains to consider a scheme that can handle a χ > 2 in order to find out the limits.In this case, it is possible to apply automatic quantum circuit encoding (AQCE) [43] or other circuit decomposition techniques [36][37][38][39][40][41][42][43][44][45][46][47] for each disentangler, isometry, and top tensor of the TN to perform divide-and-conquer quantum circuit encoding for the entire TN.
As another future research, we can perform structural optimization of TN for each random interaction.In this case, ∆ with higher accuracy is expected to be easily obtained even for TN with the current minimum degree of freedom, χ = 2. Furthermore, we aim to integrate this approach into quantum circuits constructed using the aforementioned encoding methodologies to enhance the precision of it.We anticipate that these schemes can be applied to sets of quantum gates representing each isometric tensor individually, rather than to entire circuits.There have been reports on the structure optimization for TTNs that do not include loop structures in the TN [55][56][57].However, they cannot be extend directly to MERA, including loop structures, and further research is needed.
In addition, a deep MERA (DMERA) [58] has been reported as an algorithm that mixes the roles of the disentangler and isometry in MERA and improves the performance of MERA through quantum circuit representation.We can also introduce the branching degree of freedom in the DMERA.By examining the TN microscopically in terms of the product of two-qubit gates, it is more necessary to design a TN/quantum circuit structure that matches the geometrical aspect of the entanglement contained in the target state.
Other methods not treated in this study as optimizers for updating MERA include learning rate [32], automatic differentiation [59], and Riemannian optimization [29].It is interesting to verify whether these methods perform better than the BFGS method in the all-to-all coupled non-uniform models as one of the future issues to be addressed.
Finally, we are considering the development of quantum computers and their applications [60,61], and the verification of EEVQE on actual NiSQ devices should be addressed in future studies.).The blue and orange lines represent the case of the Evenbly-Vidal method and its modification, as introduced in Sec.III.2, respectively.The branching pattern was named according to the procedure shown in Fig. 2.
erations increases monotonically with respect to h, and log 10 ∆ for the branching MERA with 10 4 iterations decreasing monotonically with respect to h.As in the random XYZ/Heisenberg model case in Sec.IV, we confirm that the EEVQE solution outperformed the random initial VQE solution up to approximately 10 2 iterations in the right panels of Fig. 13.

II. 1 .
TN structures of 1D binary MERA and branching MERA

FIG. 1 .
FIG. 1. Schematic diagrams of (a) the binary MERA for the 16-site system and (b) isometric conditions of disentangler, isometry, and top tensors.

2 FIG. 2 .
FIG.2.Pattern diagram of the binary MERA to the branching MERA in an 8-site system.Since branches can be added at each layer, we name the diagram to branch0, branch1, branch2,. . .for each number of branches.

FIG. 3 .
FIG.3.Schematic diagram of eij with (i, j) = 4, 5, and the light-shaded region shows the causal cone structures.Diagrams outside the causal cone are trivially simplified to the identity operator by using isometric conditions.
FIG. 6.Schematic diagrams of (a) quantum circuit representations of disentangler, isometry, and top tensors, and (b) quantum circuit representation of the branching MERA for N = 8n with χ = 2 n consisting of 20 SU(2 2n ) 2n-qubit gates, where blue diagrams indicate the degree of freedom of (nonbranching) MERA.

FIG. 7 .
FIG. 7. Diagrammatic representation of an effective local two-body Hamiltonian h ′(ij)mn,m ′ n ′ for MERA as in Fig.6(b) with (i, j) = (3, 6) and χ = 2.Quantum circuits labeled with numbers 2 and 4 are omitted because they do not contribute to the effective local Hamiltonian.

FIG. 9 .
FIG. 9. Iteration dependence of the random-averaged relative error ∆ in the all-to-all coupled random (a) transverse-field Ising, (b) XYZ, and (c) Heisenberg models (left panels), and comparison of ∆ in the VQE calculation using initial states with embedded MERA states after TN optimizations and random branching MERA states (right panels).

FIG. 10 .
FIG. 10.Iteration dependence of ∆ for all-to-all coupled random (a) Ising and (b) XYZ systems with (χ, N ) = (2, 16).The blue and orange lines represent the case of the Evenbly-Vidal method and its modification, as introduced in Sec.III.2, respectively.The branching pattern was named according to the procedure shown in Fig.2.

FIG. 11 .
FIG. 11.The legend χ = [•, •, •] in Fig 11 (b) and 11 (c) means the χ = [χ0, χ1, χ2] as shown in Fig.11 (a).The χ0 is fixed as 2 in case of the spin system.The results displayed in the left and right panels correspond to the conventional method and our proposed method, respectively.

5 FIG. 13 .
FIG.13.The result of the EEVQE solution for the inhomogeneous Heisenberg model as C1 with N = 16.
FIG. 12. Representation of the rainbow state showing (−l, l) valence bonds above the central link.The nearest neighbor interactions are given in terms of α = exp(−h).