Groundstate fidelity phase diagram of the fully anisotropic two-leg spin-1/2 XXZ ladder

The fully anisotropic two-leg spin-1/2 $XXZ$ ladder model is studied in terms of an algorithm based on the tensor network representation of quantum many-body states as an adaptation of projected entangled pair states to the geometry of translationally invariant infinite-size quantum spin ladders. The tensor network algorithm provides an effective method to generate the groundstate wave function, which allows computation of the groundstate fidelity per lattice site, a universal marker to detect phase transitions in quantum many-body systems. The groundstate fidelity is used in conjunction with local order and string order parameters to systematically map out the groundstate phase diagram of the ladder model. The phase diagram exhibits a rich diversity of quantum phases. These are the ferromagnetic, stripe ferromagnetic, rung singlet, rung triplet, N$\rm \acute{e}$el, stripe N$\rm \acute{e}$el and Haldane phases, along with the two $XY$ phases $XY1$ and $XY2$.


Introduction
Spin ladder systems have attracted considerable attention from both experimentalists and theoreticians alike [1,2]. One of the most striking features of the simple spin- 1 2 Heisenberg ladder is that the spin excitations are gapful (gapless) when the number of legs is even (odd) [3]. Spin ladder systems in general represent a particularly interesting class of quantum critical phenomena, exhibiting a rich variety of quantum phases [4,5,6,7,8]. Apart from a few cases [2], spin ladder systems are not exactly solvable, therefore it is necessary to develop various techniques, both analytical and numerical, to investigate their physical properties [9,10,11,12,13,14,15,16,17,18,19].
An efficient tensor network (TN) algorithm has been developed which is tailored to translationally invariant infinite-size spin ladder systems [20]. The algorithm is based on an adaptation to the ladder geometry of projected entangled pair states [21], a representation of quantum many-body wave functions. The spin ladder TN algorithm has been successfully applied to a variety of models [20,22], including the ferromagnetic frustrated two-leg ladder, the two-leg Heisenberg spin ladder with cyclic four-spin exchange and cross couplings, and the three-leg Heisenberg spin ladder with staggering dimerization. The algorithm is seen to be efficient compared to the density matrix renormalization group [12] and time evolving block decimation [23], at least as far as the memory cost is concerned.
The TN representation offers an efficient way to detect quantum phase transitions in many-body quantum systems [24,25] using the groundstate fidelity per lattice site [26,27,29,28,30,31] (for other work on the fidelity approach to quantum phase transitions, see also [32,33,34,35]). We recall that fidelity, as a measure of quantum state distinguishability in quantum information science, describes the distance between two given quantum states. It offers a powerful means to investigate critical phenomena in quantum many-body lattice systems, as demonstrated in Refs. [26,27,29,28,30,31]. For example, the groundstate fidelity per lattice site enables the characterization of drastic changes of the groundstate wave functions in the vicinity of a phase transition point. A systematic scheme to study critical phenomena in quantum many-body lattice systems has been outlined using this approach [36]. Once the groundstate phase diagram is mapped out by means of the groundstate fidelity per lattice site, the local order parameters can be derived from the reduced density matrices for a representative groundstate wave function in a symmetry-broken phase. Other phases without any longrange order can also be detected and further characterized by nontrival order parameters such as string order and pseudo-order parameters.
In this paper we use the TN algorithm developed in Ref. [20] to study a twoleg ladder model with anisotropic spin- 1 2 XXZ interactions along the legs and rungs of the ladder. A schematic phase diagram for this ladder model was mapped out previously using the sublattice entanglement of the groundstates [37]. A variant of this ladder model, with isotropic rung interactions, was shown to exhibit a rich phase diagram [9,38]. We perform an extensive numerical analysis of the fully anisotropic model to map out the groundstate phase diagram by evaluating the groundstate fidelity per lattice site from the groundstate wave functions. The phase transition points and thus the phase boundaries are detected through identifying pinch points on the groundstate fidelity surfaces, which arise from the distinct changes of the groundstate wave functions in the vicinity of the critical points as the model parameters are varied across a transition point. This requires a scan of the entire parameter space. The phase diagram unveiled in this way is shown in Figure 1. Varying the XXZ anisotropy parameter ∆ and the relative rung coupling strength J is seen to result in a rich diversity |0,0> |1,z> of phases -ferromagnetic (FM), stripe ferromagnetic (SF), rung singlet (RS), rung triplet (RT), Néel (N), stripe Néel (SN), Haldane (H) and two XY phases (XY 1 and XY 2). In particular, when |∆| > 1 the resulting Ising-type anisotropy breaks any continuous spin symmetry, in contrast to the XY -type regime for |∆| < 1. The XYtype phases without any long-range order are characterized in terms of string order parameters. In the symmetry-broken phases local order parameters are constructed from reduced density matrices for a representative groundstate wave function. In the XY -type regime varying the rung coupling J induces magnetic phases beyond standard ferromagnetism, for example, the rung-singlet and Haldane phases. For this model, we detect an additional rung-triplet phase as a result of the inherent competition between singlet formation and magnetic ordering.
The layout of the paper is as follows. The fully anisotropic spin-1 2 XXZ two-leg ladder model is described in Section 2, with the results using the TN algorithm presented in detail in Section 3. Concluding remarks, along with a discussion of previous results for this model, are given in Section 4.

The model
The fully anisotropic spin-1 2 XXZ two-leg ladder system can be considered as a pair of infinitely long spin chains coupled via rung interactions. All spin interactions are of the spin- 1 2 XXZ Heisenberg type. The Hamiltonian is defined by where Here S β α,i (β = x, y, z) is the spin-1 2 operator acting at site i on the α-th leg, J is the rung coupling between the two spins on a rung, and ∆ is the XXZ anisotropy parameter. Our aim is to map out the groundstate phase diagram of this model in the ∆-J plane.
In the FM phase, the ladder system orders ferromagnetically in both the leg and rung directions, whereas in the SF phase it orders ferromagnetically in the leg direction and antiferromagnetically in the rung direction. In contrast, in the N phase, the system orders antiferromagnetically in both the leg and rung directions, whereas in the SN phase it orders antiferromagnetically in the leg direction and ferromagnetically in the rung direction. Notice that in the FM, SF, N, and SN phases S x α,i = S y α,i = 0 and S z α,i = 0. The critical XY phases belong to the universality class of the Tomonaga-Luttinger liquid [39], which is known to exhibit a power-law decay of the spin-spin correlations with gapless excitations. The possibility for the existence of two different XY phases in a ladder system was pointed out in a slightly different model [17]. According to the Marshall-Lieb-Mattis theorem [40], the XY 1 phase is in the J < 0 region and the XY 2 phase is in the J > 0 region of the parameter space. Long-range order does not exist in the XY phases in the thermodynamic limit. Nevertheless, there is a practical way to characterize the XY phases in the context of the TN representation of the groundstate wave functions, which is achieved by defining a pseudo-order parameter arising from the finiteness of the bond dimension χ [31]. This offers a convenient means to numerically determine the phase boundaries between the XY phase and other phases. For the XY 1 and XY 2 phases, S z α,i = 0 and S x α,i 2 + S y α,i 2 = 0. The two-spin states on a rung are written as Here the subscripts 1 and 2 label the different spins on the same rung (see, e.g., Ref. [41]). The first state is the singlet and the remaining three states constitute the triplet. As can readily be seen from the Hamiltonian (1), the groundstate of the system in the strong coupling limit J → ±∞ corresponds to the RS (|0, 0 ) and RT (|1, z ) phases, respectively. It should be noted that the RS phase, the RT phase and the H phase lack long-range order in the conventional Landau-Ginzburg-Wilson sense. However, each are characterized by a suitably modified string order parameter [42]. For the RS, RT and H phases S β α,i = 0.

Results
To map out the groundstate phase diagram, we first consider the regions J > 0 and J < 0, with the anisotropy ∆ as a control parameter. We then fix ∆ and vary the rung coupling J. Illustrative examples are given in detail below for the values J = ±1. In our approach, apart from order parameters, which we also make use of, the groundstate fidelity per lattice site is used as an important tool to detect quantum phase transitions. For two different groundstates |ψ(x) and |ψ(x ′ ) in a quantum system corresponding to different values x and x ′ of a control parameter, the fidelity F (x, x ′ ) = | ψ(x)|ψ(x ′ ) | is defined as a measure of the overlap between the two states. For large but finite size N, the fidelity F scales as d N , with d interpreted as the ground state fidelity per site d = lim N →∞ F 1/N N . As a well-defined property even in the thermodynamic limit, the groundstate fidelity per lattice site d enjoys some properties inherited from the fidelity F , namely: In previous studies the fidelity per lattice site d was demonstrated to be a useful detector for different types of quantum phase transition via its singularity structure [27,28,26,29,30,31]. In this approach, a pinch point in the d(x, x ′ ) surface occurs at a quantum phase transition, i.e., at each phase transition point x c , a pinch point (x c , x c ) occurs on the surface of fidelity per lattice site d(x, x ′ ) at the crossing of two singular lines x = x c and x ′ = x c . In these studies the groundstate wavefunctions are determined using TN algorithms [21,43]. The combined groundstate fidelity and TN approach has been tailored to include ladder systems [20].
A key point in the investigation of the two-leg quantum spin ladders via the TN approach is that the computational cost scales as χ 6 , where χ is the bond dimension for each of the four necessary four-index tensors. These tensors need to be updated simultaneously, with the updating procedure closely related to the infinite PEPS [44] and translationally invariant MPS [45] algorithms. A summary of the TN approach for quantum spin ladders is given in Appendix A. For the spin-1 2 isotropic XXX two-leg ladder it was demonstrated earlier that the TN approach could accurately determine the groundstate wave function, with bond dimension χ = 6 outperforming the DMRG results [20]. Of necessity for the calculation of the order parameters, the reduced density matrix ρ can be computed directly from the TN representation of the groundstate wave functions [20]. The reduced density matrix displays different nonzero-entry structures in different phases.
Fixing the rung coupling to the value J = 1 and varying the anisotropy parameter ∆ from −2 to 2, we are able to determine the phase boundaries between the FM phase, the XY 2 phase, the RS phase and the Néel phase. The groundstate fidelity per site, d(∆ 1 , ∆ 2 ), is shown in Figure 2 as a function of ∆ 1 and ∆ 2 , calculated with bond dimension χ = 6. Here the range of the plots is restricted to highlight the fidelity surface and, in particular, the pinch points. The results in Figure 2(a) and Figure 2(c) differ very little for a larger bond dimension χ, indicating that the computational results are almost saturated for bond dimension χ = 6. However, the pinch point in Figure 2(b) shifts noticeably when χ increases from χ = 6 to χ = 12, as Figure 3 shows. Three pinch points are identified at (∆ c1 , ∆ c1 ), (∆ c2 , ∆ c2 ) and (∆ c3 , ∆ c3 ) on the fidelity surfaces, indicating that there are three phase transition points. These are located at (i) ∆ c1 = −1.00, (ii) ∆ c2 = 0.00 and (iii) ∆ c3 = 1.43. In this way we are able to identify four different phases, separated by three transition points. It remains to characterize these phases in terms of local or nonlocal order parameters.
In the FM phase, the one-site reduced density matrix yields the local order parameter, Our numerical results indicate that the ladder system is in the FM phase for ∆ < −1.00 (see Figure 3(a)).
In the Néel phase, we similarly consider the local order parameter In the Néel phase, the Néel order O N exhibits a long-range order. The evaluation of the order parameter indicates that the ladder system is in the Néel phase for ∆ > 1.43 (see Figure 3(a)).
In the RS phase, the groundstate wave function may be approximated by the product of local rung singlets. As usual, two distinct string order parameters, O odd and O even , may be used to characterize the RS phase. For two-leg ladder models these are defined by . The odd and even string order parameters are actually mutually exclusive in the RS, RT and H phases. Our results for O odd and O even indicate that the RS phase persists in the region ∆ c2 < ∆ < ∆ c3 with J = 1 (see Figure 3(a)).
In the XY phases, the TN algorithm automatically leads to infinite degenerate groundstates, arising from pseudo spontaneous symmetry breaking of the continuous U(1) symmetry, due to the finiteness of the bond dimension χ. This allows the introduction of a pseudo-order parameter that must scale to zero, in order to keep consistency with the Mermin-Wagner theorem. As suggested in Refs. [31] and [41], the pseudo-order parameters in the XY 1 and XY 2 phases may be defined as O even (12) O odd (12) N RS XY2 FM O N (6) O N (12) O FM (6) O FM (12) O 2 (4) O 2 (6) O 2 (12) O Figure 3. (a) The local order parameter O FM (χ) in the FM phase, the pseudo-order parameter O 2 (χ) in the XY 2 phase, the local order parameter O N (χ) in the Néel phase and the string order parameters O odd (χ) and O even (χ) in the RS phase as a function of the anisotropy ∆, with J = 1. Three phase transition points, ∆ c1 , ∆ c2 and ∆ c3 , are distinguished by the behavior of the order parameters. The phase transition points ∆ c1 = −1.00 and ∆ c3 = 1.43 are saturated for bond dimension χ = 6, whereas ∆ c2 shifts with increasing χ. (b) The scaling of the pseudo-order parameter O 2 (χ) for ∆ = −0.5 (see text). (c) The "pseudo-critical" point estimates ∆ c2 (χ) obtained from the pseudo-order parameter O 2 (χ) and their fitting function (see text).   (7))/E(7) J=1 Figure 5. Relative error ε in the estimates of the groundstate energy as a function of ∆ at J = 1 for increasing bond dimension χ. The relative errors are defined by is the groundstate energy estimate for given bond dimension χ n . and The pseudo-order parameter O 2 in the XY 2 phase is plotted as a function of the anisotropy ∆ for different values of the bond dimension χ in Figure 3(a). This order parameter O 2 is also plotted as a function of χ for ∆ = −0.5 in Figure 3(b). It is clearly seen that O 2 decreases as χ increases. Here the pseudo-order parameter scales to zero according to (1) and c 1 = 0.26 (1). Such scaling was adopted in previous studies of pseudo-order parameters [31,41] and ensures that they vanish for infinite bond dimension as expected.
The estimates of the "pseudo-critical" values ∆ c2 (χ) at which the pseudo-order parameter O 2 (χ) becomes zero are shown in Figure 3(c). These points are extrapolated with respect to the bond dimension χ, with an extrapolation function ∆(χ) = a 2 + b 2 χ −c 2 . The results imply the fitting constants a 2 = −0.110(2), b 2 = 0.585(7) and c 2 = 0.93(2). In the limit χ → ∞ we have ∆(∞) = a 2 , which is the estimate for the critical point ∆ c2 between the XY 2 and RS phases. The term "pseudo-critical" point is used because such estimates are obtained using the pseudo-order parameters. Nevertheless they yield estimates for real critical points.
The existence of the FM phase, the XY 2 phase, the Néel phase and the RS phase is numerically confirmed in this way, as seen from Figure 3. Specifically, the three phase transition points occur at the values ∆ c1 = −1.00, ∆ c2 ≃ −0.11 and ∆ c3 = 1.43. If the anisotropy coupling ∆ is tuned as a control parameter, the ladder system undergoes a first-order phase transition at ∆ c1 , with continuous phase transitions at ∆ c2 and ∆ c3 . The different nature of the first-order transition at ∆ c1 can be seen clearly in the fidelity surface in Figure 2(a), compared to the continuous transitions in Figure 2(b) and Figure 2(c). In general, for continuous phase transitions, the fidelity surface shows continuous behaviour near critical points, while discontinuous phase transitions show discontinuous behaviour.
In the limit J → ∞ the XY 2 phase vanishes. In addition for general J > 0 the line ∆ = −1 is the phase boundary between the FM phase ∆ < −1 and the XY 2 phase ∆ > −1. We have also confirmed that ∆ ≈ J is the phase boundary between the RS phase and the Néel phase, when ∆ → ∞ and J → ∞. There are thus three phases when J → ∞: the FM phase, the RS phase and the Néel phase.
Concerning the accuracy of our results, estimates of the groundstate energy per site with increasing bond dimension are shown in Figure 4. The estimates are seen to converge rapidly, as quantified in Figure 5. We find that almost all of the relative errors reach to order 10 −5 . This indicates the reliability of the algorithm for generating accurate groundstate wave functions for spin ladders, given a preset error tolerance. The corresponding convergence of the order parameter estimates is shown in Figure 6. As is well understood, states in gapless phases, such as the XY 1 and XY 2 phases, are more difficult to compute than states in phases with a large gap, such as the FM, N and SF phases, which converge much faster, yielding the same accuracy with a smaller bond dimension. Phases with a small gap, such as the H phase, are also more difficult to compute. In general this slower convergence in gapless (critical) versus gapped phases is a feature of numerical entanglement based approaches, such as TN methods [43]. This is because the overall strategy of the TN method is more efficient for gapped systems. By design, the finite bond dimension already induces a gap in the system. From the perspective of entanglement in quantum critical phenomena [46], gapped systems obey an area law, which can be well described using relatively small bond dimensions. On the other hand, for gapless systems, the entanglement diverges as the gap vanishes, requiring larger, in principle infinite, bond dimensions. Alternatively, for the same size bond dimension, more iteration steps are necessary for convergence in gapless compared to gapped systems. Our algorithm is seen to work sufficiently well in the phases under consideration. However, the critical XY 1 and XY 2 phases require additional extrapolation to infinite-size bond dimension.

J
In similar fashion, we determine the phase boundaries between the SF, RT (|1, z ), XY 1, H and SN phases. The groundstate fidelity per lattice site d(∆ 1 , ∆ 2 ) is shown am of the fully anisotropic two-leg XXZ n ladder 12    In the SN phase, we are similarly able to consider the local order parameter In this phase the ladder system exhibits a long-range order for ∆ 06 (see Figure 8(a)).
In the H phase, the even string order parameter even vanishes, but the odd string order parameter is nonzero. Conversely, in the RT ( , z ) phase, odd = 0 and the even string order parameter is nonzero (see Figure 8(a)).
In the XY 1 phase, the pseudo-order parameter as a function of the anisotropy ∆ for different bond dimension is shown in Figure 8(a). The pseudo-order parameter is plotted as a function of for ∆ = 0 5 in Figure 8(b). Here we have performed an extrapolation with respect to , with the fitting function ) = (1 + ), where = 0 467(7), = 0 006(4) and = 1 04(5). This is again consistent with previous studies of pseudo-order parameters. The estimates of the "pseudo-critical" points ∆ ) and ∆ ) on either side of the XY 1 phase are shown in Figure 8(c) and Figure 8  as a function of ∆ 1 and ∆ 2 in Figure 7 with bond dimension χ = 6. In this region four pinch points, (∆ c4 , ∆ c4 ), (∆ c5 , ∆ c5 ), (∆ c6 , ∆ c6 ) and (∆ c7 , ∆ c7 ), are identified on the fidelity surfaces, corresponding to the continuous phase transition points ∆ c4 = −1.26, ∆ c5 = 0.00, ∆ c6 = 0.84 and ∆ c7 = 1.06. These points separate five different phases, which we now proceed to characterize in terms of the order parameters.
In the SF phase, the one-site reduced density matrix yields the local order parameter Our results show that the ladder system exhibits a long-range order in the SF phase for ∆ < −1.26 (see Figure 8(a)).
In the SN phase, we are similarly able to consider the local order parameter In this phase the ladder system exhibits a long-range order for ∆ > 1.06 (see Figure 8(a)).
In the H phase, the even string order parameter O even vanishes, but the odd string order parameter is nonzero. Conversely, in the RT (|1, z ) phase, O odd = 0 and the even string order parameter is nonzero (see Figure 8(a)).
It follows that we have been able to compute the order parameters O SF and O SN in the SF and SN phases, the pseudo-order parameter O 1 in the XY 1 phase, the string order parameter O even in the RT (|1, z ) phase, and the string order parameter O odd in the H phase. The results are shown in Figure 8. With the anisotropy ∆ as control parameter, the ladder system undergoes four continuous phase transitions at the values ∆ c4 , ∆ c5 , ∆ c6 and ∆ c7 . In the limit J → −∞ the XY 1 phase and the H phase vanish. We note that the phase boundary between the SF phase and the RT (|1, z ) phase is located roughly near the line ∆ ≈ J in the limits ∆ → −∞ and J → −∞, while the phase boundary between the RT (|1, z ) and SN phases is located roughly near the line ∆ ∼ 1 as J → −∞. As a result there are three phases for the ladder system in the limit J → −∞, i.e., the SF, RT (|1, z ) and SN phases.
The accuracy of the results is similar to that for J > 0, with estimates of the groundstate energy per site with increasing bond dimension shown in Figure 9. The estimates also converge rapidly, as quantified in Figure 10. The corresponding (E(7)-E(6))/E (6) (E(12)-E (7))/E (7) J=-1 Figure 10. Relative error ε in the estimates of the groundstate energy as a function of ∆ at J = −1 for increasing bond dimension χ. The relative errors are defined by ε = (E(χ n+1 ) − E(χ n ))/E(χ n ) where E(χ n ) is the groundstate energy estimate for given bond dimension χ n .  convergence of the order parameter estimates is shown in Figure 11 We now turn to look more closely at the XY phase, which as we have seen, splits into the two different XY 1 and XY 2 phases identified by the local pseudo-order parameters O 1 and O 2 obtained from the groundstate fidelity per lattice site. As we shall see, for −1 ≤ ∆ ≤ 1, the phase boundary between the XY 1 and XY 2 phases is the line J = 0. Here we fix the anisotropy coupling ∆ = 0 and take J as control parameter in the range −0.5 ≤ J ≤ 0.5 to sweep across this boundary. Figure 12 shows a plot of the groundstate fidelity per site d(J 1 , J 2 ) with bond dimension χ = 6. In this case there is a clear indication of a continuous phase transition point at J = 0.00, as a pinch point occurs at this value with bond dimension χ = 6. Figure 13 shows a plot of the pseudo-order parameter O 1 in the XY 1 phase and the pseudo-order parameter O 2 in the XY 2 phase as a function of J. As already remarked, a feature of the pseudo-order parameters is that they vanish as the bond dimension χ increases.

Conclusion
We have systematically investigated the groundstate phase diagram of the infinite fully anisotropic spin-1 2 XXZ two-leg ladder model (1). This has been achieved by exploiting a TN algorithm tailored to the geometry of quantum spin ladder systems [20], allowing efficient evaluation of the groundstate fidelity per lattice site from the TN representation of the groundstate wave functions for the infinite-size quantum spin ladder. Illustrative results have been presented in Sec. III along the lines J = ±1 with varying ∆ and the line ∆ = 0 with varying J. Our computational results using the groundstate fidelity per lattice site are consistent with results obtained from the construction of local, and where appropriate, nonlocal, order parameters in the different phases.
The full phase diagram, as obtained by this approach, is shown in Figure 1. This phase diagram, featuring nine distinct quantum phases, points to the rich diversity of physics in the fully anisotropic spin-1 2 XXZ two-leg ladder model. The precise nature of the phase diagram is a significant extension on the previously obtained schematic phase diagram for this model, which was mapped out by identifying ridges and valleys in contour plots of the rescaled block entanglement entropy corresponding to possible phase boundaries [37]. With that approach, which was aimed at demonstrating the usefulness of sublattice entanglement entropy as a means of exploring quantum phase transitions across a range of different spin systems, the nature and number of the different phases of the ladder model were not discussed in detail. The phase diagram of the fully anisotropic spin-1 2 XXZ two-leg ladder model also differs in significant aspects compared to the phase diagram of the spin-1 2 XXZ two-leg ladder model with isotropic rung interactions, the precise nature of which has been the subject of debate [9,38]. With isotropic rung interactions, the model also exhibits FM, SF, RS, N, SN, H, XY 1 and XY 2 phases. However, for the fully anisotropic model considered here, there is an additional rung-triplet (RT(|1, z )) phase, with the XY 1 and XY 2 phases extending over the whole −1 ≤ ∆ ≤ 1 region in the vicinity of J ≈ 0. Moreover, the extent of the Haldane (H) phase is seen to be significantly diminished for fully anisotropic interactions. In general the difference between both the number of quantum phases and the phase diagrams of the two models highlights the key role played by anisotropies in the competition between singlet formation and magnetic ordering in quantum spin systems.

Acknowledgments
We thank Sam Young Cho, Bing-Quan Hu and Fabian Essler for helpful discussions at various stages of this work. Appendix A. Tensor network representation for the spin-1 2 two-leg ladder In this Appendix we summarise the relevant details of the TN representation for spin ladders [20]. The hamiltonian of the two-leg ladder model (1)  It is sufficient to assume that the TN representation for the wave function is translational invariant under shifts by two lattice sites along the legs. The groundstate wave function for an infinite length ladder can be described by the four different four-index tensors A s ℓrd , B s ℓrd , C s ℓru and D s ℓru . The tensor A s ℓrd is depicted in Figure A1(i), with s = 1, 2 for spin-1 2 and with ℓ, r, u, d = 1, . . . , χ, the four inner bond indices, where χ is the bond dimension. A TN representation for the groundstate wave function is shown in Figure A1(ii) with one of the two equivalent choices shown for the plaquette (unit cell).
A gradient-directed random walk method is applied to compute the groundstates. The double tensors alrd, blrd, clrũ and dlrũ are introduced, withl = (ℓ, ℓ ′ ),r = (r, r ′ ), u = (u, u ′ ) andd = (d, d ′ ). They are formed from the tensors A s ℓrd , B s ℓrd , C s ℓru and D s ℓru , and their complex conjugates, as depicted in Figure A1(i) for the double tensor a. With these double tensors, the TN for the norm of a wave function is shown in Figure A1(iii). The two different but equivalent choices for the transfer matrix T are made from the plaquettes abdc and bacd. In the above, the notationl = (ℓ, ℓ ′ ) corresponds to the two independent indices ℓ and ℓ ′ , withl = 1, 2, . . . , χ 2 . Similarly forr,ũ andd.
For a randomly chosen initial state |ψ 0 , the energy per plaquette e is The groundstate energy per plaquette also admits a TN representation as shown in Figure A1(v) and (vi) for plaquette ABDC, which absorbs the operator h acting on a plaquette for an infinite-size spin ladder. For each choice of plaquette, the dominant left and right eigenvectors of the transfer matrix T constitute the environment tensors, visualised in Figure A1(iv) and (vi) for plaquette ABDC. In this way the energy per plaquette e ABDC is obtained. The same procedure may be used to compute the energy per plaquette e BACD for an operator acting on the ith plaquette BACD.
To update the TN representation, the energy gradient is computed with respect to Four-index tensor lrd to represent a tensor network (TN) representation for the ground-state wave function for an infinite-size two-leg spin ladder, with a physical index, , and a double tensor x tensor rd x conjugate ( , with ℓ, ℓ ),˜), and ). (ii) Pictorial representation for a TN state leg and rung bonds, which are used to absorb an operator acting on the TN representation for the norm an infinite-size norm TN, which is constructed from four double tensors , and , with ℓ, ℓ ), ),˜), and on sites 1) along the legs. Here, 1 for two-leg ladder over all the possible plaquettes by taking ∈ {−∞ · · · ∞}. Assume that the TN representation for a wave ys the translational invariance under shifts by two lattice sites along the legs. Ground state wave functions for an two-leg ladder system can be described by four di x tensors rd rd ru , and ru as shown in a four-index tensor rd , with physical index and , ..., @ of the local Hilbert space @ 2 for spin 1 2 systems, and four inner bond indices , ..., , with A TN representation for the ground-state wave function is shown for an infinite-size two-leg spin ladder As wn in.?? two equivalent choices of the unit cell , and , and , and A gradient-directed random walk method is applied to commutate ground states for an infinite-size spin two-leg ladder. To be convenience, double tensors , and ℓ, ℓ ),˜),˜), and ). y form from the four-index tensors rd rd ru , and ru , and their complex conjugates, see  Fig. (iv) and (vi) . Then, the energy per unit cell ABDC is obtained. The same procedure may be used to compute the gy per unit cell BACD an operator acting on the , and Figure A1.
(i) Four-index tensor A s lrd used in the TN representation for the groundstate wave function for the infinite-size two-leg spin ladder. Also shown is the double tensor alrd formed from the tensor A s ℓrd and its complex conjugate (A * ) s ℓ ′ r ′ d ′ , withl = (ℓ, ℓ ′ ),r = (r, r ′ ) andd = (d, d ′ ). (ii) Pictorial representation for a TN state |ψ with leg and rung bonds, which are used to absorb an operator acting on the ith plaquette. (iii) TN representation for the norm of a groundstate wavefunction in the infinite ladder. Also shown is the transfer matrix T , which is constructed from the four double tensors alrd, blrd, clrũ and dlrũ, withl,r andd defined as above and u = (u, u ′ ). The contributions to the energy gradient come from three parts (for a given choice of plaquette on which the hamiltonian acts), as shown in Figure A2. To obtain the energy gradient it is necessary to sum the contributions from both plaquettes abdc and bacd.
The tensor A s ℓrd is then updated according to where δ denotes the update step size. In this procedure the tensors A s ℓrd , B s ℓrd , C s ℓru and D s ℓru are updated simultaneously. Repeating this procedure until the groundstate energy per plaquette converges leads to the groundstate wave function as approximated in the TN representation. Figure A2. The contribution to the energy gradient (A.2) for the infinite-size spin ladder consists of three parts: (i), (ii) and (iii) depending on the location of the hole cell with tensor A s ℓrd removed. In the latter two cases, there are m cells between the hole and hamiltonian cell, where m ∈ (0, 1, 2, 3, . . .). The hole cell is visualised in (iv), with the tensor A s ℓrd removed. Note that only the contribution from the plaquette abdc, as labeled here, is shown. There is also a similar contribution from the plaquette bacd.