Efficient Simulation of 1D Long-Range Interacting Systems at Any Temperature
Rakesh Achutha Donghoon Kim Yusuke Kimura Tomotaka Kuwahara





















one



















 






=260mm















 
 
 


























#1










 














































 











    ⟨⟩#1√( #1) Tr tr
        
        theoremTheoremsubtheoremSubtheoremlemmaLemmacorol[lemma]Corollaryassump[lemma]AssumptiondefinitionDefinition
        
        
        prop[theorem]Proposition
        
        claim[subtheorem]Claim
        
        conj[theorem]Conjecture
        
        
        
        
        
        
        
        
        #1#2
        
        
        
        
        
        achutha.rakesh.mat20@itbhu.ac.indonghoon.kim@riken.jpyusuke.kimura.cs@riken.jptomotaka.kuwahara@riken.jp^1
        Analytical quantum complexity RIKEN Hakubi Research Team, RIKEN Center for Quantum Computing (RQC), Wako, Saitama 351-0198, Japan
        ^2
        Department of Computer Science and Engineering, Indian Institute of Technology (Banaras Hindu University), Varanasi, 221005, India
        ^3
        PRESTO, Japan Science and Technology (JST), Kawaguchi, Saitama 332-0012, Japan^4
        RIKEN Cluster for Pioneering Research (CPR), Wako, Saitama 351-0198, Japan
        § ABSTRACT 
        We introduce a method that ensures efficient computation of one-dimensional quantum systems with long-range interactions across all temperatures. Our algorithm operates within a quasi-polynomial runtime for inverse temperatures up to β=(ln(n)). At the core of our approach is the Density Matrix Renormalization Group algorithm, which typically does not guarantee efficiency. We have created a new truncation scheme for the matrix product operator of the quantum Gibbs states, which allows us to control the error analytically. Additionally, our method is applied to simulate the time evolution of systems with long-range interactions, achieving significantly better precision than that offered by the Lieb-Robinson bound.
        Efficient Simulation of 1D Long-Range Interacting Systems at Any Temperature
            Tomotaka Kuwahara^1,3,4
            September 9, 2024
        ============================================================================
        
        Introduction.— 
        Unraveling the complex patterns of quantum many-body systems presents one of the most significant challenges in modern physics. Among these challenges, characterizing the thermal equilibrium quantum state (also known as the quantum Gibbs state) at zero and non-zero temperatures stands as a central target. These states are generally difficult to manage in high-dimensional systems as mentioned in previous studies <cit.>. However, a variety of techniques have been developed to address their properties in one-dimensional systems <cit.>. 
        Notably, techniques such as the transfer matrix method used in classical Ising models <cit.> and its generalization to 1D quantum systems <cit.> have been well-recognized. More recently, the Density Matrix Renormalization Group (DMRG) algorithm has emerged as a promising approach <cit.>, enabling the numerical construction of the matrix product operator (MPO) that encapsulates all information about the equilibrium states <cit.>. Despite its promise, providing a rigorous justification for the precision guarantee in the DMRG algorithm remains a critical open problem at both zero and non-zero temperatures <cit.>. Nevertheless, several methods have been proposed to develop a `DMRG-type' algorithm to construct the MPO with guaranteed accuracy and time efficiency <cit.>.
        
        
        
        In one-dimensional systems at non-zero temperatures, the construction of the MPO for short-range interacting systems has received extensive attention. The pioneering method, based on the cluster expansion technique by Molnar et al. <cit.>, introduced a polynomial-time algorithm. This algorithm works with a runtime scaling as (n/ϵ)^𝒪(β), where n is the system size, to approximate the 1D quantum Gibbs state by the MPO within an error ϵ. Subsequent advancements have refined the error bounds, with the leading-edge approach achieving a runtime of e^𝒪̃(β) + 𝒪̃(√(βln(n/ϵ))) <cit.>. This technique utilizes the imaginary-time adaptation of the Haah-Hastings-Kothari-Low (HHKL) decomposition, originally devised for real-time evolution <cit.>.
        
        
        
        Extending from the study of short-range interactions in one-dimensional (1D) systems, we explore the realm of long-range interactions, characterized by power-law decay. These interactions impart high-dimensional traits to systems, manifesting behaviors not observed in short-range interacting systems. Notably, phenomena such as phase transitions, typically associated with higher dimensions, can occur even in 1D systems with long-range interactions <cit.>. 
        Additionally, long-range interactions are increasingly relevant in contemporary many-body physics experiments, such as those involving ultracold atomic systems <cit.>, highlighting their practical significance <cit.>.
        
        
        
        
        Beyond the realm of short-range interacting systems, the complexity analysis of long-range interactions presents substantial challenges, especially at low temperatures. Most existing methods are tailored specifically for short-range interactions, and their direct application to long-range systems often results in a loss of efficiency <cit.>. While the cluster expansion method proves effective at high temperatures <cit.>, providing an efficient algorithm to analyze the properties of the quantum Gibbs state in systems with long-range interactions, it falls short at low temperatures. Notably, it does not offer an approximation by the MPO, leaving the computational approach at low temperatures as an unresolved issue.
        
        
        
        In this letter, we introduce an efficiency-guaranteed method for constructing the MPO based on the DMRG algorithm. Our method achieves a runtime of:
        
            e^𝒪(βln^3(n/ϵ))  (general cases),
        
        and for Hamiltonians that are 2-local, as defined in Eq. (<ref>) below, the runtime improves to:
        
            e^𝒪̃(βln^2(n/ϵ))  (2-local cases).
        
        These results demonstrate a quasi-polynomial time complexity under conditions where β=polylog(n), a regime where the quantum Gibbs states typically approximate the non-critical ground state <cit.>. 
        Furthermore, while our primary focus has been on constructing MPOs for quantum Gibbs states, our approach is equally applicable to real-time evolutions. 
        To our knowledge, this is the first rigorous justification of a classical simulation method with guaranteed efficiency for managing the real/imaginary-time dynamics of long-range interacting Hamiltonians.
        
        
        
        In terms of technical advancements, the traditional DMRG algorithm faced challenges in accurately estimating the truncation errors of bond dimensions during iterative steps <cit.> (see also Ref. <cit.>). To address this, we developed a new algorithm that captures the core mechanics, as illustrated in Fig. <ref>. Our approach comprises two main steps: (i) constructing the MPO at high temperatures with a precisely controlled error, using arbitrary Schatten p-norms, and (ii) merging high-temperature quantum Gibbs states into their low-temperature counterparts. The success of the first step allows us to apply techniques from Refs. <cit.> in the second step.
        Hence, the primary challenge lies in the initial construction of the MPO. We begin by constructing the quantum Gibbs state in independent small blocks and iteratively merge them into larger ones. 
        This approach retains the spirit of the original DMRG algorithm developed by White <cit.>, focusing on local interactions and iterative improvements.
        A key innovation of our algorithm is the efficient control of the precision in the merging process, quantified by the general Schatten p-norm (see Proposition 2 below).
        
        
        
        
        
        
        
        Setup.— We consider a one-dimensional quantum system composed of n qudits, each in a d-dimensional Hilbert space. The total set of sites in this system is denoted by Λ, i.e., Λ = {1, 2, 3, …, n}. We define a k-local Hamiltonian H for this system by:
        
            H = ∑_|Z| ≤ k h_Z,   max_i ∈Λ∑_Z : Z ∋ ih_Z≤ g,
        
        where Z ⊂Λ, and |Z| represents the cardinality of Z, meaning the number of sites involved in each interaction term h_Z. The term ⋯ denotes the operator norm, which is the maximum singular value of the operator. For any arbitrary subset L ⊆Λ, the subset Hamiltonian H_L is defined to include only those interaction terms involving sites in L, specifically:
        
            H_L = ∑_Z:Z ⊂ L h_Z.
        
        
        
        
        
        
        To elucidate the assumptions underpinning our analysis, let us consider a decomposition of the total system Λ into subsets A and B such that Λ = A ⊔ B. Here, subset A is defined as [i, i']∩Λ with 1 ≤ i < i' ≤ n. We assume that the norm of the boundary interaction between A and B, denoted by ∂ h_A, is finite:
        
            ∂ h_A≤g̃,   ∀ A,   ∂ h_A := H - (H_A + H_B),
        
        where g̃ represents an 𝒪(1) constant. This condition is more general than the assumption of power-law decay of interactions. If we specifically consider the power-law decay condition in the form:
        
            ∑_Z: Z ∋{i, i'}h_Z≤J/|i - i'|^α i≠ i',
        
        then the condition in (<ref>) is satisfied as long as α > 2
        (see Supplementary materials <cit.>).
        
        
        
        
        Throughout the paper, we analyze the quantum Gibbs state as follows:
        
            ρ_β := 1/Z_βe^-β H,    Z_β := tr(e^-β H).
        
        We aim to approximate the Gibbs state ρ_β by a MPO, which is generally described as follows:
        
            M = ∑_s_1, s_2, …, s_n = 1 
             s_1', s_2', …, s_n' = 1^d    tr(M_1^[s_1, s_1'] M_2^[s_2, s_2']⋯ M_n^[s_n, s_n']) 
               |s_1, s_2, …, s_n⟩⟨ s_1', s_2', …, s_n'|,
        
        where each {M_j^[s_j, s_j']}_j,s_j,s_j' is a D × D matrix, with D referred to as the bond dimension.
        As measures of the approximation, we often use the Schatten p-norm as follows:
        O_p :=[(|O|^p)]^1/p
        with O an arbitrary operator and |O| = √(O^†O). For p=1, it is equivalent to the trace norm, and for p=∞, it corresponds to the operator norm.
        
        
        
        
        We denote the bond dimension of the MPO form of the Hamiltonian by D_H, which is at most of n^k d^k in general <cit.>. 
        In particular, if we restrict ourselves to a 2-local Hamiltonian in the form of:
        
            H    = ∑_i < i'1/|i - i'|^α∑_ξ, ξ' = 1^d^2 J_ξ,ξ' P_i,ξ⊗ P_i',ξ',
        
        where {P_i,ξ}_ξ=1^d^2 denotes the operator bases at site i and J_ξ,ξ' are the corresponding coefficients, the Hamiltonian can be approximated by an MPO with bond dimension D_H = c_H ln^2(n/ϵ), where c_H = 𝒪(1), up to an error of ϵ (see Supplementary material <cit.>).
        
        We note that any subset Hamiltonian H_L is also described by an MPO with the bond dimension D_H of the global Hamiltonian.
        
        
        
         Main Result.— 
        Our task is to develop a time-efficient method to approximate the quantum Gibbs state using an MPO with a fixed bond dimension. 
        The main result concerning the time-efficient construction of the MPO is as follows:
        
         Theorem. 
        For an arbitrary β, we can efficiently compute the MPO M_β that approximates the unnormalized Gibbs state up to an error ϵ:
        
            e^-β H - M_β_p ≤ϵe^-β H_p,   ∀ p.
        
        The bond dimension of M_β and the time complexity are given by 
        e^C_0 βln^3(n/ϵ),
        where C_0 is an 𝒪(1) constant. 
        Typically, we consider ϵ = 1/(n), which leads to a time complexity of e^𝒪(β) ln^3(n). 
        If we consider the 2-local Hamiltonian in Eq. (<ref>), 
        the time complexity qualitatively improves to e^C_0 ln^2(n/ϵ) lnln (n/ϵ).
        From the theorem, we can efficiently simulate the quantum Gibbs state in long-range interacting systems <cit.>.
        
        
         
        
        
        
        In the following, we detail the explicit calculation steps for constructing the MPO M_β, which is also illustrated in Fig. <ref>:
        
          * Initial State Preparation: We begin with a product state of the quantum Gibbs states of small blocks, each containing 2 sites, at an inverse temperature β_0. This temperature β_0 is chosen to be smaller than 1/(24gk^2). The exact MPO representation of the initial state has a bond dimension at most d^2, where d is the dimension of the local Hilbert space.
            
            
          * Merging Operation: Adjacent quantum Gibbs states on blocks are merged using a merging operator Ψ. This operator is approximated by Ψ̃ using truncated expansions as described in Eq. (<ref>), up to m_0-th order, where m_0 = ln(n/ϵ). Consequently, the bond dimension of the MPO form of Ψ̃ becomes D_H^m_0, with D_H representing the bond dimension of the Hamiltonian.
            
            
          * Approximation of High-Temperature Gibbs State: By repeating the merging processes log_2(n) times, we achieve the global quantum Gibbs state e^-β_0 H. This high-temperature quantum Gibbs state is approximated by the MPO M_β_0, which has a bond dimension D_H^m_0 ln(n).
        
            
            
          * Approximation of Low-Temperature Gibbs State: The high-temperature quantum Gibbs state is connected (β/β_0) times to yield the MPO (M_β_0)^β/β_0, serving as the approximation of the desired quantum Gibbs state e^-β H.
        
        In the next section, we will focus on analytically estimating the precision error and the time complexity resulting from these approximations.
        
        
        Finally, we discuss the preparation of quantum Gibbs states on a quantum computer. 
        In conclusion, the time complexity (or circuit depth) for this preparation is also given by the quasi-polynomial form, i.e., e^C_0 βln^3(n/ϵ). 
        To achieve this, we first purify the Gibbs state into the form of a thermofield double state <cit.>, which can also be well-approximated by a Matrix Product State (MPS) with the same bond dimension as that of M_β <cit.>. 
        For the quantum circuit representation of the MPS, we refer to Ref. <cit.>, where it is generally shown that an arbitrary MPS with bond dimension χ can be constructed using a circuit with depth n ×poly(χ). By setting χ = e^C_0 βln^3(n/ϵ) for the quantum Gibbs state's preparation, we achieve the desired circuit depth.
        
        
        
        Proof of main theorem.— 
        We here provide the outline of the proof, and we defer the details to Supplementary materials <cit.>.
        We follow Steps 1 to 4 above and estimate the precision error depending on the bond dimension. 
        We denote the Gibbs state of the block L_s^(q) as e^-β_0 H_s^(q), where H_s^(q) serves as an abbreviation for the subset Hamiltonian H_L_s^(q) = ∑_Z : Z ⊂ L_s^(q) h_Z on block L_s^(q). Under this notation, e^-β_0 H_1^(q_0) = e^-β_0 H represents the Gibbs state of the entire system, encompassing all interactions within the system.
        
        We then estimate the error arising from the merging operations. 
        In general, we define a merging operator Ψ that connects two subsets, A and B, as follows:
        
            Ψ = e^-β_0 H_AB e^β_0 (H_A + H_B),
        
        which implies the relation e^-β_0 H_AB = Ψ e^-β_0 H_A e^-β_0 H_B. Here, Ψ facilitates the merging of the Gibbs states of A and B to form the combined Gibbs state on the larger block AB, as illustrated in Fig. <ref>(a). 
        In this context, H_AB = H_A + H_B + ∂ h_A, where ∂ h_A represents the boundary interactions between systems A and B, as defined in Eq. (<ref>). 
        
        To construct an MPO with a small bond dimension for Ψ, we first approximate Ψ using a polynomial expansion in terms of β_0. 
        However, the independent approximations of e^-β_0 H_AB and e^β_0 (H_A + H_B) require polynomial degrees of at least (β_0H_AB)^1/2 and (β_0H_A + β_0H_B)^1/2, respectively <cit.>. This results in a demand for sub-exponentially large bond dimensions for an accurate approximation.
        
        
        
        
        The key idea here is that by combining e^-β_0 H_AB and e^β_0 (H_A + H_B) the truncation order is significantly reduced for a good approximation of Ψ. 
        We adopt the following polynomial approximation of the merging operator (<ref>): 
        
            Ψ  = ∑^m_0_m=0β_0^m∑_s_1+s_2=m (-1)^s_1H_AB^s_1(H_A+H_B)^s_2/s_1!s_2!,
        
        where we apply the Taylor expansion for each of e^-β_0 H_AB and e^β_0 (H_A + H_B) and truncate the terms with β_0^m for m>m_0. 
        Then, we prove the following precision error:  
        Proposition 1.  If β_0 is smaller than 1/(24gk^2), we achieve the approximation of 
        
            Ψ - Ψ≤δ_0  m_0≥log_2 c_0/δ_0 ,
        
        with c_0=e^g̃/(6gk^2) which is an 1 constant. 
        
        
        
        To derive the MPO form of Ψ, we first create individual MPOs for each term in the expansion (<ref>). These individual MPOs are then combined to obtain the MPO representation of Ψ.
        Given the MPO form of the Hamiltonians H_AB and H_A + H_B, the bond dimension required for the term H_AB^s_1(H_A + H_B)^s_2 is at most D_H^s_1D_H^s_2 = D_H^m, where s_1 + s_2 = m. The total bond dimension required to represent all terms up to m_0 in the expansion is calculated as:
        
            ∑_m=0^m_0 (m+1) D_H^m ≤ (m_0+1)^2 D_H^m_0 =: D̃_δ_0.
        
        Thus, we construct the MPO with the bond dimension D̃_δ_0 to approximate the merging operator Ψ for any arbitrary adjacent subsets A and B.
        
        
        After constructing the MPO approximation of the merging operator, we proceed to find the MPO for the combined Gibbs state. 
        We in general consider the merging of e^-β_0 H_2s-1^(q-1) and e^- β_0 H_2s^(q-1), which is given by 
        e^-β_0 H_s^(q)= Ψ_s^(q-1) e^-β_0 H_2s-1^(q-1) e^-β_0 H_2s^(q-1). 
        Our purpose is to estimate the approximation by Ψ̃_s^(q-1) M_2s-1^(q-1) M_2s^(q-1) when approximate MPOs M_2s-1^(q-1) and M_2s^(q-1) are given for e^-β_0H_2s-1^(q-1) and e^-β_0 H_2s^(q-1), respectively. 
        For simplicity, we denote Ψ_s^(q-1) as Ψ and Ψ̃_s^(q-1)  as Ψ̃, by omitting the indices s and q.  
         
        Proposition 2. Let M_2s - 1^(q - 1) and M_2s^(q - 1) represent the MPOs approximating e^-β_0 H_2s - 1^(q-1) and e^-β_0 H_2s^(q-1), respectively, such that
        
            e^-β_0 H_i^(q-1) - M_i^(q-1)_p ≤ϵ_q-1e^-β_0 H_i^(q-1)_p
        ∀ p with i = 2s-1,2s.
        Using the approximate merging operator Ψ with Ψ - Ψ≤δ_0, the merged quantum Gibbs state e^-β_0 H_s^(q) can be approximated as:
        
            e^-β_0 H_s^(q) - M_s^(q)_p ≤   ϵ_qe^-β_0 H_s^(q)_p,
        
        with ϵ_q = a_2 δ_0 + a_1 ϵ_q-1,  
        where the MPO M_s^(q) is constructed as M_s^(q) = Ψ M_2s-1^(q-1) M_2s^(q-1), and a_1 and a_2 are 𝒪(1) constants with a_1 > 1. 
        
        
        This shows that ϵ_q and ϵ_q-1 are related by a linear equation, leading to a recursive relation on q. By solving this recursive relation, we find:
        ϵ_q  ≤ a_2 δ_0 q a_1^q-2,
        using the fact that ϵ_1 = 0, because the MPO description of the initial state is exact.
        We then obtain the MPO approximation M_β_0 of the Gibbs state e^-β_0 H, which results in the following error ∀ p:
        
            e^-β_0 H - M_β_0_p ≤ϵ'_0 e^-β_0 H_p
        
        with ϵ'_0:= a_2 δ_0 q_0 a_1^q_0-2.
        The MPO M_β_0 = M_1^(q_0) is constructed through q_0 approximate merging processes Ψ, starting from the Gibbs states in layer 1. Consequently, the bond dimension increases by a factor of D̃_δ_0^q_0, where D̃_δ_0 was defined as the bond dimension of Ψ in Eq. (<ref>).
        If we consider the error Ψ - Ψ≤δ_0 = ϵ/poly(n), then the bond dimension becomes D̃_δ_0^q_0 = D_H^ln^2 (n/ϵ), as q_0 = log_2(n) and D̃_δ_0 = D_H^ln (1/δ_0).
        
        
        
        
        
        Finally, we use the MPO M_β_0 for the high-temperature Gibbs state e^-β_0 H to approximate the target low-temperature Gibbs state e^-β H, as illustrated in Fig. <ref>(b). We combine the (β/β_0) MPOs M_β_0 to approximate the low-temperature Gibbs state, resulting in:
        (M_β_0)^β/β_0≈ e^-β H.
        Consequently, the bond dimension is multiplied (β/β_0) times, leading to a final required bond dimension of the order D_H^βln^2 (n/ϵ).
        
        To derive the error bound, we use the inequality from <cit.>, which states the following: for fixed positive integers p_1 and p_2, if 
        e^-β_0H - M_β_0_p_1p_2≤ϵ' e^-β_0H_p_1p_2,
        then 
        e^-p_1β_0H - M_β_0^p_1_p_2≤ (3e/2)p_1ϵ' e^-p_1β_0H_p_2.
        We have already derived the MPO approximation for a general Schatten p-norm in (<ref>). 
        Therefore, we can obtain  
        
            e^-β H - (M_β_0)^β/β_0_p ≤ϵ_βe^-β H_p,
        
        where 
        ϵ_β = 5(β/β_0)  ϵ'_0.
        Since a_1^q_0 = poly(n) with q_0 = log_2(n), the choice of δ_0 = ϵ/poly(n) yields the main precision bound in Eq. (<ref>), where we defined M_β = (M_β_0)^β/β_0 and assume β<n without loss of generality <cit.>.
        
        By estimating each of the algorithm steps, we can see that the runtime is also upper-bounded by D_H^βln^2 (n/ϵ) (see Supplementary materials <cit.>).
        We thus prove the main theorem.  □Extension to Real-Time Evolution–
        In our method, we did not rely on the fact that the inverse temperature β is a real number. Therefore, it is possible to generalize our approach to cases where β is a complex number. By considering the case where β = it, we can address the simulation of real-time evolution. 
        
        Given an MPO M_t that approximates e^-iHt, we can estimate the error as:
        e^-iHt|Ψ⟩ - M_t|Ψ⟩≤e^-iHt - M_t_∞ for an arbitrary quantum state |Ψ⟩. 
        By following our algorithm, we can construct an MPO M_t for e^-iHt that satisfies the error bound e^-iHt - M_t_∞≤ϵe^-iHt_∞=ϵ. The bond dimension and time complexity required are given by a quasi-polynomial form of e^tln^3 (n/ϵ). 
        Compared to methods based on the Lieb-Robinson bound, our approach significantly reduces the time complexity. Assuming the Lieb-Robinson bound as [O_i(t), O_j]≤ t^η/r^α for α > 2 <cit.>, the time cost to calculate the average value of a time-evolved local observable O_i(t) scales as e^(t^η/ϵ)^1/(α-1).
        
        
        
        Summary and outlook.—
        In this letter, we have presented an algorithm that achieves quasi-polynomial time complexity in approximating the quantum Gibbs state of 1D long-range interacting systems using MPOs. This was accomplished by adopting a DMRG-type method, as depicted in Fig. <ref>. 
        We note that achieving quasi-polynomial complexity is likely the best attainable result in general <cit.>. 
        However, there is hope to achieve polynomial time complexity for the 2-local cases, as described in Eq. (<ref>). Identifying the optimal time complexity for simulating long-range interacting systems remains an important open question.
        
        Another intriguing avenue for exploration is whether our method can be extended to cases where the power-law decay is slower than r^-2 (α < 2), where the assumption of finite boundary interactions in Eq. (<ref>) may no longer hold. In specific cases, this condition can still be recovered. For instance, in fermion systems with long-range hopping and short-range fermion-fermion interactions, the condition in Eq. (<ref>) holds for α > 3/2 <cit.>. Another interesting class is the non-critical quantum Gibbs state, where power-law clustering is satisfied. In such systems, the MPO approximation is expected to hold for α > 1, since the mutual information for any bipartition obeys the area law <cit.>.
        
        
        
         
         
         
         
         
         
         
        
         
         
        
        All the authors acknowledge the Hakubi projects of RIKEN. 
        T. K. was supported by JST PRESTO (Grant No.
        JPMJPR2116), ERATO (Grant No. JPMJER2302),
        and JSPS Grants-in-Aid for Scientific Research (No.
        JP23H01099, JP24H00071), Japan.
        Y. K. was supported by the JSPS Grant-in-Aid for Scientific Research (No. JP24K06909). 
        This research was conducted during the first author’s internship at RIKEN, which was supervised by the last author.
        
        
        
        section0
        
         Supplementary Material for "Efficient Simulation of 1D Long-Range Interacting Systems at Any Temperature"
        
        
        Rakesh Achutha^1,2, Donghoon Kim^1, Yusuke Kimura^1 and Tomotaka Kuwahara^1,3,4
        
        ^1 Analytical quantum complexity RIKEN Hakubi Research Team, RIKEN Center for Quantum Computing (RQC), Wako, Saitama 351-0198, Japan
        ^2 Department of Computer Science and Engineering, Indian Institute of Technology (Banaras Hindu University), Varanasi, 221005, India
        ^3 PRESTO, Japan Science and Technology (JST), Kawaguchi, Saitama 332-0012, Japan
        ^4 RIKEN Cluster for Pioneering Research (CPR), Wako, Saitama 351-0198, Japan
        
        
        § SEVERAL BASIC STATEMENTS
        
        
         §.§ Operator norm of the boundary interaction for long-range  interacting systems when α > 2
        
        
        We assume in the main text that the boundary interaction is bounded:
        
            ∂ h_A≤g,     ∀ A ⊂Λ,
        
        where the boundary interaction ∂ h_A is defined as the interaction acting on both subsystems A and B, derived by subtracting the Hamiltonians H_A and H_B, which act exclusively on A and B, from the total Hamiltonian,
        
            ∂ h_A := ∑_Z: Z ∩ A ≠∅, Z ∩ B ≠∅ h_Z = H - (H_A + H_B).
        
        We demonstrate that this assumption is inherently satisfied for long-range interactions with a power law decay of α > 2:
        
        Considering the power-law decay of interactions in the form 
        
            ∑_Z: Z ∋i,i'h_Z≤J/|i-i'|^α i≠ i',
        
        we prove that g̃ in Eq. (<ref>) is 1 as long as α>2.
        Proof of Lemma <ref>.
        Let us consider two adjacent subsystems, A and B. 
        From the definition of ∂ h_A, we obtain 
        
            ∂ h_A   ≤∑_Z: Z ∩ A ≠∅, Z ∩ B ≠∅h_Z≤∑_i ∈ A∑_i' ∈ B∑_Z: Z ∋ i, i'h_Z≤∑_i ∈ A∑_i' ∈ BJ/|i-i'|^α≤∑_i=1^∞∑_j=1^∞J/(i+j - 1)^α = ζ(α - 1),
        
        where ζ(s) is the Riemann zeta function. For α > 2, we can choose g̃ = ζ(α - 1) = 𝒪(1).
        
        This completes the proof. 	□
        
         §.§ MPO approximation of the 2-local Hamiltonians
        
        
        We generally obtain the MPO forms of k-local Hamiltonians which have the bond dimensions at most of n^k d^k.  
        In this section, we show that the bond dimension is improved in the specific cases of the 2-local Hamiltonians.
        
        For a 2-local Hamiltonian of the form
        
            H    = ∑_1 ≤ i < i' ≤ n1/|i - i'|^α∑_ξ, ξ' = 1^d^2 J_ξ,ξ' P_i, ξ⊗ P_i', ξ',
        
        where {P_i, ξ}_ξ = 1^d^2 are operator bases on site i and J_ξ,ξ' are constants, there exists a Hamiltonian H̃ that approximates H with an error bounded by H - H̃≤ϵ_H, and admits an MPO representation with bond dimension
        
            D_H̃ = c_Hln^2(n/ϵ_H),
        
        where c_H is an 1 constant. 
         
         Remark.
        Here, the MPO form for the Hamiltonian H̃ is not exactly the same as that of the original Hamiltonian H.
        We hence need to take into account the error between H and H̃ for the simulation of the quantum Gibbs states (see Sec. <ref>).
        
         
        Proof of Lemma <ref>. 
        We approximate r^-α as an exponential series in accordance with Ref. <cit.>, as follows
        We cannot directly employ Ref. <cit.> since it treats only the regime r<1.:
        
            r^-α - ∑_s∈ℤ e^α s x e^-e^s xr ≤ϵ r≥ 1,  
               x= 2π/ln(3)+αln[1/cos(1)] + ln(1/ϵ).
        
        We then consider the finite truncation as s∈ℤ_m= [-m,m]. 
        We obtain 
        
            ∑_s∈ℤ e^α s x e^-e^s xr - ∑_s∈ℤ_m e^α s x e^-e^s xr ≤∑_s<-m e^α s x +  ∑_s>m e^α s x -e^s x.
        
        We assume that m satisfies
        
            α m x -e^m x≤ -α m x    ⟶   e^m x≥ 2α m x      ⟶    
            m x ≥ -  ProductLog-1/2α
               ⟶    m ≥2/xln(2α) ,
        
        where ProductLog(x) is equivalent to the Lambert W function W(x), defined by x=W(x)e^W(x).  
        
        Under the assumption, we reduce the inequality (<ref>) to 
        
            ∑_s∈ℤ e^α s x e^-e^s xr - ∑_s∈ℤ_m e^α s x e^-e^s xr ≤2∑_s>m e^-α s x = 2e^-α x m/e^α x-1 =  2/e^α x-1ϵ/2α^2α
        
        with 
        
            m=⌈2/xln(2α /ϵ) ⌉ .
        
        Under the above choice, we obtain 
        
            r^-α  - ∑_s∈ℤ_m e^α s x e^-e^s xr ≤ϵ+  2/e^α x-1ϵ/2α^2α≤ζϵ ,
        
        where ζ is an 1 constant for α≥ 2. 
        
        Based on this approximation, we define the Hamiltonian H̃ as:
        
            H̃   =  ∑_1 ≤ i < i' ≤ n∑_s∈ℤ_m e^α s x   e^-e^s x|i - i'| ∑_ξ, ξ' = 1^d^2 J_ξ,ξ' P_i, ξ⊗ P_i', ξ'.
        
        We then obtain 
        
            H-    H̃≤ζϵ∑_1 ≤ i < i' ≤ n∑_ξ, ξ' = 1^d^2 J_ξ,ξ'P_i, ξ⊗ P_i', ξ'≤J̅ n^2ζϵ ,
        
        with 
        
            J̅ := ∑_ξ, ξ' = 1^d^2 |J_ξ,ξ'|  ,
        
        where we set the operator bases such that P_i, ξ=1 for ∀ i and ∀ξ. 
         
        To ensure an approximation error of ϵ_H for H-    H̃, we let J̅ n^2ζϵ=ϵ_H 
        (or ϵ=ϵ_H/(J̅ζ n^2)), which implies 
        
            m= ⌈2/xln(2α /ϵ) ⌉≤ c'_Hln^2(n/ϵ_H),
        
        where c'_H is a constant of order 𝒪(1).
        
        Moreover, it has been shown that each Hamiltonian of the form ∑_1 ≤ i < i' ≤ n e^α s x  e^-e^s x |i - i'| ∑_ξ, ξ' = 1^d^2 J_ξ,ξ' P_i, ξ⊗ P_i', ξ' 
        with exponential decay admits an MPO representation with a bond dimension of 𝒪(d^2) <cit.>. Therefore, the MPO representation of H̃ can be constructed by combining the MPOs of each term, resulting in a total bond dimension of:
        
            D_H̃   = c_Hln^2(n/ϵ_H)
        
        for d=1, where c_H is an 𝒪(1) constant.
        This completes the proof. □
        
        § PROOF OF MAIN PROPOSITIONS AND TECHNICAL DETAILS
        
        
         §.§ Proof of Proposition 1 in the main text
        
        
        We now present the proof of Proposition 1 in the main text.
        For the convenience of readers, we recall the setup. 
        We define a merging operator Ψ that connects two subsets A and B as follows:
        
            Ψ = e^-β_0 H_AB e^β_0 (H_A + H_B),      H_AB = H_A + H_B + ∂ h_A,
        
        which leads to the relation
        
            e^-β_0 H_AB = Ψ  e^-β_0 H_A e^-β_0 H_B.
        
        Consider the polynomial expansion of the merging operator Ψ in powers of β_0, given by:
        
            Ψ = e^-β_0 H_AB e^β_0 (H_A + H_B)   = e^-β_0 H_AB e^β_0 (H_AB- ∂ h_A)
               = ∑_m=0^∞β_0^m ∑_s_1 + s_2 = m1/s_1! s_2! (-H_AB)^s_1 (H_AB - ∂ h_A)^s_2.
        
        We truncate the above series to m_0 terms and denote the resulting truncated series as Ψ̃:
        
            Ψ̃   = ∑_m=0^m_0β_0^m ∑_s_1 + s_2 = m1/s_1! s_2! (-H_AB)^s_1 (H_AB - ∂ h_A)^s_2.
        
        Then, if β_0 is smaller than 1/24gk^2, we achieve the approximation
        
            Ψ - Ψ̃≤δ_0  m_0 ≥log_2 c_0/δ_0 ,
        
        where c_0 =e^g̃/(6gk^2) is a constant of 𝒪(1).
        Proof of Proposition <ref>. 
        We derive the approximation of the truncation of the expansion. 
        To this end, we first reformulate each term in the expansion (<ref>)  by employing an alternative expansion, specifically using the interaction picture:
        
            e^β_0 (H_AB - ∂ h_A) = e^β_0 H_AB𝒯( e^-∫_0^β_0∂ h_x   dx),
        
        where we define ∂ h_x := e^- x H_AB∂ h_A e^x H_AB and 𝒯 is the time-ordering operator. We then obtain
        
            Ψ   = e^-β_0 H_AB e^β_0 (H_AB- ∂ h_A)
               = 𝒯( e^-∫_0^β_0∂ h_x   dx) = ∑_s = 0^∞ (-1)^s∫_0^β_0 dx_1∫_0^x_1 dx_2⋯∫_0^x_s-1 dx_s ∂ h_x_1∂ h_x_2⋯∂ h_x_s.
        
        By applying the decomposition
        
            ∂ h_x = ∑_p=0^∞(-x)^p/p!ad_H_AB^p (∂ h_A) =: ∑_p=0^∞(-x)^p/p!∂ h^(p),
        
        the m-th order term in the expansion of Ψ in powers of β_0, which includes the contribution from β_0^m in Eq. (<ref>), is expressed as
        
            Ψ_0   := 1, 
            Ψ_m    := ∑_s=1^m ∑_p_1 + p_2 + ⋯ + p_s = m - s (-1)^s ∫_0^β_0 dx_1∫_0^x_1 dx_2⋯∫_0^x_s-1 dx_s (-x_1)^p_1/p_1!⋯(-x_s)^p_s/p_s! ∂ h^(p_1)⋯∂ h^(p_s)
        
        for ∀ m ≥ 1.
        As both Eq. (<ref>) and Eq. (<ref>) are expansions in powers of β_0, the term Ψ_m defined above corresponds to the mth order term in Eq. (<ref>):
        
            Ψ_m =  β_0^m ∑_s_1 + s_2 = m1/s_1! s_2! (-H_AB)^s_1 (H_AB - ∂ h_A)^s_2,
        
        which we aim to upper-bound. 
        
        The norm of Ψ_m is upper-bounded as
        
            Ψ_m   ≤∑_s=1^m β_0^m-s∑_p_1 + p_2 + ⋯ + p_s = m - s∫_0^β_0 dx_1  ∫_0^x_1 dx_2⋯∫_0^x_s-1 dx_s ∏_j=1^s 1/p_j!∂ h^(p_j)
               ≤β_0^m ∑_s=1^m 1/s!∑_p_1 + p_2 + ⋯ + p_s = m - s∏_j=1^s 1/p_j!∂ h^(p_j).
        
        From Theorem 2.1 in <cit.>, if H is a k-local and g-extensive Hamiltonian and A is a r-local operator, then
        
            [H, A]≤ 6gkr A.
        
        Since H_AB and ∂ h_A are k-local Hamiltonians, ad_H_AB^p-1 (∂ h_A) is at most kp-local. Using the inequality (<ref>), we get
        
            ∂ h^(p) = ad_H_AB^p (∂ h_A)   ≤ 6gk(kp) ad_H_AB^p-1 (∂ h_A) .
        
        By iterating the above process for ad_H_AB^p' (∂ h_A) with p' ≤ p-2 and using (<ref>), we obtain
        
            ∂ h^(p)≤ (6gk^2)^p p! g̃ =: C^p p! g̃,
        
        where we define C := 6gk^2. Substituting Eq. (<ref>) into Eq. (<ref>), we get
        
            Ψ_m   ≤β_0^m ∑_s=1^m C^m-sg̃^s/s!∑_p_1 + p_2 + ⋯ + p_s = m - s 1 = β_0^m ∑_s=1^m C^m-sg̃^s/s!(m-ss).
        
        Using (m-ss) = m - 1s≤ 2^m ,
        
            Ψ_m   ≤ (2Cβ_0)^m ∑_s=1^m1/s!g̃/C^s≤ (2Cβ_0)^m e^g̃/C.
        
        Given the assumption β_0 ≤1/24gk^2 = 1/4C, it follows that
        
            Ψ_m≤ 2^-m e^g̃/C,
        
        which exhibits exponential decay with respect to m. 
        Thus, the truncation error resulting from approximating Ψ by Ψ̃ is bounded by
        
            Ψ - Ψ̃≤∑_m=m_0+1^∞Ψ_m≤∑_m = m_0 + 1^∞ 2^-m e^g̃/C = e^g̃/C 2^-m_0.
        
        Letting c_0 =e^g̃/C = e^g̃ / (6gk^2), to achieve the desired error Ψ - Ψ̃≤δ_0, it is sufficient to choose m_0 such that Ψ - Ψ̃≤ c_0 2^-m_0≤δ_0, i.e.,
        
            m_0 ≥log_2 c_0/δ_0.
        
        This completes the proof of Proposition <ref>. □
        
        
        
        
        
        
        
        
        
         §.§ Proof of Proposition 2 in the main text
        
        
        In this subsection, we prove Proposition 2 in the main text, which concerns the approximation of the merged quantum Gibbs state.
        
        
        Let M_2s - 1^(q - 1) and M_2s^(q - 1) represent the MPOs approximating e^-β_0 H_2s - 1^(q-1) and e^-β_0 H_2s^(q-1), respectively, such that
        
            e^-β_0 H_2s-1^(q-1) - M_2s-1^(q-1)_p    ≤ϵ_q-1e^-β_0 H_2s-1^(q-1)_p, 
            e^-β_0 H_2s^(q-1) - M_2s^(q-1)_p    ≤ϵ_q-1e^-β_0 H_2s^(q-1)_p,
        
        for an arbitrary Schatten p-norm. 
        Using the merging operator Ψ̃_s^(q-1), which provides an approximation to Ψ_s^(q-1) with the error bounded by Ψ_s^(q-1) - Ψ̃_s^(q-1)≤δ_0,
        the merged quantum Gibbs state e^-β_0 H_s^(q) can be approximated by the newly constructed MPO M_s^(q) based on M_s^(q) = Ψ_s^(q-1) M_2s-1^(q-1) M_2s^(q-1) as
        
            e^-β_0 H_s^(q) - M_s^(q)_p ≤ϵ_qe^-β_0 H_s^(q)_p,
        
        with the error
        
            ϵ_q = a_2 δ_0 + a_1 ϵ_q-1.
        
        Here, a_1 and a_2 are 𝒪(1) constants with a_1 > 1, defined as
        
            a_1    := 12 e^g̃/(4gk^2)  ,      a_2 := 2 e^g̃/(24gk^2) .
        Proof of Proposition <ref>. The merging operator is defined as Ψ_s^(q-1) = e^-β_0 H_s^(q) e^β_0 (H_2s - 1^(q-1) + H_2s^(q-1)), and the approximate merging operator Ψ̃_s^(q-1) is given by Eq. (<ref>).
        With the MPO defined as M_s^(q) = Ψ̃_s^(q-1) M_2s-1^(q-1) M_2s^(q-1), we proceed with the following inequality:
        
            e^-β_0 H_s^(q) - M_s^(q)_p    = Ψ_s^(q-1) e^-β_0 H_2s-1^(q-1) e^-β_0 H_2s^(q-1) - Ψ̃_s^(q-1) M_2s-1^(q-1) M_2s^(q-1)_p 
               ≤Ψ_s^(q-1) - Ψ̃_s^(q-1) e^-β_0 H_2s-1^(q-1) e^-β_0 H_2s^(q-1)_p 
                    + Ψ̃_s^(q-1)e^-β_0 H_2s-1^(q-1) e^-β_0 H_2s^(q-1)  - M_2s-1^(q-1) M_2s^(q-1)_p 
               ≤Ψ_s^(q-1) - Ψ̃_s^(q-1)_∞e^-β_0 H_2s-1^(q-1) e^-β_0 H_2s^(q-1)_p 
                    + Ψ̃_s^(q-1)_∞e^-β_0 H_2s-1^(q-1) e^-β_0 H_2s^(q-1) - M_2s-1^(q-1) M_2s^(q-1)_p.
        
        From Ref. <cit.>, we can prove for arbitrary operators  and 
            e^+ -e^_p ≤ e^e^_p .
        
        By applying the inequality and substituting -β_0 H_s^(q) for  and β_0 H_s^(q) - β_0 H_2s^(q-1) -β_0 H_2s-1^(q-1) for , we obtain ≤g̃ from assumption (<ref>) and
        
            e^-β_0 H_2s-1^(q-1) e^-β_0 H_2s^(q-1) - e^-β_0 H_s^(q)_p ≤ e^β_0 g̃e^-β_0 H_s^(q)_p.
        
        The triangle inequality leads to
        
            e^-β_0 H_2s-1^(q-1) e^-β_0 H_2s^(q-1)_p 
            ≤e^β_0 g̃+1e^-β_0 H_s^(q)_p
            ≤ 2 e^β_0 g̃e^-β_0 H_s^(q)_p .
        
        From the above inequality and the condition Ψ_s^(q-1) - Ψ̃_s^(q-1)_p ≤δ_0, we reduce inequality (<ref>) to 
        
            e^-β_0 H_s^(q) - M_s^(q)_p 
                   ≤  2δ_0  e^β_0 g̃e^-β_0 H_s^(q)_p + Ψ̃_s^(q-1)_∞e^-β_0 H_2s-1^(q-1) e^-β_0 H_2s^(q-1) - M_2s-1^(q-1) M_2s^(q-1)_p.
        
        
        
        To estimate the second term on the RHS of (<ref>),
        we use the following mathematical lemma:
        
         
        If two operators O_1 and O_2 act on disjoint supports, with their respective approximations Õ_1 and Õ_2 defined on the supports of O_1 and O_2, satisfying 
        
            O_1 - Õ_1_p    ≤ϵO_1_p, 
            O_2 - Õ_2_p    ≤ϵO_2_p,
        
        with ϵ<1, we then obtain 
        
            O_1 O_2 - Õ_1 Õ_2_p ≤ 3 ϵO_1 O_2_p.
        Proof of Lemma <ref>. When O_1 and O_2 operate on disjoint supports, we have O_1 O_2_p = O_1_p O_2_p. By leveraging this fact along with Eqs. (<ref>) and (<ref>), we derive
        
            O_1 O_2 - Õ_1 Õ_2_p    ≤O_1 (O_2 - Õ_2)_p +  (O_1 - Õ_1) Õ_2_p 
               = O_1_p O_2 - Õ_2_p + O_1 - Õ_1_p Õ_2_p 
               ≤ϵO_1_p O_2_p + ϵO_1_p Õ_2_p.
        
        The relation Õ_2_p≤O_2_p + O_2 - Õ_2_p≤ (1 + ϵ) O_2_p and ϵ^2≤ϵ leads to
        
            O_1 O_2 - Õ_1 Õ_2_p    ≤ϵO_1_p O_2_p + ϵ (1 + ϵ) O_1_p O_2_p  = (2 ϵ + ϵ^2) O_1_p O_2_p ≤ 3 ϵO_1_p O_2_p.
        
        This completes the proof of Lemma 3. 
        	□ 
        
        
        Since e^-β_0 H_2s-1^(q-1) and e^-β_0 H_2s^(q-1) act on disjoint supports (i.e., L_2s-1^(q-1)∩ L_2s^(q-1) = ∅), applying the Lemma 3 yields:
        
            e^-β_0 H_2s-1^(q-1) e^-β_0 H_2s^(q-1) - M_2s-1^(q-1) M_2s^(q-1)_p    ≤ 3 ϵ_q-1e^-β_0 H_2s-1^(q-1) e^-β_0 H_2s^(q-1)_p 
               ≤ 6 ϵ_q-1 e^β_0 g̃e^-β_0 H_s^(q)_p,
        
        where we use the upper bound (<ref>) in the second inequality. 
        
        We then upper-bound the norm of  Ψ̃_s^(q-1), which is the truncated series up to m_0 terms in the expansion of Ψ_s^(q-1) as in Eq. (<ref>).
        By using the upper bound (<ref>), which gives Ψ_m≤ 2^-m e^g̃/C, we obtain  
        
            Ψ̃_s^(q-1)≤∑_m = 0^m_0Ψ_m≤∑_m=0^m_0 2^-m e^g̃/C≤ 2 e^g̃/C=2 e^g̃/(6gk^2) .
        
        We thus derive the upper bound 
        
            Ψ̃_s^(q-1)_∞e^-β_0 H_2s-1^(q-1) e^-β_0 H_2s^(q-1) - M_2s-1^(q-1) M_2s^(q-1)_p
            ≤  12 e^g̃/(6gk^2)+β_0 g̃ϵ_q-1e^-β_0 H_s^(q)_p .
        
        By combining the above inequality with (<ref>), we have 
        
            e^-β_0 H_s^(q) - M_s^(q)_p 
                   ≤ 2δ_0  e^β_0 g̃ + 12 e^g̃/(6gk^2)+β_0 g̃ϵ_q-1e^-β_0 H_s^(q)_p  
               ≤ 2δ_0  e^g̃/(24gk^2) + 12 e^g̃/(4gk^2)ϵ_q-1e^-β_0 H_s^(q)_p   ,
        
        where we use β_0 ≤ 1/(24gk^2).  
        Thus, we establish the desired inequality (<ref>), along with Eqs. (<ref>) and (<ref>).
        This completes the proof of Proposition 2. 
        	□
        
         §.§ Explicit estimation of the bond dimension
        
        
        From Proposition 2, we see that ϵ_q follows a recursive relation, with ϵ_1 = 0 since the blocks in layer 1 have exact MPOs. Using ϵ_1=0, the recursive relation yields the following relations:
        
            ϵ_q    = a_2 δ_0 + a_1 ϵ_q-1 = a_2 δ_0 (1 + a_1 + a_1^2 + ⋯ + a_1^q-2) + a_1^q-1ϵ_1 = a_2 δ_0 (a_1^q-1 - 1/a_1 - 1) + 0  ≤ a_2 δ_0 q a_1^q-2.
        
        Given that e^-β_0 H_1^(q_0) = e^-β_0 H and M_β_0 := M_1^(q_0), Proposition 2 implies
        
            e^-β_0 H - M_β_0_p    ≤ϵ_q_0e^-β_0 H_p.
        
        Thus, using the inequality (<ref>) and q_0= log_2 (n), we obtain the MPO approximation of the high-temperature Gibbs state e^-β_0 H given by:
        
            e^-β_0 H - M_β_0_p    ≤ a_2 δ_0 q_0 a_1^q_0-2e^-β_0 H_p
                 ≤ a_2 δ_0  n^log_2 (2a_1)e^-β_0 H_p.
        
        
        Starting from the approximate MPO M_β_0 for the high-temperature Gibbs state, we construct the MPO for the low-temperature Gibbs state. By concatenating the MPO M_β_0(β / β_0) times, we obtain (M_β_0)^β / β_0≈ e^-β H.
        To derive the error bound for the low-temperature Gibbs state, we begin by considering the approximation in the (pQ)-norm, where Q = β / β_0. Given that we have already established the MPO approximation for any general Schatten p-norm, it follows that
        
            e^-β_0 H - M_β_0_pQ   ≤ a_2 δ_0  n^log_2 (2a_1)e^-β_0 H_pQ.
        
        We now use the inequality from <cit.>, which states:
        For fixed positive integers p_1 and p_2, if e^-β_0 H - M_β_0_p_1 p_2≤ϵ' e^-β_0 H_p_1 p_2, then 
        
            e^-p_1 β_0 H - M_β_0^p_1_p_2≤3e/2 p_1 ϵ' e^-p_1 β_0 H_p_2.
        
        Applying this to our case, with p_1 =Q= β / β_0 and p_2 = p, we get:
        
            e^-β_0 Q H - (M_β_0)^Q_p    ≤3e/2 Q a_2 δ_0  n^log_2(2a_1)e^-β_0 Q H_p ≤ 5 (β/β_0) a_2 δ_0  n^log_2(2a_1)e^-β_0 Q H_p.
        
        Now we choose
        
            δ_0 = β_0ϵ/5β a_2  n^log_2 (2a_1).
        
        Since a_2, β_0, and a_1 are 𝒪(1) constants, we have δ_0 = ϵ/poly(n) under the assumption β < n. By selecting δ_0 as in Eq. (<ref>), we obtain the final desired inequality,
        
            e^-β H - M_β_p    ≤ϵe^-β H_p.
        
        
        In the following, we denote the bond dimensions of the MPOs M_β_0 and M_β by D_β_0 and D_β, respectively.
        The MPO M_β_0 is constructed by (q_0-1) approximate merging processes of Ψ̃ starting from the layer 1.  
        The bond dimension of the MPO for the Gibbs states in layer 1 is d^2. This bond dimension increases by D̃_δ_0 with the application of the merging operator Ψ̃ at each subsequent layer, where D̃_δ_0 is defined as the bond dimension of Ψ̃ in Eq. (<ref>).
        Hence, the final bond dimension after (q_0-1) mergings becomes d^2 D̃_δ_0^q_0-1.
        From Eq. (<ref>), we have D̃_δ_0 = (m_0 + 1)^2 D_H^m_0≤ D_H^2m_0, where m_0 is chosen as 
        
            m_0 = log_2(c_0 / δ_0)
            = log_25 c_0 β a_2  n^log_2 (2a_1)/β_0ϵ
         according to Proposition 1, with a_0 > 1 ensuring that m_0 is an integer.
        This yields
        
            D_β_0   = d^2 D̃_δ_0^q_0-1≤D̃_δ_0^q_0≤ D_H^2 m_0 q_0.
        
        Using Eq. (<ref>)  and q_0= log_2 (n), we get
        
            2 m_0 q_0    ≤ 2log_25 c_0 β a_2  n^log_2 (2a_1)/β_0ϵlog_2 (n) 
               = 2/ln^2(2)ln5 c_0 β a_2  n^log_2 (2a_1)/β_0ϵln (n) 
               ≤ 5 ln5 c_0a_2/β_0+ log_2 (4a_1)ln(n/ϵ)ln (n) 
               ≤ 5 ln5 c_0a_2/β_0+ log_2 (4a_1)ln^2 (n/ϵ)=:b_1 ln^2 (n/ϵ)  (n/ϵ)≥ e,
        
        with
        
            b_1:=  5 ln5 c_0a_2/β_0+ log_2 (4a_1) ,
        
        where we use β < n and define an 𝒪(1). 
        We recall the definitions of  β_0, c_0, a_1 and a_2 as follows:
        
            β_0 ≤1/24gk^2,   c_0 =e^g̃/(6gk^2) ,   a_1:= 12 e^g̃/(4gk^2)  ,    a_2 := 2 e^g̃/(24gk^2) .
        
        Hence, the order of D_β_0 becomes D_H^𝒪(ln^2(n/ϵ)), which directly leads to the following bond dimension D_β for the MPO M_β corresponding to the low-temperature Gibbs state:
        
            D_β   = D_β_0^β / β_0≤ D_H^b_1(β/β_0 ln^2(n/ϵ))≤ D_H^(b_1/β_0)βln^2(n/ϵ) (n/ϵ)≥ e.
        
        Hence the order of the bond dimension of M_β is D_H^βln^2 (n / ϵ) since b_1, β_0 are 𝒪(1) constants. 
        
        
        
         §.§ Estimation of the time complexity
        
        
        
        We finally estimate the time complexity of constructing the MPO M_β. We begin with the construction of M_β_0. Given two MPOs M_1 and M_2 with bond dimensions D_1 and D_2, respectively, the time cost to calculate M_1 M_2 is given by n (D_1 D_2)^2 d^3.
        The total number of merging processes is given by 
            2^q_0-1 + 2^q_0-2 + ⋯ + 1 = 2^q_0 - 1 < n,
         and each merging process requires a time cost of at most n D_β^4 d^3, 
        where we use D_β_0 is smaller than D_β, i.e.,  the final bond dimension of the MPO M_β. 
        Thus, the total time cost for constructing M_β_0 is bounded by n^2 D_β^4 d^3.
        
        For each multiplication of M_β_0 and M_β_0^s with s ≤β/β_0, the time cost is also upper-bounded by n D_β^4 d^3. Therefore, the time cost to calculate (M_β_0)^β/β_0 is at most n (β/β_0) D_β^4 d^3.
        In total, the time cost to construct M_β is upper-bounded by
        
            (n^2 +  nβ/β_0) D_β^4 d^3 = D_H^βln^2 (n/ϵ).
        
        By applying D_H = n^k d^k in general cases and D_H = c_H ln^2(n/ϵ) in 2-local cases, we confirm the desired time complexity.
        
        
        
        
        
        
        
        
        
        
        
        
        
        
        
         §.§ MPO construction for the 2-local cases
        
        In the 2-local case, the original total Hamiltonian H can be approximated by H̃ as an exponential series, following Sec. <ref>, such that H - H̃≤ϵ_H for any arbitrary ϵ_H. The MPO for the approximate Hamiltonian H̃ is then constructed according to Sec <ref>, satisfying
        
            e^-βH̃ - M̃_β_p ≤ϵ/3e^-βH̃_p,
        
        where we assume ϵ≤1. 
        Here the bond dimension for M̃_β is given by D_H^(b_1/β_0)βln^2(3n/ϵ) with D_H=c_H ln^2(n/ϵ_H).
        
        We estimate the error between e^-β H and e^-βH using the Ref. <cit.>, which states that we can prove for arbitrary operators  and 
            e^+ -e^_p ≤ e^·e^_p .
        
        By applying the inequality and substituting -β H for  and βH - β H for , we get
        
            e^-β H - e^-βH_p    ≤ e^βH-HβH-H·e^-β H_p ≤ e^βϵ_Hβϵ_H e^-β H_p ≤ 3βϵ_He^-β H_p,
        
        where we choose ϵ_H such that ϵ_H≤ 1/β. 
        For a given β, if we choose ϵ_H ≤ϵ/ (6β)[From this inequality, the condition βϵ_H≤ 1 is satisfied because ϵ≤ 1.], we obtain
        
            e^-β H - e^-βH_p    ≤ϵ/2e^-β H_p.
        
        Using inequalities (<ref>) and (<ref>),
        
            e^-βH- M_β_p ≤ϵ/3e^-βH̃_p ≤ϵ/3e^-βH-e^-β H_p + e^-β H_p   ≤ϵ/3ϵ/2 + 1e^-β H_p
               ≤ϵ/2e^-β H_p ,
        
        where we use ϵ≤ 1 to derive ϵ/2 + 1≤ 3/2. 
        
        Therefore, using inequalities (<ref>) and (<ref>), we finally obtain the following error
        
            e^-β H- M_β_p ≤e^-β H-e^-βH_p + e^-βH-M_β_p ≤ϵe^-β H_p.
        
        This gives us the approximation MPO for the final Gibbs state e^-β H for the 2-local case.
