Efficient Adiabatic Preparation of Tensor Network States

We propose and study a specific adiabatic path to prepare those tensor network states that are unique ground states of few-body parent Hamiltonians in finite lattices, which include normal tensor network states, as well as other relevant nonnormal states. This path guarantees a gap for finite systems and allows for efficient numerical simulation. In one dimension, we numerically investigate the preparation of a family of states with varying correlation lengths and the one-dimensional Affleck-Kennedy-Lieb-Tasaki (AKLT) state and show that adiabatic preparation can be much faster than standard methods based on sequential preparation. We also apply the method to the two-dimensional AKLT state on the hexagonal lattice, for which no method based on sequential preparation is known, and show that it can be prepared very efficiently for relatively large lattices.

We propose and study a specific adiabatic path to prepare those tensor network states that are unique ground states of few-body parent Hamiltonians in finite lattices, which include normal tensor network states, as well as other relevant non-normal states. This path guarantees a gap for finite systems and allows for efficient numerical simulation. In 1D we numerically investigate the preparation of a family of states with varying correlation lengths and the 1D AKLT state and show that adiabatic preparation can be much faster than standard methods based on sequential preparation. We also apply the method to the 2D AKLT state on the hexagonal lattice for which no method based on sequential preparation is known, and show that it can be prepared very efficiently for relatively large lattices.
Matrix Product States (MPS) [1,2], and more generally, Projected Entangled-Pair States (PEPS) [3], capture the physical properties of systems obeying the entanglement area law [4]. PEPS contain a rich set of manybody states [5] such as the cluster state [6], toric codes [7], GHZ state [8] and W state [9] in quantum information, or the AKLT states [10,11], valence-bond states [12] and string net states [13] in condensed matter physics. There is thus increasing interest in finding ways of preparing them in quantum computers or quantum simulators, either for quantum information applications like computing [14], metrology [15], communication and networking [16], or as variational states for the study of manybody quantum systems [17].
MPS are most naturally prepared sequentially [18], which requires a time that scales linearly in the number of sites N . In higher dimensions, for PEPS, this is not possible in general [19]. However, certain subclasses of PEPS can be generated sequentially in linear time [20][21][22][23]. Sequential preparation has been used in various platforms to experimentally prepare MPS and PEPS [24][25][26][27].
Besides quantum circuits, adiabatic algorithms are also widely used to prepare many-body states on quantum devices [28]. By smoothly tuning the Hamiltonians along a gapped path that connects a trivial state to the target state, quasi-adiabatic evolution for a time T produces a state very close to the target state. Adiabatic algorithms have been proposed to prepare PEPS [29][30][31][32], and in particular, Ref. [31] proved that it is possible to prepare a generic family of them, so-called normal PEPS [1,2], in time T = O (polylogN ) with a specific method that switches on and off certain Hamiltonian terms adiabatically and provided there exists a gap along the whole path that is lower bounded by a constant. A method to compute such a lower bound based on semidefinite programming has been presented in Ref. [32]. While those methods provide rigorous proofs for the asymptotic limit N ≫ 1, it is not clear how they perform in practice, in particular for the intermediate sizes available in the near term. For such cases, there is no guarantee that they provide any advantage with respect to sequential methods.
In this paper, we propose a specific adiabatic path to prepare PEPS in any dimension that are unique ground states of local frustration-free Hamiltonians and analyze its performance. Our path guarantees the existence of a gap (for finite systems), and in contrast to Ref. [31,32] yields Hamiltonian with substantially smaller support for preparing the 2D AKLT state on the hexagonal lattice. Moreover, it extends to certain non-normal PEPS [33], which allows us to prepare the AKLT states in arbitrary geometries.
Since the ground states along our path are always PEPS, we are able to simulate relatively large systems, which we use to numerically determine the performance of the algorithm. In 1D, we consider the family of MPS introduced in Ref. [34], which allows us to investigate how the efficiency of the algorithm depends on correlation length. We also consider the paradigmatic 1D AKLT state. We obtain that for system sizes up to N = 5000, the preparation can be much more efficient than sequential preparation, with T ∼ polylog N in the regime we study [35]. In 2D, our adiabatic path overcomes several difficulties and allows us to simulate the adiabatic preparation of the 2D AKLT state on the hexagonal lattice up to N ∼ 10 × 10. Our results indicate that adiabatic preparation is very efficient also in higher dimensions.
PEPS.-PEPS can be built by applying local commuting operators to a product state of maximally entangled pairs in a lattice [3,36]. Let us consider a regular lattice denoted by a graph G, with edges E and sites V. The coordination number of site v ∈ V is n v , i.e. each site contains n v virtual qudits. Defining local operators {Q v } that map the D-level virtual qudits on site v ∈ V to a d-level physical site, the PEPS is expressed as [see Fig. 1(a) for 1D case and Fig. 1(c) for 2D hexagonal lat- Here D is the bond dimension of the PEPS, and d is the physical dimension. For instance, MPS can be viewed as 1D PEPS with n v = 2 virtual qudits per site [c.f. Fig. 1(a)]. The matrix representation of {Q v } in the bulk for MPS then reads The operators on the boundary each act on a single qudit [37]. By blocking neighboring sites, we can enlarge the physical dimension such that d ≥ D nv . In this case, without loss of generality, we can apply a polar decomposition to write {Q v } as positive-semidefinite operators with d = D nv , which holds for arbitrary PEPS up to a layer of local isometries [33]. [1,2]. If the operators obtained after blocking a finite number of sites are invertible, the PEPS is called normal.
In this paper, we aim to prepare a large class of PEPS that are unique ground states of local frustration-free Hamiltonians. This includes all normal (and thus all injective) PEPS [1,2], but also other relevant states like the AKLT states (possibly non-normal [33,36]), where a much simpler parent Hamiltonian is known [10,11]. In particular, we consider the following parent Hamiltonian [34,38] [c.f. Fig. 1(a)] where Π ker projects on the kernel of ρ e , which is the reduced density matrix of neighboring sites around the edge e ∈ E [39]. Note that ∥Π ker [ρ e ]∥ = 1, thus the time is unit-less in this paper. The parent Hamiltonian H [Eq. (3)] for injective PEPS has a unique ground state [38], which implies a nonzero gap that may depend on the system size N . Moreover, H for 1D injective MPS is guaranteed to be gapped also in the thermodynamic limit [1]. Finally, H for the AKLT states is equivalent to the known two-body parent Hamiltonian [10,11].
Examples.-We study two paradigm examples of PEPS in this paper. The first example is a family of MPS of bond dimension D = 2 [34]. In this case the graph G corresponds to a chain formed by N = 2N p qubits forming N p pairs [c.f. Fig. 1(a)]. After blocking each neighboring two sites, we arrive at the injective form of the MPS family for g ̸ = 0 (with d = D 2 = 4), where the matrices in Eq. (2) are given through ... (c) 2D PEPS on the hexagonal lattice. The green circles denote the operators {Qv}, and the connected dots denote maximally entangled virtual qudit pairs. Each term of the parent Hamiltonian H acts on neighboring sites (shown as red rectangles). The size of the lattice is Lx × Ly. In our numerics, we focus on the 2D AKLT state on the hexagonal lattice with cylinder boundary conditions (illustrated through gray lines).
The corresponding parent Hamiltonian [Eq. (3)] acts only on nearest neighbors, but with each site containing two qubits. We will study the preparation of states with g ∈ (−1, 0), which interpolates between the cluster state (g = −1) and the GHZ state (g = 0). For g < 0, the correlation length ξ of the MPS family can be obtained as [c.f. Fig. 1 Thus by tuning g, we can explore the effect of correlation length on the performance of the adiabatic algorithm. Note that g ∈ (−1, 0) already covers all states with g < 0, since the tensors {A(g)} in Eq. (4) can be mapped to {A(1/g)} by a gauge transformation [40]. The other example we consider is the 1D AKLT state of spin S = 1 and the 2D AKLT state of spin S = 3/2 in the hexagonal lattice [c.f. Fig. 1(c)] [10,11]. AKLT states can be formed by first having a product state of singlets consisting of virtual qubits that connect neighboring sites of the lattice, then projecting the virtual qubits at each site v to their symmetric subspace. AKLT states can be written as D = 2 PEPS [Eq. (1)], and we promote the virtual qubits into physical ones, such that the operators {Q v } are already positive-semidefinite without blocking [33].
Adiabatic algorithm.-We propose an adiabatic path parametrized by s, which connects a product state of maximally entangled pairs |ψ(0)⟩ ≡ e∈E |Φ + ⟩ e [41] to the target PEPS |ψ(1)⟩ ≡ |ψ⟩. We choose the instantaneous ground states |ψ(s)⟩ in this path to be always PEPS of bond dimension D, with [see Eq. (1)] In the following, we study the adiabatic preparation of these examples using our path.
Preparation of the MPS family.-First, we computed the minimal gap ∆ min during the adiabatic path [c.f. Eq. (6)] for the MPS family [c.f. Eq. (4)] in Fig. 1 One can see that ∆ min decreases as the correlation length ξ increases, which suggests that the adiabatic algorithm should perform better when the correlation length is smaller.
Assuming that the state we obtained after the evolution is |ϕ(T )⟩, its fidelity F compared to the target state |ψ(1)⟩ is |⟨ψ(1)|ϕ(T )⟩| 2 . For all fixed T ≥ 0, we numerically find that F decays exponentially with system size N [33], which allows us to define an error density κ(T ) that is independent of system size. Thus we can write where c(T ) is an error that comes from the boundaries of the system and is independent of N . This indicates that during the adiabatic dynamics, the errors in different regions of the chain change almost uniformly. The error density κ(T ) can be obtained by fitting the fidelity of preparing the same state of various system sizes N and a fixed time T using the scaling Eq. (7) [33]. In Fig. 2(a) we show κ(T ) for the MPS family as a function of T , and it features two regimes. When T is small, the dynamics is not adiabatic, and we see κ(T ) already starts to decay quickly. When T becomes larger, κ(T ) enters a regime of almost exponential decay, which we fit with The decay rate γ decreases with increasing correlation length ξ [see Fig. 2(c)], and the boundary term |c(T )| shows a similar behavior as κ(T ) [33]. Note that, due to the finite smoothness of the interpolation function s(t/T ) 1D we use, we expect that κ ∼ 1/poly(T ) when T → ∞ [33,[43][44][45][46]. One can extend the range of the exponential decay by making the interpolation function smoother, however at the expense of reducing the decay rate γ [33]. In Fig. 2(b), we show the dependence of time T required to prepare the given target state with fidelity F = 0.99 on system size N up to N ∼ 5000. The results agree well with the simple expression Eq. (7), which lead to T ∼ polylog N in this regime (which is the relevant regime experimentally). We also compare the adiabatic preparation to sequential preparation [18], which we assume takes a time T = N [see Fig. 2(b)]. One sees that the adiabatic algorithm outperforms the sequential preparation method in terms of the preparation time T when the system size N is larger than a threshold value N c (ξ) [47]. We find numerically that [c.f. Fig. 2(c)] N c almost grows exponentially with ξ, which indicates that when the correlation length is smaller than a fraction of the system size, the adiabatic algorithm prepares the MPS family [c.f. Eq. (4)] faster than the sequential method.
Finally, we show the time T to prepare states of system size N = 100 as a function of the correlation length ξ in Fig. 2(d). We observe T ∝ ξ (T ∝ ξ 2 ) when ξ is smaller (larger) than the length of the lattice site, which is ξ = 2 since each site contains two qubits. This shows that the size of each lattice site sets another length scale in the system, and one can also see such behavior for the decay rate γ shown Fig. 2(b).
Preparation of AKLT states.-Now we study the preparation of the 1D and 2D AKLT states using our adiabatic path. In 1D, it has been proposed to prepare the AKLT state sequentially [18], dissipatively [48], using measurements [49], or by parallelly fusing multiple AKLT chains [48], which has the best known preparation time T = O(log 2 N ). In Fig. 2(a,b) we show the results for adiabatic preparation of the 1D AKLT state using our adiabatic path. As expected, T ∼ polylog N up to N = 5000.
Preparation of the 2D AKLT states are much less explored. For the case of hexagonal lattice, this state have a gapped parent Hamiltonian [50,51], and there is indirect evidence suggesting that this state can be adiabatically prepared with T = O(N ) [52]. The general protocol [31,32] predicts the preparation time T = O (polylogN ) when N ≫ 1, but it faces the following challenges: First, the construction of the parent Hamiltonian there requires the target PEPS to be normal, which does not work for non-normal PEPS such as the 2D AKLT state on the square lattice [36], or leads to a Hamiltonian for the 2D AKLT state on the hexagonal lattice that acts on large clusters [53], making it difficult to implement in current devices. More importantly, it is difficult to simulate an adiabatic evolution in 2D since the time cost of classical simulation algorithms typically has heavy dependence on the bond dimension of the un-derlying PEPS [54].
Our adiabatic path overcomes the above problems (partially because we promote the virtual qubits to physical ones [33]). The Hamiltonian [c.f. Eq. (3)] along the whole path is two-body and gapped (note that each site contain 3 qubits). Moreover, the instantaneous ground state [c.f. Eq. (6)] is always a PEPS of bond dimension D = 2 [55].
We classically simulate the preparation of the 2D AKLT state on the hexagonal lattice with cylinder boundary condition, and N ≡ L x × L y sites [c.f. Fig. 1(c) In Fig. 3 we show the preparation time T needed to reach a fidelity F = 0.99 for different system sizes N . In the case of N = L × L, since each hexagon is of size 2 × 1 [c.f. Fig. 1(c)], for even or odd L the boundary affects T differently. We see that overall T is practically short (T ∼ 10), and increases only mildly with system size N . We also show T for preparing this state of lattice size N = L x × L by fixing different L x = 4, 6, 8, and observe similar behavior. In particular, the lattice geometry does not strongly affect the preparation time T , as it takes a similar T to prepare the state on a 4 × 20 lattice or 9 × 9 lattice.
We expect that T will still be practically short when we increase the system size further than that shown in Fig. 3. Moreover, the general nature of the proposed path suggest that this method may be used to efficiently prepare a large variety of (high-dimensional) PEPS. For example, in the SI [33] we provide numerical evidence that the 2D AKLT state on the square lattice can be prepared with this method.
Outlook.-We have proposed and studied a specific adiabatic path to prepare a large family of MPS and PEPS, which applies to the normal PEPS and other relevant PEPS like the AKLT states.
It is worth checking if the 2D AKLT state of N ≫ 100 can still be efficiently prepared using quantum devices, and exploring the performance of this adiabatic path to prepare other (potentially non-normal [36]) PEPS. By studying the gap during the adiabatic path in the thermodynamic limit [32,56], it is possible to probe the asymptotic behavior of the adiabatic algorithm and further improve its performance with adiabatic ramp rate that adapts to the magnitude of the gap. Moreover, in Ref. [31] it is shown that adiabatic preparation can also be implemented efficiently on digital quantum computers, which also applies to our results. It is also important to design efficient physical realization of the proposed adiabatic path, which require engineering fewbody Hamiltonians. Hamiltonain engineering can be realized in various platforms like the superconducting qubits [57], ion traps [58] and Rydberg atomic arrays [59], which potentially allows to realize large classes of Hamiltonians [60][61][62]. Finally, one can study the effect of noise on the adiabatic state preparation. In the presence of noise, we expect the adiabatic method provides an even bigger advantage over the sequential methods for preparing short-range correlated states, since more error accumulates during the sequential preparation (which takes a longer time).
Acknowledgments.-We thank Norbert Schuch and Yilun Yang for their insightful discussions. The research is part of the Munich Quantum Valley, which is supported by the Bavarian state government with funds from the Hightech Agenda Bayern Plus. We acknowledge funding from the German Federal Ministry of Education and Research (BMBF) through EQUAHUMO (Grant No. 13N16066) within the funding program quantum technologies -from basic research to market, and the European Union's Horizon 2020 research and innovation program under Grant No. 899354 (FET Open Su-perQuLAN). The numerical calculations were performed using the ITensor Library [63].
Note added.-During completing of this manuscript we became aware of a related protocol that probabilistically creates the 1D and 2D AKLT states using constdepth circuits and post-selection, with a success rate that exponentially decays with the system size N [64].  (8)] in a wide range of evolution time T . Since in we need to reach larger system size N ∼ 5000 in 1D compared to N ∼ 100 in 2D, we choose the function s(t, T )1D to be more smooth than s(t, T )2D at t = 0 and t = T . [47] Note that the sequential method prepares the state deterministically, while the adiabatic method always prepares the state approximately (with high fidelity), which is enough for practical purposes. Moreover, there could be a constant overhead to implement the Hamiltonian dynamics using quantum circuits, and this does not affect our statement that there will be a regime (N ≫ 1) where the adiabatic preparation will be more efficient.
v is a positive-semidefinite operator of dimension D nv × D nv . Thus one can rewrite the Eq.(1) in the main text as Thus up to a single layer of local isometries v∈V U v , one can take the operators {Q v } to be positive-semidefinite. Note that, to physically implement the isometries, we require d ≥ D nv . This can be realized by blocking neighboring sites, or promoting the virtual qub(d)its on each site into physical ones in the case of AKLT states. Therefore, during the preparation of the target PEPS with operators {Q v }, first, we use the adiabatic path to prepare the state with operators {Q ′ v }, then apply a single layer of local isometries {U v } at the end of the adiabatic preparation. As it takes a constant time T iso = O(1) to implement the isometries, in the main text we focus on the adiabatic evolution time T , and the total state preparation time is T tot = T + T iso .
Note that, for AKLT states, one does not need to implement local isometries after the adiabatic preparation (i.e. T iso = 0). Thus the adiabatic preparation time T shown in Fig.2(b) and Fig.3 in the main text for AKLT states are the total preparation time of the state.

A. 1D TEBD simulation
We use the TEBD algorithm [S1] to study the one-dimensional adiabatic dynamics. Given a local Hamiltonian H(t) = Np i=1 h i (t) with t = nτ , we approximate each step U t = e −iH(t)τ by following Trotter-Suzuki approximation [S2] and apply such unitaries to MPS to evolve the state. The precision of this algorithm can be tuned by changing the Trotter step τ and a cutoff parameter δ (in ITensor [S3]) that controls the precision of applying a unitary to MPS. We typically observe the convergence of results with a Trotter time step τ = 0.04 and δ = 1 × 10 −10 for the data shown in this paper. Note we are able to use TEBD because the state during the adiabatic evolution [c.f. Eq.(6) in the main text] remains close to the ground state of the parent Hamiltonian, which is an MPS of bond dimension D = 2.

B. Additional numerical results
In Fig.1(b) in the main text, we show the minimal gap ∆ min during the adiabatic path for preparing states in the MPS family with different g. This is obtained using the density-matrix renormalization group (DMRG) method [S4]. We present more data on the gap in Fig. S1(a). In Fig. S1(b) we show that the gap ∆ min computed with different system sizes does not change with N .
and we see the collapse of data of different N for all time T . Therefore, Fig. S1(c,d) demonstrate the correctness of Eq. (7) in the main text. In Fig. S1(e) we plot c(T ) for various states in the MPS family and the AKLT state. As c(T ) represents errors on the boundary of the chain, it also shows an exponential decay with T (thus |c(T )| → 0), similar to the behavior of κ(T ) in Fig.2(a) in the main text. To understand this behavior, here we employ an exact diagonalization (ED) method to study the evolution of the same adiabatic dynamics [c.f. Eq.(6) in the main text] with a small system size N = 8, which can enter very small error regimes that κ(T ) ≪ 1. Thanks to the scaling between fidelity F and N [c.f. Eq.(7) in the main text], the error density obtained for such a small system size could already provide a qualitative understanding of the behavior that also applies to larger systems.
Previous work has shown that, for a few-body system, the location where the behavior of κ(T ) transition from the exponential decay to power-law decay is mainly determined by the smoothness of s(t/T ) at the beginning (t = 0) and end (t = T ) of the adiabatic evolution [S8]. So here we use the incomplete Beta functions s(t/T ) = θ k (t/T ) in the adiabatic path [c.f. Eq.(6) in the main text], which are given through with B λ (a, b) ≡ λ 0 y a−1 (1 − y) b−1 dy. This function satisfies k-th order smoothness at λ = 0 and λ = 1, namely d n θ k dλ n λ=0,1 = 0, ∀n < k.
Thus by controlling the order k, one can systematically control the smoothness of the interpolation function.  One can see an exponential decay of the infidelity in a short time and then transition to a power-law decay regime at a late time. By increasing k, the exponential decay region can be increased, and the late-time power-law decay behaviors are with α ≈ 2k +2 shown in the inset of Fig. S2(a). This behavior agrees with that in Ref. [S7, S8]. Also in Fig. S2(b), we show the evolution of the infidelity for preparing the MPS family with various g using a fixed interpolation function (k = 3) [c.f. Eq. (S5)]. For states with different correlation lengths, we see the transition from the exponential decay regime to the power-law decay regime generally happens when the infidelity is very small. Therefore, by choosing a s(t/T ) that is smooth enough, we expect the exponential decay of the error density κ(T ) can last to practically relevant (and potentially very large) system sizes. Note that, the above analysis does not mean that we should choose a too smooth s(t/T ). As Ref. [S8] has shown, the exponential decay rate γ [c.f. Eq. (8) in the main text] decreases with the smoothness of s(t/T ). Moreover, when choosing s(t/T ) to be infinitely smooth at t = 0 and t = T , it has been shown for the case of the Gevrey class 1 + α with α > 0 [S9] that the error density κ Gev (T ) [c.f. Eq.(7) in the main text] evolves as [S10] κ Gev (T ) which is a slower-than-exponential decay since α > 0. Therefore, in practice, one needs to choose the level of smoothness of s(t/T ) appropriately to minimize the overall infidelity.

B. 2D TDVP simulation
We order the 2D hexagonal lattice sites as a one-dimensional chain as shown in Fig. S3(a). In this case, the 2D local Hamiltonian [c.f. Eq.(3) in the main text] translates to a long-range 1D Hamiltonian, and we use the time-dependent variational principle (TDVP) method [S14, S15] to study the time evolution of the system since it potentially has better performance than TEBD method for long-range interacting systems [S16]. Also, note that the bond dimension D of this effective MPS grows exponentially with the linear dimension L as D = O[exp(L)], therefore our numerical method does not scale up to very large system sizes. The adiabatic path with a ground state as D = 2 PEPS [c.f. Eq.(6) in the main text] is thus crucial for simulating the adiabatic dynamics of this system up to a system size N ∼ 10 × 10 [c.f. Fig.3 in the main text]. This is demonstrated in Fig. S3(b), that one can see the evolution of the bond dimension D during the adiabatic evolution with various total evolution time T . When T is small, the bond dimension increases significantly during the evolution. By increasing T , the dynamics become more adiabatic, thus we observe the overall bond dimension gets significantly reduced.
C. Preparation of the 2D AKLT state on the square lattice To explicitly demonstrate that our method can prepare non-normal high-dimensional PEPS (in contrast to the adiabatic path in Ref. [S10, S17]), in Fig. S4(b) we show the evolution time T needed to prepare the AKLT state of spin