Non-Hermitian effects of the intrinsic signs in topologically ordered wavefunctions

Negative signs in many-body wavefunctions play an important role in quantum mechanics because interference relies on cancellation between amplitudes of opposite signs. The ground-state wavefunction of double semion model contains negative signs that cannot be removed by any local transformation. Here we study the quantum effects of these intrinsic negative signs. By proposing a generic double semion wavefunction in tensor network representation, we show that its norm can be mapped to the partition function of a triangular lattice Ashkin-Teller model with imaginary interactions. We use numerical tensor-network methods to solve this non-Hermitian model with parity-time symmetry and determine a global phase diagram. In particular, we find a dense loop phase described by non-unitary conformal field theory and a parity-time-symmetry breaking phase characterized by the zeros of the partition function. Therefore, our work establishes a connection between the intrinsic signs in the topological wavefunction and non-unitary phases in the parity-time-symmetric non-Hermitian statistical model. While negative signs in wave functions are known to play a relevant role in many-body quantum mechanics, their physical consequences on the nature of topologically ordered phases of matter is still unclear. Here, the authors use tensor network methods to show that the intrinsic negative signs of the double semion model are responsible for the non-unitary phases and non-unitary phase transitions found in the phase diagram.

I n quantum mechanics, negative signs in the many-body wavefunctions play an important role, as the interference relies on a cancellation between amplitudes of opposite signs of the wavefunctions. Recently, topologically ordered phases of matter with anyonic excitations have attracted considerable attention due to the potential application of topological quantum computation [1][2][3][4][5] . A prototype topologically ordered state is given by the toric code (TC) model 1 , and its ground-state wavefunction is an equal weight superposition of all closed domain-wall loops. A closely related topologically ordered state is the double semion (DS) model, whose ground-state wavefunction is a superposition of all closed domainwall loops weighted by a π phase multiplying the number of closed loops 3,6 . Due to the presence of such non-local phase factors, lowenergy excitations of the DS phase contain semions with opposite chirality and bosonic bound states of two semions, in contrast to the anyonic excitations (charge-vertices, magnetic-fluxes, and their fermionic bound states) of the TC model. Actually it has been speculated that the negative signs in the DS wavefunction are intrinsic and cannot be removed by any local transformation 7,8 , making it the first prototype example of a topological ground state wavefunction with intrinsic signs. However, the norm of the DS wavefunction is identical to that of TC 9,10 . Therefore, important open questions arise what physical consequences are caused by these intrinsic signs and whether unknown anyon condensed phases are resulted in adjacent to the DS phase.
Meanwhile, it has been shown that non-Hermitian and paritytime (PT) symmetric quantum systems provide a useful description of dissipative quantum systems 11,12 . The fundamental property of these systems is that an unusual spontaneous PTsymmetry breaking (PTSB) phase with spectral singularity occurs. In correlated many-body PT-symmetric systems, an exotic universality class of quantum criticality occurs due to the presence of negative/complex Boltzmann weights 13 . In statistical mechanics, a well-known example is the Yang-Lee edge singularity 14,15 , where an imaginary magnetic field in the hightemperature Ising model was demonstrated to trigger an exotic phase transition described by non-unitary conformal field theory (CFT).
In this work, we attempt to establish a connection between the intrinsic sign problem in the topological wavefunction and the PTsymmetric statistical model with negative Boltzmann weights, and provide insight into the quantum effects of the intrinsic signs. By proposing a generic DS wavefunction with two distinct tuning parameters for a two-dimensional hexagonal lattice model, we explore the possible topological phase transitions out of the DS phase. In the tensor network representation, we express the non-local phase factors in the DS wavefunction in terms of auxiliary spins on the dual triangular lattice. By integrating out the physical degrees of freedom, we map this wavefunction norm to the partition function of a twodimensional triangular lattice Ashkin-Teller model with imaginary magnetic fields and imaginary three-spin triangular face interactions. To solve this PT-symmetric statistical model with negative Boltzmann weights, we employ the corner-transfer-matrix renormalization group (CTMRG) tensor-network method [16][17][18][19] . A global phase diagram is fully determined. Adjacent to the DS phase, we find a gapped dilute loop phase with condensed bosonic anyons, a gapless dense loop phase described by a non-unitary CFT, and a PTSB phase characterized by the zeros of the partition function in the thermodynamic limit. The PTSB transitions described by exceptional points 20 are continuous, which have no Hermitian critical counterparts. In this sense, the derived phase diagram for the DS phase is distinctly different from that for the TC phase. Further, the properties of the emerged PTSB phase cannot be reproduced by the partition function with positive Boltzmann weights. Although no dissipation is involved, the intrinsic sign problem in the topological wavefunction is thus closely connected with non-unitary phases in a non-Hermitian PT-symmetric statistical model.

Results
The DS tensor-network wavefunction. The DS model is defined by the following Hamiltonian which acts on quantum spin-1/2 operators living on the edges of a two-dimensional hexagonal lattice 3 Here the vertex term A v and the plaquette term B p shown in Fig. 1a are given by where E(v) is the set of edges around a vertex v, E(p) is the set of inner edges and e EðpÞ is the set of outer edges around a hexagon p. Operator A v = −1 is associated to a semion or anti-semion excitation, while operator B p = −1 is a bosonic excitation of a semion-antisemion pair. When all vertex terms A v = +1 restrict the Hilbert space of states to the zero-flux subspace, the ground state of the Hamiltonian is stabilized by B p = 1 for all plaquettes: where " j i is the eigenvector of σ z with eigenvalue σ = +1, and N is the total number of spins. Actually, the loop representation provides a very intuitive way to understand this wavefunction, in which σ = −1 and σ = +1 states are interpreted as the presence a The vertex term A v (green) and plaquette term B p (yellow) of the DS model, where the black dots are the locations of the physical degrees of freedom. We have used the Pauli matrices σ x , σ z , and defined the operator g ¼ i ð1þσ z Þ=2 . b A typical loop configuration consists of two closed loops (red lines), where the physical spin-down states σ = −1 form closed loops. The red/black circles represent the physical spin-down/up state, and σ is the quantum number of σ z . c The auxiliary spin (blue arrows) configuration corresponds to the loop (red lines) configuration, where the loops are regarded as the domain-walls of auxiliary spins. d The auxiliary spins (blue arrows) s p , s q are defined on the sites p, q and the physical spin (black/red circles) σ p,q is defined between the sites p and q, which satisfy s p s q = σ p,q . The domain-walls (red lines) are assigned orientations according to the rule that down/up auxiliary spins are on the left/right. A left-turn/right-turn of the domain wall is associated with a phase factor e ±iα , where the parameter α is fixed as e i6α = −1.
or absence of loops on the edges, as shown in Fig. 1b. Expanding this product, we have the ground state as a superposition of all closed-loop configurations with alternative signs 3,6 where #L is the number of contractible closed loops in the configuration L. Moreover, it is very useful to express the DS wavefunction in terms of tensor-network state (TNS) 21 . The loops can be regarded as the domain-walls of local auxiliary Ising spins, which form a dual triangular lattice, as shown in Fig. 1c. These auxiliary Ising spins contain the crucial information about anyonic excitations 22,23 . Importantly, when the domain-walls of the Ising paramagnet are assigned orientations according to the rule displayed in Fig. 1d, the non-local signs in the wavefunction can be locally expressed in terms of these auxiliary Ising spins. In the hexagonal lattice, the difference between the left/right-turnings for a closed contractible orientated loop must be six, so we have e i6α = −1. Then the many-body entangled state between the physical spins and auxiliary spins is given bỹ where 〈pq〉 denotes the nearest-neighbor dual lattice sites, 〈pqr〉 stands for the minimal triangular faces. The DS tensor-network wavefunction is obtained by summing over all the auxiliary spins: Ψ 0 j i ¼ P fsg sjΨ . The above DS wavefunction is just the fixed-point wavefunction with zero correlation length. To study the phase transitions out of the DS phase, we need a generically deformed DS wavefunction with tuning parameters, and we can employ the wavefunction approach to reveal its essential physics. When subjected to two different magnetic fields h 0 x and h 0 z , the DS model with the additional Zeeman terms are no longer exactly solvable. If these additional terms are treated perturbably, the ground-state wavefunction is obtained When the wavefunction corrections are expressed as the operator product with the modified magnetic field parameters, we have which can be regarded as a generic DS wavefunction in an expanded parameter space. Actually, the similar deformation has been used to express the generic TC wavefunction [23][24][25] and Fibonacci quantum-net wavefunction 26  and can be regarded as a partition function of a statistical model: where s p and t p are the auxiliary Ising spins in the ket and the bra layers, the 2-spin and 4-spin couplings are given by and the parameter α has been fixed as e i6α = −1. This is a twodimensional triangular lattice Ashkin-Teller model with imaginary magnetic fields and imaginary three-spin triangular face interactions. The imaginary terms of Eq. (9) originated from the negative signs in the DS wavefunction, and this statistical model is PT-symmetric, i.e., invariant under the combined operation of parity symmetry P (s p → −s p and t p → −t p ) and time-reversal symmetry T (i → −i). This PT-symmetry ensures that all eigenvalues of the transfer matrix operator of the partition function are either real or complex conjugate pairs 29 .
Since this model has the symmetry h x → −h x , the tuning parameters are limited to 0 ≤ h ≤ 1 and −π/2 ≤ θ ≤ π/2 in the following. In general, this model cannot be solved analytically except the following two special limits: 1. When h x = 0, we have J 4 → ∞, the Ising spins in two layers are locked together, and the imaginary terms vanish, leading to an inter-layer partial order 〈st〉 = ±1. The norm of the DS and TC wavefunctions become identical 9,10 and Eq. (9) is reduced to a two-dimensional triangular lattice Ising model: For h z > 0, it is known that a ferromagnetic critical point 30 exists at h z = (3 1/4 −1)/(3 1/4 + 1) and denoted by H. This critical point separates the intra-layer ferromagnetic phase and intra-layer paramagnetic phase, corresponding to the dilute loop phase and the DS phase, respectively. For the antiferromagnetic coupling h z < 0, however, there is no phase transition up to the multicritical point (h z = −1) due to the presence of spin frustration 31 .
2. When h = 1 and −π/2 < θ < π/2, the inter-layer coupling J 4 = 0 reduces the model Eq. (9) to two decoupled single-layer Ising spin models with J ¼ 1 4 log 1Àsin θ 1þsin θ . In the loop representation, the partition function of H s is just the single-layer O(n = −1) loop model: where l L is the combined length of the closed contractible loops. By mapping Eq. (11) to the antiferromagnetic Potts model, two exactly solvable critical points (A and E) were found [32][33][34] ) is characterized by a non-unitary CFT with a central charge c = −3/5 or an effective central charge ) corresponds to the fixed point of the dense loop phase 35 , which is described by a non-unitary CFT with a central charge c = −7 or an effective central charge c * = 1. However, when the two-spin coupling J is very small, the imaginary terms in Eq. (11) become dominate, yielding zeros of the partition function. Namely, the weights for different configurations in the partition function exactly cancel and the free energy is no longer continuous. In the following, we will employ the tensor-network numerical methods and find another phase transition point L at e 2J ≈ 1.18. Further, the zeros of the partition function are actually distributed in a region ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 À ffiffi ffi 3 p p ≤ e 2J <1:18, forming a non-unitary spontaneously PTSB phase.
Tensor-network numerical calculations. To establish the global phase diagram of the generic DS wavefunction, we will apply the numerical CTMRG method [16][17][18][19] to the PT-symmetric statistical model Eq. (9) with negative Boltzmann weights. In order to obtain the local tensors of TNS, each auxiliary spin s p is represented by a circle in Fig. 2a, and then the generic wavefunction is expressed as a double-line TNS, whose local tensors are defined on the vertices of the hexagonal lattice. For convenience, we transform the hexagonal lattice into the square lattice by combining two local tensors with nearest-neighbor triangle faces into one local tensor Q. Contracting the physical indices of Q and Q * leads to the local double-layer tensor Q shown in Fig. 2b, and the wavefunction norm is expressed by where "tTr" denotes the tensor contraction over all auxiliary indices, T is the column-to-column transfer matrix operator displayed in Fig. 2c, and L x is the number of columns. To determine the phase diagram, we need to calculate the dominant eigenvalues of transfer operator T, which contains the crucial information about the phase transitions. Since the spectrum can be complex, we sort the eigenvalues according to their moduli. As a comparison, we first calculate the phase diagram for the generic TC wavefunction. Since the imaginary terms are absent, Eq. (9) reduces to the two-dimensional Ashkin-Teller model on a triangular lattice, and the transfer operator is Hermitian. So the powerful variational uniform matrix-product-state (VUMPS) method 19,36,37 can be applied. In Fig. 3a, we show the correlation length ξ as a function of h z along the h x = 0.1 cut, and a sharp peak appears at h z = 0.1378, indicating a continuous phase transition from the inter-layer partial ordered phase (〈s〉 = 〈t〉 = 0, 〈st〉 ≠ 0) to the ferromagnetic phase (〈s〉 ≠ 0, 〈t〉 ≠ 0). In Fig. 3b, ξ along the h z = −0.1 cut is displayed, and a sharp peak around h x = 0.3196 represents a continuous phase transition from the partial ordered phase to the paramagnetic phase (〈s〉 = 〈t〉 = 〈st〉 = 0). Moreover, in Fig. 3c, a less sharper peak appears at h z = 0.4590 along the h z = 0.7−h x cut, implying the phase transition from the ferromagnetic phase to the paramagnetic phase.
It has been proved that 23 the partial ordered phase, ferromagnetic phase, and paramagnetic phase of the Ashkin-Teller model correspond to the TC phase, Higgs phase, and confining phase, respectively. So a global phase diagram can be derived as displayed The local double-line tensor Q on the square lattice is obtained by contracting two local tensors on the hexagonal lattice. The local doublelayer tensor Q consists of two layers, where s and t denote the auxiliary spins in ket and bra layers, and σ is the physical degrees of freedom. c The one-dimensional column-to-column transfer matrix operator T of the double-layer tensor network.  in Fig. 4a. The TC phase is enclosed by the Higgs and confining phases, corresponding to the electric and magnetic anyon condensed phases, respectively. Both condensed phases are gapped, and there is a critical line IJ between them, characterized by the Kosterlitz-Thouless transition. Meanwhile, the critical lines JH and JKF belong to the two-dimensional Ising universality class. The points I and J can be exactly determined as ðh x ; h z Þ ¼ ð ffiffi ffi 3 p =2; 1=2Þ and (h x , h z ) = (2/5, 1/5), because the first corresponds to the critical point of two-decoupled Ising models (J 4 = 0) and the second to the critical point of a four-state Potts model (J = J 4 ).
To establish the corresponding phase diagram of the generic DS wavefunction displayed in Fig. 4b, we have to employ the CTMRG method because the transfer matrix operator is non-Hermitian. The correlation length and the dominant eigenvalues are calculated and shown in Fig. 5a-c to determine the phase boundaries. The associated critical properties are deduced by estimating the central charges (Fig. 5d, e). In Fig. 5a, the correlation length is shown as a function of h z along the cut h x = 0.5. As h z decreasing from unity, a sharp peak first appears at h z ≈ 0.271, indicating a continuous phase transition from the dilute loop phase to the DS phase. The peak positions are nearly the same as various bond dimensions. When h z is further decreased, a hump appears around h z ≈ −0.470 and gradually becomes a peak with increasing bond dimension. After this hump, the model enters into a dense loop phase. By fitting the scaling relation of the entanglement entropy 38 at the point (h x , h z ) = (0.5, −0.8) inside the dense loop phase, we estimate an effective central charge c * ≈ 2.00 as shown in Fig. 5d. Considering that the exactly fixed point E is inside this region and is described by the non-unitary CFT with a central charge c = −7 × 2, we can conclude that the gapless dense loop phase is characterized by the non-unitary CFT with an effective central charge c * = 2.
Moreover, the correlation length ξ is also calculated along the h x -axis and shown in Fig. 5b. As h x is increased, a peak appears around h x = 0.785 and is enhanced by larger bond dimensions.
For h x > 0.785, the system enters into a special phase, where the largest eigenvalues of the transfer matrix are given by a complex conjugate pair and the resulting many-body phase spontaneously breaks the PT-symmetry 11,13 . To clearly see the PTSB phase, several leading eigenvalues of the transfer matrix along the circle line h = 0.85 are displayed in Fig. 5c. For θ > 0.19π, the model is in the gapped dilute loop phase, in which the dominant eigenvalues of the transfer matrix are real and positive. At θ ≈ 0.19π, the arguments of the leading eigenvalues reveal an exceptional point, where the corresponding eigenvectors coalesce. After this exceptional point, the arguments of a pair of complex dominate eigenvalues vary continuously in the PTSB phase. Moreover, the absolute value of the leading eigenvalue (red line) shows a cusp at θ ≈ −0.07π between the PTSB phase and dense loop phase, indicating a discontinuous change of the free energy. For θ < −0.07π, the dominate eigenvalue becomes real in the dense loop phase and the arguments of the second and third largest eigenvalues become ±2π/3, which are related to the momenta carried by the corresponding eigenvectors.
With these numerical results, a global phase diagram of the generic DS wavefunction can be fully established. From our numerical results, the CF line describes a continuous phase transition without any symmetry changes, while the BC line is a continuous phase transition with spontaneous PTSB. But both phase transitions are characterized by a non-unitary CFT with an effective central charge c * ≈ 1.01, and one of the fitting is displayed Fig. 5e. The corresponding central charge may correspond to the c = −2 non-unitary CFT. Meanwhile, the phase transition line AB is described by the same non-unitary CFT with the central charge c = −6/5 as the critical point A, and the line BH is the same as the two-dimensional Ising phase transition. More importantly, it should be mentioned that the continuous PTSB transitions are composed of the exceptional points, where the corresponding leading eigenvectors coalesce. On the other hand, the zeros of the partition function are distributed along the line AL, and the transition line CL corresponds to discontinuous ones.

Discussion
It should be emphasized that, the three non-topological phases adjacent to the DS phase in Fig. 4b are associated with the physics of two decoupled two-dimensional triangular lattice Ising spin model Eq. (11). On the decoupled line h = 1, the single-layer spin model allows us to perform exact diagonalization on a larger lattice size. For the two-column transfer operator TT t with a periodic boundary condition, the arguments of the eigenvalues related to the lattice momenta vanish. As shown in Fig. 6a, the absolute values of the eigenvalue spectrum reveal a discontinuous phase transition point L at (θ ≈ −0.05π), where many eigenvalues of the transfer matrix cross and the partition function changes sign. In Fig. 6b, the arguments of the eigenvalue spectrum have clearly demonstrated the presence of the exceptional points at A and E, where the spectra become completely real. This is not surprising because the exactly solvable non-Hermitian PT-symmetric models described by non-unitary CFTs always have the entire real spectra 39 . Since the transfer matrix cannot be transformed into a Hermitian operator, the nature of the critical points A and E is non-Hermitian and non-unitary. Although both A and E are exceptional points, there is an essential difference: the dominant eigenvectors coalesce at the point A, while the coalescing eigenvectors at point E are not the dominant ones. And the logarithms of moduli of the eigenvalues cross at both the points L and A, which are two boundary points of the PTSB phases.
Compared to the phase diagram of the TC phase, distinct physics is introduced by the intrinsic signs in the DS wavefunction. Since the bosonic anyon excitations of the DS phase are related to the charge anyons of the TC phase, the gapped dilute loop phase corresponds to the Higgs phase, while the dense loop phase corresponds to the confining phase even though it is a gapless phase and described by a non-unitary CFT. But the PTSB phase is an emerging non-unitary phase, where the auxiliary Ising spins oscillate spatially and the pairs of the bosonic anyons condense. This emerging phase may be caused by the breakdown of Mermin-Wagner theorem, which usually forbids spontaneous breaking of continuous symmetries in dimensions d ≤ 2. However, when the unitarity is broken in statistical systems with negative Boltzmann weights, the Mermin-Wagner theorem is no longer valid. The properties of the PTSB phase deserves further investigation. Moreover, there is a logarithmic attraction between the charge anyons along the transition line IJ because of its unitary CFT description, while along the transition line AB of the DS there is logarithmic repulsion between the pairs of bosonic anyons due to the negative scaling dimensions of the auxiliary spins from the nonunitary CFT. It should be noticed that both topological phase transitions out of the DS phase into the dense loop phase or the PTSB phase are non-unitary, while the phase transition into the dilute loop phase is still unitary. And the phase transition between the PTSB phase and the dilute loop phase is also non-unitary.
In summary, the quantum effects of the negative signs in a generic DS wavefunction have been explored using tensor network representation. The norm of this wavefunction is mapped to the partition function of a two-dimensional triangular lattice Ashkin-Teller model with imaginary magnetic fields and imaginary three-spin triangular face interactions. Using the CTMRG method, the global phase diagram has been fully established. Adjacent to the DS topological phase, we found a gapped dilute loop phase, a gapless dense loop phase described by non-unitary CFT, and a PTSB phase filled with zeros of the partition function in the thermodynamic limit. So our study has proved that the negative signs in the topological DS wavefunction are intrinsic and deeply connected to the negative Boltzmann weights of a PTsymmetric statistical model, shedding light on the understanding of topologically ordered phases of matter. Further, the intrinsic sign of a topologically ordered wavefunction can be understood as an obstruction to mapping a quantum topologically many-body state to a classical partition function. In two dimensions, there is a broad class of topologically ordered phases suffering such a sign problem, including twisted quantum doubles and even more generic string-net models 40 . The DS state is the first prototype one, and the similar investigations can be carried out to identify novel correlated many-body phases of matter with inherent quantum nature, whose physical properties cannot be reproduced by a partition function with positive Boltzmann weights.

Methods
The phase diagrams in Fig. 4 are determined from numerical tensor-network methods. The key challenge in the numerical tensor-network method is the contraction of the tensor network generated by Q. In general, the contraction of an infinitely large tensor network is implemented by finding the approximate environments with controllable errors. In our numerical calculation, we use the infinite matrix product state (iMPS) to approximate the environment for Hermitian system, and their main procedures are summarized in Fig. 7a, b. For a non-Hermitian system, the corner transfer matrix (CTM) is used to approximate the environment and their main ideas are explained in Fig. 7c-e.
The transfer matrix operator for the norm of the generic TC wavefunction is a Hermitian operator, so we can calculate its environments using iMPS with bond dimension D and optimize the iMPS using the VUMPS algorithm 19,36,37 in Fig. 7a. The (j + 1)th dominant eigenvalue λ j of the column-to-column transfer operator can be calculated from the local tensors of iMPS, as shown in Fig. 7b. Then the correlation length ξ can be obtained: ξ ¼ À1=log ðjλ 1 =λ 0 jÞ.
Since the transfer matrix operator for the norm of the generic DS wavefunction is non-Hermitian, its environments cannot be worked out using variational method. Instead, we solve the partition function using the CTM method, and the environment is approximated using edge tensors and corner tensors with bond dimension D. These tensors can be optimized using the CTMRG algorithm [16][17][18][19] in Fig. 7c. The correlation length can also be calculated, and the (j + 1)th dominant eigenvalue λ j is defined in Fig. 7d. The entanglement entropy S can be obtained by S ¼ À Tr ðρlog ρÞ, where the reduced density matrix ρ is defined in Fig. 7e.
We should emphasize that the CTMRG calculations do not converge when all dominant eigenvalues are complex. However, due to the double layer structure of the partition function, it can be derived that the most dominant eigenvalue of the double layer transfer operator is always real and positive. In the PTSB phase, the single layer transfer operator has a pair of dominant eigenvalues |λ|e ±iϕ , but the equal moduli dominant eigenvalues of the corresponding double-layer transfer operator are given by |λ| 2 and |λ| 2 e ±i2ϕ , where |λ| 2 has two-fold degenerate. When the couplings between two layers are turned on, there are small splits between the Fig. 6 The spectrum of the two-column transfer operator TT t for the single layer spin model. With the periodic boundary condition, we calculate the spectrum of the single layer spin model with the circumference L y = 12, where λ j is the (j + 1)th dominant eigenvalue and θ is the spin angle. a The red lines represent Àlog ðjλ j =λ 0 jÞ, where each red line corresponds to an eigenvalue, and the spectrum becomes nearly continuous in the filled areas. The black dashed lines correspond to the phase transition points L and A in the phase diagram Fig. 4b. b The blue dots are the arguments of the spectrum argðλ j Þ that become nearly continuous in the filled areas. The vertical dashed lines indicate two exactly solvable points E and A in the phase diagram Fig. 4b. Fig. 7 Main procedures of the numerical tensor-network methods. a In the variational uniform matrix-product-state (VUMPS) algorithm, the environment (light blue background) of an infinite tensor network is approximated by the infinite matrix-product-state (iMPS), which consists of local tensors (light blue squares) with a bond dimension D. The iMPS is optimized by maximizing the tensor network with the variational method. b The eigen-equation T j j > ¼ λ j j j > of the column-to-column transfer operator T is approximated by the local tensors, where the green box is the eigenvector j j > with the (j + 1)th largest eigenvalue λ j . c In the cornertransfer-matrix renormalization group (CTMRG) algorithm, the environment (light blue background) of an infinite tensor network is approximated by the edge fixed point tensors (light blue squares) and corner fixed point tensors (light blue snip single corner squares) with a bond dimension D. d The eigen-equation T 2 j j > ¼ λ j j j > of the column-tocolumn transfer operator T 2 is approximated by the edge tensors and local double tensors. e The tensor network representation of the reduced density matrix ρ in the CTMRG algorithm, where Z is the normalization factor. The red dashed line emphasizes that the up/down cut bonds are regarded as the row/column of the reduced density matrix. moduli of these dominant eigenvalues, leading to one positive dominant eigenvalue. Thus the CTMRG calculations can slowly converge in the PTSB phase.

Data availability
The data that support the findings of this study are available from the corresponding author on request.