Tensor network simulation of the (1+1)-dimensional $O(3)$ nonlinear $\sigma$-model with $\theta=\pi$ term

We perform a tensor network simulation of the (1+1)-dimensional $O(3)$ nonlinear $\sigma$-model with $\theta=\pi$ term. Within the Hamiltonian formulation, this field theory emerges as the finite-temperature partition function of a modified quantum rotor model decorated with magnetic monopoles. Using the monopole harmonics basis, we derive the matrix representation for this modified quantum rotor model, which enables tensor network simulations. We employ our recently developed continuous matrix product operator method [Tang et al., Phys. Rev. Lett. 125, 170604 (2020)] to study the finite-temperature properties of this model and reveal its massless nature. The central charge as a function of the coupling constant is directly extracted in our calculations and compared with field theory predictions.


I. INTRODUCTION
The (1+1)-dimensional nonlinear σ-model (NLSM) has played important roles in both high energy and condensed matter physics. The NLSM shares various common features with the (3 + 1)-dimensional non-Abelian gauge theories, such as the asymptotic freedom [1], dynamical generation of mass gap [2], solitons [3,4], and nontrival θ vacua. The NLSM can be further generalized to the 1/N expandable CP N−1 model [5,6] which is believed to be relevant to the study of the strong CP problem [7]. Hence, a thorough understanding of the nature of the NLSM can undoubtedly give much insight into the study of the non-Abelian gauge theories in 3 + 1 dimensions.
From the condensed matter side, the (1+1)-dimensional NLSM arises in the context of the Haldane's conjecture [8,9]: the (1+1)-dimensional O(3) NLSM with θ = 0 (θ = π) is the low-energy, long-wavelength effective theory for the quantum antiferromagnetic Heisenberg chain with integer (halfinteger) spin S . Since the NLSM with θ = 0 is known to have exponentially decaying correlations [2], the integer-spin antiferromagnetic (AF) Heisenberg chain is conjectured to be gapped. Meanwhile, based on known results for the spin-1/2 AF Heisenberg chain, Haldane conjectured [9] that the O(3) NLSM with θ = π topological term is a massless theory. Later, by mapping the NLSM with θ = π term to a modified quantum rotor model, Shankar and Read [10] claimed that this model should be massless for all values of the coupling constant. In the strong-coupling limit, as pointed out by Affleck and Haldane [11,12], the O(3) NLSM with θ = π is equivalent to the SU(2) 1 Wess-Zumino-Novikov-Witten (WZNW) conformal field theory (CFT) [13][14][15], whose central charge is c = 1. In the weak-coupling limit, this model corresponds to two massless bosons, thus has a central charge c = 2. In the renormalization group framework, the NLSM with θ = π flows from the unstable fixed point at the weak-coupling limit * wanglei@iphy.ac.cn † hong-hao.tu@tu-dresden.de to the stable fixed point at the strong-coupling limit, and, according to c-theorem, the central charge varies monotonically between these two limiting cases. From the numerical side, in the study of lattice field theories, the Monte Carlo algorithm has been a standard approach since the beginning of this field. However, in many cases, the Monte Carlo approach is hindered by the sign problem-more specifically, for example, in the case of NLSM, the straightforward Monte Carlo simulation encounters the sign problem when the θ-term is nonzero. Although several approaches, such as the Meron-cluster Monte Carlo algorithm [16][17][18] and the analytic continuation approach based on imaginary θ simulation data [19][20][21][22], have been successfully developed to overcome the sign problem for this specific case, these methods are rather specific and cannot be easily extended and applied to other systems. On the other hand, in recent years, the tensor network algorithms have achieved rapid development, and have been increasingly applied to the numerical simulation of lattice field theories [23,24]. Unlike the Monte Carlo approach, the tensor network methods are free from the sign problem, and thus can hopefully be applied to many problems where the Monte Carlo simulation are hindered or even prohibited. It is then meaningful to develop and test tensor network algorithms for the NLSM with θ-terms.
In this work, we perform a tensor network simulation of the NLSM with θ = π topological term. Inspired by Ref. [10], we work in the Hamiltonian formulation and map the NLSM with θ = π term to a modified quantum rotor model where the quantum rotors are decorated with magnetic monopoles. By representing the modified quantum rotor model in the basis of magnetic monopoles, we obtain its matrix representation which automatically enables MPS-based simulations of this model. Making use of the recently proposed continuous matrix product operator (cMPO) method [25], we simulate the finite-temperature properties of the modified quantum rotor model and present clear numerical evidence for its massless nature. Moreover, we also obtain the central charge as a function of the coupling constant, and compare the result with the field theoretical predictions. This paper is organized as follows. In Sec. II, we introduce the Hamiltonian formalism, i.e., the quantum rotor model for the NLSM, respectively for θ = 0 and θ = π cases. We also introduce the matrix representation for the quantum rotor models. In Sec. III, we introduce the cMPO approach and its application in the quantum rotor model. In Sec. IV, we show the numerical results which prove the massless nature of the NLSM with θ = π. Finally, Sec. V summarizes the results and provides some outlook. In Appendix A, we discuss the matrix representation of the modified quantum rotor model introduced by Sec. II in detail. In Appendix B, we present a proof for the mapping between NLSM and quantum rotor models for both θ = 0 and θ = π cases. Appendix C includes some details in the numerical simulation.

II. HAMILTONIAN FORMULATION
In (1+1) dimensions, the Euclidean action of the O(3) NLSM is given by where n is a unit vector that rotates in the three-dimensional space, g is the dimensionless bare coupling constant, and x, τ represent the Euclidean space coordinates. Due to the possible existence of instantons in this model, one can extend the action by adding a θ-term where θ is periodic in 2π. The partition function of this model is written in the path integral formulation as Z = Dn exp(−S [n]), where, throughout this work, the functional integration Dn is defined with respect to the real unit vector field. In this section, we will briefly review the Hamiltonian formulation of the NLSM with θ = 0 and θ = π.
A. Hamiltonian formulation for NLSM with θ = 0 For the NLSM with θ = 0, it has been well established that the Hamiltonian formulation is given by the one-dimensional quantum O(3) rotor model on the lattice [10,[26][27][28][29] whereL j andn j respectively represent the angular momentum operator and rotor operator on site j, a is the lattice spacing, and K > 0 is a constant. The operators satisfy the following commutation relations: where µ, ν, λ = x, y, z. In the low-energy, long-wavelength limit, the field theoretical description of this rotor model is just the O(3) NLSM with θ = 0 and 1/g 2 = K.
The eigenstate of the rotor operatorn is parametrized by continuous angle variables, which is not convenient for tensor network simulations. Instead, a discrete basis is preferred, for which the eigenbasis of the angular momentum operators serves as a natural choice:L 2 |l, m = l(l + 1)|l, m and L z |l, m = m|l, m . The quantum numbers l and m take integer values, l = 0, 1, 2 . . . and m = −l, −l + 1, . . . , l. In this basis, the kinetic term in the Hamiltonian (3) becomes diagonal. The rotor couplings can be rewritten asn i ·n j = ν∈{0,±}n ν in −ν j withn ± = (n x ± in y )/ √ 2 andn 0 =n z , whose matrix representation in the angular momentum basis can be obtained by takingn ν as spherical tensor operators [29]. With the matrix representation of the Hamiltonian, it is then straightforward to represent the partition function as a tensor network or perform ground/excited state calculations via MPS-based methods.
B. Hamiltonian formulation for NLSM with θ = π In the presence of a nonvanishing θ-term, it is a nontrivial task to incorporate it in the quantum rotor model formulation. As pointed out in Ref. [10], by adding a magnetic monopole with magnetic charge q = 1/2 at the center of the rotor and setting nearest-neighbor couplings to be antiferromagnetic, the low-energy effective theory becomes a NLSM with θ = π. In the Hamiltonian formulation, the presence of the magnetic monopole modifies the definition of the angular momentum operatorL, and the Hamiltonian becomes [10] In spatial coordinates, the modified angular momentum operator is defined byL = n × (−i∇ − A) − n, where n is the unit vector pointing at the direction of the rotor, and A is the vector potential describing the magnetic field generated by the magnetic monopole (see, e.g., Ref. [30] for more details). The low-energy physics of this model is described by the NLSM with θ = π, whose coupling constant satisfies 1/g 2 = K.
To find the matrix representation of the Hamiltonian (5), we make use of the eigenbasis of the modified angular momentum operator -the monopole harmonics [30,31]. It is known thatL still satisfies the angular momentum commutation relations [30], and the eigenbasis of monopole harmonics can be labeled by well-defined angular momentum quantum numbers (l, m), which satisfies (L ) 2 |q, l, m = l(l + 1)|q, l, m and (L ) z = m|q, l, m . Here, the quantum number q denotes the magnetic charge at the center of the rotor and hence takes the value q = 1/2, which distinguishes itself from the ordinary spherical harmonics with q = 0. In the presence of the magnetic charge q = 1/2, l and m can only take half-integer values, l = 1/2, 3/2, 5/2, . . . and m = −l, −l + 1, . . . , l. Making use of this angular momentum eigenbasis, the matrix representation of the Hamiltonian in Eq. (5) is similar to the case with θ = 0 -the kinetic term is diagonal, and the matrix representation for the rotor coupling term can be evaluated by using the properties of spherical tensor operators, the details of which are included in Appendix A. Based on the properties of the monopole harmonics, we also provide a proof for the mapping between the lattice Hamiltonian and the continuous theory of the NLSM with θ = π in Appendix B.
In practical simulations, we need to truncate the physical Hilbert space at each site. The most natural choice is to choose a maximally allowed angular momentum quantum number l max and drop the states with l > l max . For the Hamiltonian in Eq. (5), one can infer that this truncation scheme is effective only when the constant K is small, i.e., near the strong coupling limit. An interesting limit is K → 0, where one can choose l max = 1/2, and the modified quantum rotor model reduces to the S = 1/2 antiferromagnetic Heisenberg chain. For large values of K, in principle, one has to use large enough l max to obtain quantitatively accurate results.

III. TENSOR NETWORK APPROACH TO THE MODIFIED QUANTUM ROTOR MODEL
From the numerical side, we make use of the recently developed cMPO method [25] to study the finite temperature properties of the (modified) quantum rotor model in Eq. (5). The reason for using this approach is twofold. First, since the theory is expected to be massless for all choices of the coupling constant, working at the finite temperature can help reduce the requirement on the bond dimensions compared to ground-state simulations in the thermodynamic limit [32,33]. It also allows us to adjust the temperature for different choices of the coupling constant. Second, the cMPO approach works in the continuous time limit, which automatically eliminates the discretization error in the imaginary time direction.
In this section, we will briefly review the cMPO approach, and introduce the cMPO formulation for the modified quantum rotor model defined in Eq. (5). We will also discuss two key properties of the cMPO for this model: (i) Hermiticity, which enables a direct global optimization during the simulation; (ii) Symmetry, which is inherited from the Hamiltonian and allows us to reduce the computational cost in numerical simulations.

A. Brief review of the cMPO approach
The cMPO approach is based on the observation that there exists a compact MPO representation for the infinitesimal time evolution operator exp(− Ĥ ) when we only consider up to the first order in [34]. The neglected higher order terms of will not incur any errors since we will take the → 0 limit. From this MPO one can build the tensor network representation for the partition function Z = Tr e −βĤ (see Fig. 1). The local tensor T in the MPO can be expressed as where Q is an operator-valued scalar, I is the identity operator, L and R are operator-valued vectors (not to be confused with the angular momentum operatorL), and P is a matrix consisting of operators. The operators contained in Q, L, R, and P are operators acting on the physical Hilbert space, which all come from the Hamiltonian: Q corresponds to local terms, L and R encode nearest-neighbor interactions, and P comes from longer-range interactions. The physical dimension is thus the dimension of the physical Hilbert space at each site. The virtual bond dimension D = d + 1, where d is the dimension of vectors L and R. Next, as shown in Fig. 1, the tensor network for the partition function Z can be formed by stacking β/ layers of the MPOs together. In the thermodynamic limit, i.e., L → ∞, this tensor network can be efficiently contracted using the idea of transfer matrix [35][36][37]. The transfer matrix T refers to the column of tensors in the tensor network, which is also an MPO, and we can approximate its dominant eigenvector with an MPS, which is referred to as the boundary MPS (see Fig. 1). More specifically, here, as we take the continuous time limit → 0, the MPO representation for the transfer matrix becomes continuous (hence the name cMPO), and the corresponding boundary MPS becomes a continuous MPS (cMPS). Since the cMPO is uniform with the periodic boundary condition (along the imaginary time direction), it is natural to use a uniform cMPS |ψ , parametrized by the local tensor as the boundary cMPS. Similar to Eq. (6), Q ψ and R ψ in Eq. (7) correspond to an operator and a vector of operators, respectively. The operators contained in T ψ are parametrized by matrices, whose dimension is the bond dimension of the cMPS. From the boundary cMPS, one can further extract the thermodynamic properties of the system. To obtain the boundary cMPS |ψ , one can directly mini-mize the free energy density if the cMPO T is Hermitian. In general, T is non-Hermitian, and |ψ has to be optimized by the power method, i.e., by repeatedly acting T on a trial solution for |ψ and compressing its bond dimension. At each iteration step, the compression of T|ψ into the cMPS |φ (with a smaller bond dimension) is again a variational optimization process that maximizes the fidelity F = φ|T|ψ / φ|φ , where we have dropped a constant factor 1/ ψ|T † T|ψ for simplicity.

B. cMPO formulation for the modified quantum rotor model
For the modified quantum rotor Hamiltonian in Eq. (5), the local tensor T is given by One can easily identify the contents of Q, L, R, and P in Eq. (6). It is worth mentioning that this cMPO is Hermitian.
To see this, we can perform a unitary transformation described byÛ = exp(iπ(L ) y ] to all the operators contained in the local tensor T . Using the commutation relations betweenL andn, we find On the one hand, according to Eq. (10),Û switches the contents of L and R, and since P is empty, it effectively switches the left and right bonds of the local tensor T [see Eq. (6) and Fig. 1]. When viewed as a large matrix, the cMPO T becomes its own transpose (and also its own Hermitian conjugate, as T is a real matrix) after this unitary transformation. On the other hand, since the physical bonds of the local tensors in cMPO are all contracted (see Fig. 1), this unitary transformation is merely a gauge transformation and leaves the cMPO unchanged. Therefore, the cMPO T is Hermitian, which allows us to directly optimize the boundary cMPS by variationally minimizing the free energy in Eq. (8).
In our calculations, before variationally minimizing the free energy, we perform a few power method steps to obtain a good initialization for the variational optimization. To help stabilize the power method procedure, we apply the unitary transformation exp[iπ(L ) z ] on every second site, such that the rotor couplings in x and y directions in the Hamiltonian (5) become ferromagnetic [38], For this "rotated" Hamiltonian, the local tensor now reads The Hermiticity of the cMPO can be proven in an analogous way.
C. U(1) symmetry of the cMPO Although the Hamiltonian in Eq. (11) is SO(3) symmetric, we shall just use its U(1) subgroup (i.e., rotational invariance around the z axis) in our tensor network simulations. Mathematically, the U(1) symmetry is generated bŷ where ν = 0, ±. Equation (13) not only proves the U(1) invariance of the Hamiltonian, but also indicates that the U(1) symmetry can be encoded into the cMPO. Combining Eq. (12) and Eq. (13), we get T U( ) where, the same as Fig. 1, the vertical and horizontal lines respectively correspond to the physical and virtual indices of the local tensor T . U(θ) is the matrix representation of U(θ). On the virtual bonds, V(θ) = exp(iθZ), where Z = diag(0, 1, −1, 0). As shown in Fig. 1, the physical indices of local tensors in the cMPO are all contracted, and the left-hand side of Eq. (14) is thus a gauge transformation that leaves the cMPO invariant. Meanwhile, the right-hand side of Eq. (14) gives rise to a global U(1) rotation, which is described by V(θ) on each virtual bond. Therefore, we can take V(θ) as a symmetry transformation of the cMPO. Based on this fact, we can correspondingly construct a U(1)-symmetric boundary cMPS, which contains block structures and helps lower down the computational cost. The construction of the U(1)invariant boundary cMPS follows the general rules to construct symmetric MPS [39][40][41][42]. More specifically, we build the boundary cMPS from the local tensors satisfying up to a phase factor. Here, T ψ represents the cMPS local tensor, and U ψ (θ) denotes the U(1) rotation acting on the internal indices of the cMPS. The degeneracy sectors on the cMPS vertical bond are determined dynamically, the details of which are included in Appendix C 1.

IV. RESULTS
In this section, we describe our numerical results, which provide numerical evidence of the massless nature of the NLSM with θ = π. More specifically, we demonstrate the results for the free energy density and the bipartite entanglement entropy of the boundary cMPS and compare them with the predictions of CFT.
In our numerical simulation, we calculate the finite temperature properties of the modified quantum rotor model. We perform simulations from K = 1.0 to K = 6.0 in order to cover a fairly large range of values for the coupling constant. The range of temperatures varies with K. For each choice of K, the temperature T ranges from K/300 to K/60, since the energy scale increases with K. The maximal angular momentums are chosen to be l max = 1/2, 3/2, 5/2. The bond dimension of the boundary cMPS ranges among χ = 12, 18, 24, 30, and we extrapolate the results to an infinite bond dimension (see Appendix C 2 for details).
Our code implementation is publicly available at [43].

A. Universal correction to free energy
For one-dimensional quantum systems described by CFT, a universal finite-size correction to the free energy appears at low temperature [44,45] where f 0 is the free energy density at zero temperature, c is the central charge, and v is the effective "velocity of light" in the theory. Equation (16) predicts that the specific heat is linear in T . This property is in sharp contrast to that of gapped systems, where the free energy manifests an exponential scaling at low temperatures. From our calculation results, we verify the massless nature of the system by comparing with Eq. (16). Figure 2 shows the free energy density as a function of (T/K) 2 for different choices of K and the linear fitting of the data to Eq. (16). The linear fitting is performed within the range K/300 ≤ T ≤ K/100, the detailed results of which are shown in Table I. For all values of K that we consider, we observe a clear linear dependence of the free energies on T 2 , which perfectly coincides with the prediction of Eq. (16). For small values of K, the results quickly converge with respect to l max . As K increases, the results for different l max 's gradually deviate from each other, since the angular momentum basis becomes less effective as the system approaches the weak coupling limit. Nonetheless, they show a tendency of convergence and still serve as qualitative evidence for the massless nature of this model.

B. Temporal entanglement entropy
For systems described by CFT, the temporal direction and the spatial direction are equivalent, and the boundary cMPS can be viewed as the ground state of the "temporal Hamiltonian" which is described by the same CFT. Therefore, the bipartite entanglement entropy of the boundary cMPS is another important indicator of the massless nature of the system, as it satisfies [46][47][48][49]    where c is the central charge, and S 0 is a nonuniversal term. By fitting the entanglement entropy with respect to the inverse temperature β, we can extract the central charge of the system. The bipartite entanglement entropy in the uniform temporal cMPS is calculated as follows. We divide the cMPS by an equal bipartition, and calculate the von Neumann entropy between the two subsystems A and B (see Fig. 3), From the bipartition, the boundary cMPS can naturally be written as where m, n denote virtual indices through the cuts. To transform {|φ A(B) mn } to an orthonormal basis, we consider the matrix M A(B) m n ,mn = φ A(B) m n |φ A(B) mn and its eigen-decomposition One can easily verify that and the reduced density matrix is given by Combining this equation with Eq. (18), one can evaluate the entanglement entropy. In Fig. 4, for different values of K and l max , we show the bipartite entanglement entropy of the boundary cMPS as a function of ln(Kβ) and the linear fitting of the data to Eq. (17). The linear fitting results are shown in Table II. As the parameter K increases, results for different l max gradually deviate from each other, which is similar to the results of the free energy. Nonetheless, these results still show a tendency to converge as l max increases, and, at least qualitatively, from these results we can confirm the linear relation between the bipartite entanglement entropy and ln β, which verifies Eq. (17). Furthermore, with Eq. (17), we can obtain the central charge c from the linear regression of the bipartite entanglement data. As pointed out in Refs. [10,11], in the renormalization group framework, the NLSM with θ = π flows from the unstable fixed point at the weak-coupling limit (g = 0 and central charge c = 2) to the stable WZNW fixed point at the strong-coupling limit (g = ∞ and c = 1), and the central charge varies monotonically between the two fixed points according to the c-theorem. Figure 5 shows the numerical results of the central charge. For the special case of l max = 1/2, the system reduced to the antiferromagnetic Heisenberg model, which is well known to be described by the SU(2) 1 WZNW model [11,50]. As shown in Fig. 5, for l max = 1/2, and for all values of K, the numerical result indeed gives c ≈ 1, and it can get closer to c = 1 if one pushes the calculation to lower temperatures. For l max = 3/2, 5/2, the numerical result for central charge shows a monotonic tendency with respect to K and approaches c ≈ 1 as K becomes small. As is similar to the previous results, the result reaches good convergence at l max = 5/2 when K is small (K ≤ 3.0), but only serves as qualitative evidence for larger values of K.

V. SUMMARY AND OUTLOOK
In summary, we have numerically studied the (1+1)dimensional O(3) nonlinear σ-model with θ = π term using our recently developed cMPO method. We work in the Hamiltonian formulation-the modified quantum rotor model decorated with magnetic monopoles, and derive its matrix presentation in the monopole harmonics basis. From this matrix representation, we obtain the cMPO representation of the modified quantum rotor model, and study its finite-temperature properties with the cMPO method. We calculate the free energy density of the system and the bipartite entanglement entropy of the boundary cMPS, and compare their scaling with the predictions of CFT, from which we confirm the massless For tensor network simulations, we have to truncate the monopole harmonics basis and only consider the states with relatively small angular momentum. This truncation scheme works well near the strong coupling limit since the high angular momentum sectors are suppressed, and the result can easily get converged with the maximal angular momentum l max . However, as the system goes toward the weak coupling limit, the angular momentum is no longer suppressed, and the monopole harmonics basis becomes less effective and moderate values of l max often do not lead to a converged result. An intermediate improvement based on our formalism is to implement the non-Abelian symmetry in tensor network simulations [51][52][53] so that one can push the calculation to larger l max . For future investigations, it is also important to look for a better basis near the weak coupling limit in order to truncate the local physical Hilbert space in a more effective way.
Another interesting direction is to consider the possibility of a Hamiltonian formulation for the NLSM with θ 0 or π. The current formulation relies on the monopole harmonics basis, which is tied to the quantization of the magnetic monopole charge and does not seem to have a straightforward generalization to θ ∈ (0, π) cases. Undoubtedly, the capability to efficiently simulate the (1+1)-dimensional NLSM with a full range of θ will open a broad way to better understand this paradigmatic field theory.
In this Appendix, we derive the matrix representation for the modified quantum rotor model in Eq. (5). Apparently, the kinetic term is diagonal in the angular momentum basis q, l, m| and we only need to focus on the coupling term between neighboring rotors. Then i ·n j term can be represented aŝ wheren ± = (n x ±n y )/ √ 2. The matrix representation ofn ± andn z can be evaluated by noticing their relations with the spherical tensor operators [29]. More specifically, we havê whereŶ 1,M (M = 0, ±1) is the spherical tensor operator of rank 1. The matrix elements of these operators can be computed as integrals over monopole harmonics [31] q, l 1 , where Appendix B: Proof for the mapping between the quantum rotor model and NLSM In this Appendix, we present a proof for the mapping between the quantum rotor model and the NLSM based on the path-integral formalism. Our proof covers both θ = 0 and θ = π cases.

θ = 0 case
For the θ = 0 case, we write the partition function of the ordinary quantum rotor model in Eq. (3) as a path integral where the imaginary time τ ∈ [0, β] is divided into N intervals, i.e., ∆τ = β/N, and n (k) ≡ (n (k) 1 , n (k) 2 , . . . , n (k) L ) represents the rotor configuration at the kth time interval. The periodic boundary condition n (N) = n (0) along the imaginary time direction is imposed. Substituting the explicit form ofĤ into Eq. (B1), we obtain Within Eq. (B2), we look at n (k+1) j |e −∆τL 2 j /2Ka |n (k) j first. By inserting a complete set of angular momentum basis into it, we get The summation over m can be handled by the addition formula of the spherical harmonics, where P l is the Legendre polynomial, and γ (k) j represents the angle between n (k) j and n (k+1) j . In the ∆τ → 0 limit, γ (k) j becomes infinitesimal, and Therefore, where we have introduced ∂ τ n (k) j = γ (k) j /∆τ. From Eq. (B6) to Eq. (B7), the summation over l is carried out by the Euler-MacLaurin formula.
Next, we evaluate exp[∆τ(K/a) L j=1 n (k) j · n (k) j+1 ] in Eq. (B2). Here we assume that the lattice spacing a is small, and n (k) where δ (k) j represents the angle between n (k) j and n (k) j+1 . Along this line, we get Substituting this equation into Eq. (B1), and taking the limit a, ∆τ → 0, we find where Dn(x, τ) ≡ lim a→0 lim ∆τ→0 j,k (Kae ∆τK/a /2π∆τ), and the Lagrangian density L is given by For the NLSM with θ = π term, similarly as Eq. (B1), we write the partition function of Eq. (5) as where n (k) ≡ (n (k) 1 , n (k) 2 , . . . , n (k) L ) denotes the kth configurations of the quantum rotors, and n (N) = n (0) . At each time slice, we have n (k+1) |e −∆τĤ |n (k) which is similar to Eq. (B2) except the modified angular momentum operator and the antiferromagnetic coupling between the neighboring rotors.
In Eq. (B14), we first focus on the kinetic terms n (k+1) j |e −∆τ(L j ) 2 /2Ka |n (k) j , which can be evaluated as where |q, l, m represents the eigenbasis of the modified angular momentum operator (where q = 1/2). Using the addition formula of the monopole harmonics [30,31], the summation over m can be carried out, and we get l m=−l n (k+1) j |q, l, m q, l, m|n (k) where P µ,ν n is the Jacobi polynomial, γ (k) j represents the angle between n j (k) and n j (k + 1). Ω (k) j is the area of the spherical triangular formed by n j (k), n j (k + 1), and the north-pole axis on the unit sphere. As γ (k) Substituting Eq. (B17) into Eq. (B14), and performing the summation over l with the Euler-MacLaurin formula, we arrive at The continuum in the imaginary time direction can be straightforwardly taken. To take the continuum limit in the spatial direction, we split the n (k) j into the slowly varying part m (k) j and the rapidly varying part l (k) j , where a is the lattice spacing. Substituting Eq. (B20) into terms in Eq. (B19), we get /a, and the higher-order of a are neglected. The surface term Ω (k) j requires a more careful treatment. Since the neighboring rotors tend to have an antiparallel alignment, we separate the summation of surface terms in pairs. We look at each of these pairs which corresponds to the surface area of the ribbon formed by the trajectories of n (k) 2r and −n (k) 2r−1 on the unit sphere. We calculate ∆S (2r) as Substituting Eqs. (B21), (B22), (B24) into Eq. (B19), and integrating out the fast fields l (k) j , we get Finally, taking the continuum limit, we arrive at where the Lagrangian density is given by In the cMPO approach, the boundary cMPS is uniform, and has periodic boundary condition (PBC). For the U(1)symmetric boundary cMPS, one has to determine the degeneracy sectors on the vertical bonds, which further determines the block structures in the cMPS local tensors. However, unlike the open-boundary MPS simulations, it is difficult to determine an optimal choice for these degeneracy sectors. In our simulation, we dynamically determine the degeneracy sectors during a power method process.
As is discussed in Sec. III B, we initialize the variational optimization of the boundary cMPS by a few steps of power method process. To start with, we construct an initial cMPS by using the cMPO tensors at the boundary sites. Using the right boundary cMPS as an example, the local tensor T (0) ψ of the initial cMPS |ψ (0) can be obtained from the local tensor T [see Eqs. (6) and (7)], The U(1) degeneracy sectors for |ψ (0) are simply derived from those for the cMPO. Next, we repeatedly act the cMPO T on the cMPS, and compress the cMPS variationally if the bond dimension exceeds the target bond dimension χ. By considering the truncation scheme illustrated in Fig. 6, we initialize the variational compression process and determines the degeneracy sectors. Suppose we are at the nth power method step, and we act the cMPO T on the cMPS |ψ (n) . We construct a reduced density matrix ρ (n) by cutting one horizontal bond of T|ψ (n) in the tensor network corresponding to ψ (n) |T † T|ψ (n) , as shown in Fig. 6 (a) and (b). The reduced density matrix ρ (n) also has a block-diagonalized structure due to the U(1) symmetry, and it shares the same degeneracy structure with T|ψ (n) . After diagonalizing ρ (n) , by keeping χ eigenvalues with the largest abstract values, we can obtain an isometry and insert it into T|ψ (n) , which can be used as the starting point for the variational compression [see Fig. 6 (c) (d)]. The degeneracy sectors for ψ (n+1) can be automatically determined during this process. Finally, after a few power steps described above, we can further optimize the boundary cMPS by variationally minimizing the free energy globally.

Extrapolation to the infinite bond dimension
For each set of parameters, the bond dimension for the boundary cMPS ranges among χ = 12, 18, 24, 30. We eliminate the bias brought by the finite χ and extrapolate the data to the infinite bond dimension.
Our extrapolation procedure follows that given in Ref. [29]. First, we perform a linear fitting with respect to 1/χ for