O(3) nonlinear sigma model in 1+1 dimensions with matrix product states

We numerically study the spectral properties, the entanglement and the zero-temperature phase structure at nonvanishing chemical potential of the O(3) nonlinear sigma model. Using matrix product states, a particular kind of one-dimensional tensor network state, we show that we are able to reach the asymptotic scaling regime and to reproduce the analytical predictions for the mass gap at vanishing chemical potential. In addition, we study the scaling of the entanglement entropy towards the continuum limit obtaining a central charge consistent with 2. Moreover, our approach does not suffer from the sign problem and we also explore the phase structure of the model for nonzero chemical potential and map out the location of the transitions between different charge sectors with high precision.


I. INTRODUCTION
Nonlinear sigma models in 1 þ 1 dimensions with their spin/rotor degrees of freedom have a long history as interesting models in their own right and as benchmark models for four-dimensional gauge theories. Analogies include a negative beta function, which leads to the phenomenon of asymptotic freedom as well as nonperturbative effects like mass generation and topology. Moreover, the sign problem that hampers importance sampling at nonzero chemical potentials [1] also occurs in sigma models. Dual variables have been shown to solve the sign problem in CP(N − 1) as well as OðNÞ sigma models [13], and subsequent simulations of the O(3) model [14,15] revealed a second-order phase transition where the lowest mass in the spectrum equals the chemical potential, and with a dynamical critical exponent consistent with 2 as well as two-particle phase shifts in agreement with the analytical S-matrix [16].
In this study we numerically examine the O(3) model in one spatial dimension in its Hamiltonian lattice formulation. One of our motivations is that in the Hamiltonian framework the introduction of a chemical potential does not cause additional difficulties in numerical simulations; it simply multiplies the Hermitian charge operator. Thus, we can directly address the zero-temperature physics, e.g., the ground state in each sector of fixed charge. We note in passing that in principle also nonzero temperature physics can be addressed in the Hamiltonian setup. A difficulty arises from the fact that the Hilbert spaces per site, which we represent in angular momentum eigenfunctions (spherical harmonics), are infinite dimensional. Hence, we have to truncate them to a finite dimension, which we do by restricting the total angular momentum at each site. Nevertheless, the total dimension grows exponentially with the volume, i.e., the number of sites, thus rendering numerical approaches dealing with the full Hilbert space infeasible. Here we use matrix product states (MPS)-a particular kind of one-dimensional tensor network (TN) state-which efficiently parametrize a subspace of weakly entangled states. In particular, MPS techniques allow for computing ground states and low-lying excitations for many physically relevant Hamiltonians, and the effort only grows polynomially in the size of the system and the tensors [17,18].
In this work we apply MPS to the asymptotically free O(3) nonlinear sigma model and explore its ground state and the mass gap. Compared to previous TN studies of the model [48][49][50][51][52], we show that we are able to reach the asymptotic scaling regime, and, even for modest truncations, we do reproduce the analytical predictions, in particular in the scaling regime towards the continuum limit. This is rather nontrivial since in this regime large angular momenta are not suppressed, which will eventually render our truncations insufficient to faithfully describe the physics of the model when the gap closes exponentially.
Furthermore, we have access to the (bipartite) entanglement entropy in the ground state, and we investigate its scaling towards the continuum limit. As a consequence of the asymptotic freedom, the gap closes in this limit and the model becomes critical. We confirm the theoretically predicted logarithmic divergence with the correlation length [53] and extract the value for the central charge which we find to be consistent with 2.
Taking advantage of the fact the MPS methods do not suffer from the sign problem, we also explore the zerotemperature phase structure of the model at nonvanishing chemical potential. We observe the expected transitions between the different charge sectors of the Hamiltonian and are able to precisely locate the transition points.
The rest of the paper is organized as follows. In Sec. II we introduce the Hamiltonian lattice formulation and the observables we are using, and describe our MPS approach. Subsequently, we present our numerical results in Sec. III. Finally we discuss our findings in Sec. IV. Technical details of our basis and numerical extrapolations are provided in Appendixes A and B.

II. MODEL AND METHODS
The simplest representation of the O(3) model is by its continuum Euclidean action, where n is a real three-component unit vector, g is the dimensionless (bare) coupling and ν is summed over one spatial dimension and Euclidean time extending from 0 to the inverse temperature. As usual, the partition function is given by a path integral of expð−SÞ over n-field configurations. On a lattice with spacing a, the two-dimensional integral over the derivative terms is replaced by a 2 times the ðx; νÞsum over hopping terms −2nðxÞnðx þ aνÞ=a 2 (approaching −n∂ 2 ν n in the limit of vanishing lattice spacing), with the prefactor β replacing 1=g 2 . Thanks to asymptotic freedom, the continuum limit is achieved by β → ∞.
The action with chemical potential coupled to rotations in the first two components of n is obtained by replacing ∂ ν → ∂ ν − δ ν;0 μτ 2 with τ 2 being the second Pauli matrix (see, e.g., Ref. [54]) and, thus, is no longer real.

A. Hamiltonian lattice formulation
To study the model numerically with MPS, we switch to the Hamiltonian formulation. It can be deduced from the action above, but is also well known in the literature [55][56][57][58][59][60][61]. The potential term, corresponding to the spatial gradient, is discretized as described above. The kinetic term per site is that of a three-vector with fixed radius, i.e., momentum squared without radial component, which is nothing but the square of the orbital angular momentum. Hence, the lattice Hamiltonian for a system with N sites reads In the expression above, L 2 k is the angular momentum operator acting on site k (note the inverse β in this kinetic term), and μ the chemical potential coupling to the total charge Q which is nothing but the sum of third components of the angular momenta. Notice that contrary to conventional Monte Carlo approaches to the Lagrangian formulation we work with open boundary conditions for convenience in our MPS simulations [62], and the sum for the potential term only ranges to N − 1. Equation (2) describes a linear chain of coupled quantum rotors, and operators acting on the same site fulfil the well-known commutation relations The Hamiltonian conserves Q because the kinetic term obviously commutes with L z at every site, and the commutator with the hopping term yields ½n α k n α kþ1 ; þ ðα⇌γÞÞ, which vanishes as a result of the antisymmetry of the Levi-Civita symbol.
A suitable basis for such a system is the tensor product of the common eigenfunctions jlmi of L 2 and L z for each site (whose spatial representations are the spherical harmonics), Note that the commutation relations from Eq. (4) imply that the value of l is not bounded from above and hence the basis is infinite dimensional. While the first two Hamiltonian terms are diagonal in this representation, the potential is more complicated. Defining the combinations n AE ¼ ðn x AE in y Þ= ffiffi ffi 2 p , we can rewrite the hopping part as n þ k n − kþ1 þ n − k n þ kþ1 þ n z k n z kþ1 . In Appendix A we give the matrix elements hlmjn AE jl 0 m 0 i and hlmjn z jl 0 m 0 i in terms of Wigner-3j symbols. Selection rules are such that involved total angular momenta obey jl 0 − lj ¼ 1, while m 0 agrees with m AE 1 or m. In the Hamiltonian these expressions have to be used at neighboring sites k and k þ 1.
As we elaborate below, these states are very useful in the strong-coupling limit at small β, where the angular momentum term dominates. Near the continuum at large β, however, the hopping term tends to align the rotors. For this picture, independent angular momentum states are far from being effective, and it would be very useful to find a more suitable basis.

B. Observables and expectations
In the Hamiltonian framework, the charge Q is simply the sum over all m quantum numbers along the chain. The mass gap is essentially given by the energy difference between the first excited and the ground state, am ¼ aΔE Here we take into account the correction η due to the anisotropy of spatial and temporal couplings in the Hamiltonian formalism [63], to be able to make contact with the continuum prediction. From the one-loop renormalization one finds η¼1=ð1−1= βπÞ. In the Lagrangian formulation, Eq. (1), with periodic boundary conditions the continuum limit of the mass gap m in units of the lattice spacing a reads [59,63] am where we have used the two-loop expression for the Callan-Symanzik beta function to compute Λ L . An advantage of the MPS approach to be used is that it allows for easy access to the entanglement entropy of the system. Dividing the system into two contiguous subsystems A and B, we can efficiently compute the entanglement entropy for the bipartition which is given by the von Neumann entropy of the reduced density matrix In the expression above jψi is the state of the system and tr B refers to the partial trace over the subsystem B. In our simulations, we analyze the entropy in the ground state and choose the subsystem to be half of the chain. For a fixed value of β and vanishing μ the Hamiltonian from Eq. (2) is local and gapped. Thus we expect the entanglement entropy S for a subsystem A to grow with its surface area [17]. Since for a one-dimensional chain the latter consists of two points, S A is expected to saturate upon increasing the size of A.
Moreover, from Eq. (6) we see that as we approach the continuum limit β → ∞ the gap closes, and the system becomes critical. Hence, in the asymptotic scaling regime we expect the entropy for the reduced density matrix describing half of the system to diverge logarithmically as S ¼ ðc=6Þ logðξ=aÞ þk, where ξ=a is the correlation length in lattice units and c the central charge of the underlying conformal field theory describing the critical point andk a (nonuniversal) constant [53]. Using the fact that the correlation length is inversely proportional to the gap, ξ=a ∝ 1=am, we get together with Eq. (6) where k sums up the constant terms.
C. Spectrum at strong coupling In the limit of small β the system prefers zero angular momentum. Therefore, the ground state has energy E 0 ¼ 0. The first excited state (at μ ¼ 0) has l ¼ 1 at a single site and thus E 1 ¼ 1=β þ OðβÞ, which is also the leading order gap at small β (the next to leading order correction is −2β=3 [55,56]).
If we enforce a nonzero charge, the angular momentum cannot vanish everywhere, since at least a single site has to have a nonzero quantum number m and, thus, the associated l must be nonzero, too. It is not hard to see that it is energetically favorable to induce a certain charge Q (smaller than or equal to N) by assigning the minimal quantum numbers ðl; mÞ ¼ ð1; 1Þ to Q sites and ðl; mÞ ¼ ð0; 0Þ to the complementary ones [64]. The leading kinetic term is not sensitive to the spatial arrangement of those occupied sites, which yields ð N Q Þ degenerate "unperturbed eigenstates", but the subleading hopping term knows about neighbors, and thus breaks this degeneracy: Let L ¼ f0; 1g denote the relevant quantum numbers ðl; mÞ ¼ fð0; 0Þ; ð1; 1Þg. Then the matrix elements of the hopping term, hL k L kþ1 jn k n kþ1 jL 0 , if an excitation is created at site k and annihilated at site k þ 1, or vice versa. This can be read as a Hubbard model.
For periodic boundary conditions, it is possible to analyze the strong-coupling limit further analytically, at least for a small number of sites. In the fixed charge sector Q ¼ 1, for instance, degenerate perturbation theory on the L ∈ f0; 1g states reveals that the ground state is the sign coherent superposition ðj10…0i þ j010…0i þ Á Á Á þ j0…01iÞ= ffiffiffiffi N p with entanglement entropy log 2. At Q ¼ 2 the ground-state coefficients increase with the distance between the two occupied sites. For a small number of sites N ¼ 4, 6, 8, 10 we have determined the half-chain entropies to grow like 1.202, 1.300, 1.331, 1.344.

D. Phase structure for nonvanishing chemical potential
The Hamiltonian conserves the total charge Q; hence it is block diagonal and each block can be labeled with the corresponding eigenvalue q of the charge operator. For convenience, let us rewrite Eq. (2) as where aW aux sums up the potential term and the kinetic part, and is independent of μ. If we now restrict aH to a block characterized by q, the charge operator is simply given by q1 and the Hamiltonian further simplifies to From this equation we can see that the ground-state energy inside a block with charge q is given by with aE 0 ½aW aux j q being the minimum eigenvalue of the operator aW aux restricted to block q and independent of aμ. Thus, inside a sector the ground-state energy scales linearly with the chemical potential and the slope is given by q. In particular, for zero charge the energy does not depend on aμ at all, which is also called Silver blaze property [68]. Moreover, since Q is proportional to the identity inside a specific charge sector, the ground state for each block is independent of aμ and simply given by the ground state of aW aux j q . While the energies aE 0;q ðμÞ can be measured in our approach, in the physics of grand-canonical ensembles only the chemical potential is fixed, and the system minimizes these energies over all charge sectors. The global minimum of the energy selects a certain charge sector that depends on aμ. At certain values aμ c the energy levels of two sectors with different charges cross and it is energetically favorable to go from one sector with charge q to another one with chargeq, which generically is the successor,q ¼ q þ 1. The value of aμ c can be determined by equating aE 0;q ðμ c Þ ¼ aE 0;qþ1 ðμ c Þ, which yields Thus, we expect the total charge of the ground state to exhibit discontinuous changes as we increase the chemical potential.
In particular, the first transition, from 0 to unit charge, is expected to happen at the energy gap, aμ c ¼ aΔE (related to the mass am as discussed in Sec. II B). This follows from the equation above, if in the absence of μ, i.e., for aW aux , the ground state has Q ¼ 0, aE 0 ½aW aux ¼ aE 0 ½aW aux j 0 , whereas in the first excited state the massive particles form a triplet and Q ¼ 1, being the lowest state in that sector, aE 1 ½aW aux ¼ aE 0 ½aW aux j 1 . As a result, the right-hand side of Eq. (12) is nothing but the energy gap. The next transition can be used to extract two particle energies in finite volumes and from them phase shifts [14].
In the infinite-volume limit all these transitions merge to one curve for the charge density as a function of μ, and only the first critical μ c ¼ m survives to mark-at zero temperature-a second order quantum phase transition [15]. Similar observations have been made in the O(2) model at small temperature [69].

E. Truncation in angular momentum
In our Hamiltonian, Eq. (2), the spectrum of L 2 (and L z ) is unbounded; thus the Hilbert space even for a single lattice site is infinite dimensional. In any numerical approach working on a finite basis one therefore has to apply a truncation. We choose to keep all basis states at all sites with Together with all possible m quantum numbers this yields a local Hilbert space with dimension d ¼ Note that the normalization of n and the commutator relations from Eq. (4) are violated when representing n in such a truncated basis. As we show in Appendix A, this only concerns expectation values where the highest sector with l max is involved. Nevertheless, the truncated Hamiltonian still conserves the total charge Q. Thus, for large enough values of l max the violations are expected to be negligible, which will be checked a posteriori.

F. Numerical approach
In order to compute the low-lying spectrum of the Hamiltonian from Eq. (2) we use the MPS ansatz. For a system with N sites on a lattice with open boundary conditions the ansatz reads is a Ddimensional complex row (column) vector. The parameter D, called the bond dimension of the MPS, determines the number of variational parameters in the ansatz and limits the amount of entanglement that can be present in the state (see Refs. [70][71][72][73] for detailed reviews). For our system, i k should be read as a super index for ðl k ; m k Þ with the summation ranges discussed above.
In our simulations we are interested in the ground state and the energy gap of the model. The MPS approximation for the ground state can be found variationally by minimizing the energy. To this end, one iteratively updates the tensors M i k k , one at a time, while keeping the others fixed [18]. The optimal tensor in each step is found by computing the smallest eigenvalue [74,75] of the effective Hamiltonian describing the interactions of site k with its environment. After obtaining the ground state, the first excited state can be found in a similar fashion by projecting the Hamiltonian on a subspace orthogonal to the ground state and running the same algorithm with the projected Hamiltonian [21,76].

A. Spectral properties
Let us first focus on the case of vanishing chemical potential. In order to benchmark our approach, we study the ground-state energy and the energy gap for a range of values of β ∈ ½0.65; 1.8 and different system sizes N ¼ 40, 80, 120. In addition, to probe for the effects of truncating the local bases to a finite angular momentum, we explore l max ¼ 1, 2, 3, 4. Moreover, we have another source of error due to the limited bond dimension available in the numerical simulations. This error can be controlled by repeating the calculation for every combination of ðβ; l max ; NÞ for several bond dimensions D ∈ ½80; 160 and extrapolating to the limit D → ∞ (see Appendix B for details on the extrapolation procedure).
In Fig. 1 we show the extrapolated results for the groundstate energy density aE 0 =N for various values of l max and N. In general, we observe that the errors due to the finite bond dimension in our simulation are negligible and there is almost no dependence on the system size, as the values for aE 0 =N for N ¼ 40, 80 and 120 are essentially identical. For the simplest nontrivial truncation, l max ¼ 1, we see that the ground-state energy density notably differs from those obtained for larger values. In contrast, there is hardly any difference between the values obtained with l max ≥ 2, only for larger values of β the results for l max ¼ 2 deviate slightly from those for l max ¼ 3, 4.
While the ground-state energy density is fairly insensitive to the truncation and does not show strong finite-size effects, the situation is noticeably different for the energy gap. Figure 2 reveals that there is a significant difference between results for different truncations with l max ≤ 3, and only for the largest two values of l max our results are in agreement. Although finite-size effects are small for l max ¼ 1, 2, we see from Figs. 2(a) and 2(b) that our numerical values for the gap in those cases are not compatible with the asymptotic scaling predicted by Eq. (6). These data rather seem to approach a constant value as the slope decreases with increasing values for β. On the contrary, for l max ¼ 3, 4 there are much stronger finite-size effects. While for all our system sizes the data seem to enter the asymptotic scaling around β ≈ 1.  system size, N ¼ 120, we recover the asymptotic scaling up to the largest value of β ¼ 1.8 we study, as can be seen in Figs. 2(c) and 2(d). In addition, we observe that in this case our error bars are growing as the value of β increases. This is a direct consequence of the asymptotic scaling: since the gap closes exponentially with increasing β and we use the same range of D values in all our simulations, the MPS approximations we obtain are becoming progressively worse. Nonetheless, the errors are reasonably small up to β ¼ 1.6 and we reliably recover the asymptotic scaling between 1.2 ≤ β ≤ 1.6.
B. Scaling of the entanglement entropy towards the continuum limit Taking advantage of the fact the we have easy access to the entanglement entropy in the ground state, we can also explore its scaling towards the continuum limit. Our results for the gap show that only for our largest system size the finite-size effects are small enough that we recover the asymptotic scaling regime over such a large range of values for β. Hence, for the following we focus on N ¼ 120 and study the half-chain entropy and its scaling as we approach the continuum limit. Our results are shown in Fig. 3. For l max ¼ 3, 4 we clearly observe an almost linear scaling for large β, as predicted in leading order by Eq. (8), and there is hardly any difference between the data for both truncations. The entropies for l max ¼ 2 still show an approximately linear behavior towards the continuum limit; however, the values for larger β are significantly smaller than those obtained for l max ¼ 3, 4. For completeness, we also show the data for the simplest nontrivial truncation, l max ¼ 1. As one can see, these do not show the linearly divergent behavior in β but rather seem to approach a constant value, thus giving an indication that the model does not become critical as we approach the continuum. In particular, this is also in agreement with our observation in Fig. 2(a) that the gap is not closing.
Fitting our data for l max ≥ 2 to the theoretical prediction, we can extract the central charges, too. Looking again at our results for the gap in Figs. 2(c) and 2(d), we see that we enter the asymptotic scaling regime around β ¼ 1.2 and that our errors are reasonably small up to β ¼ 1.6. Hence, we choose this range to fit our data for the entropy to Eq. (8). Additionally, we estimate our systematic error by comparing the results from different fitting intervals (see Appendix B for the details). The results for the central charge obtained by this procedure are shown in Table I. For l max ¼ 2, we see that our data points for β > 1.6 are progressively below our fit, thus indicating that the slope further decreases as we go closer to the continuum limit. Hence, our value of c ¼ 1.66 for that case only seems to be an upper bound. Similar to l max ¼ 1 this might give an indication that the model does not become critical, consistent with our observation in Fig. 2(b) that the gap does not follow Eq. (6) but rather seems to tend to a constant value. On the contrary, for l max ¼ 3, 4 our fit describes the data in the entire range for β we study and we obtain values of c ≈ 2. Note that OðNÞ models in the perturbative regime, for which the nonlinearity does not play a role, describe N − 1 massless bosons. These turn into N massive bosons at low energies for N ≥ 3. Thus, the measured central charge of the ground state points to the two perturbative degrees of freedom.
These observations together with the gap provide a comprehensive picture of the effect of the truncation. Although the ground-state energy density is rather insensitive to the truncation in angular momentum, we observe that for l max ¼ 1, 2 the gap does not seem to close, and the model does not become critical for β → ∞, in particular for l max ¼ 1. Only for l max ¼ 3, 4 we capture the relevant features and recover the theoretical prediction for the asymptotic scaling of the gap as well as the divergence in the entropy according to Eq. (8) as we approach the continuum limit.

C. Phase structure in the presence of a chemical potential
Contrary to conventional Monte Carlo methods, our MPS approach does not suffer from the sign problem and  2.04 AE 0.14 we can also investigate the phase structure of the model at nonvanishing chemical potential. To this end we explore the ground-state energy and the total charge as a function of aμ for β ¼ 1.0, 1.1, 1.2 and N ¼ 40, 60, 80, which is around the beginning of the asymptotic scaling for aμ ¼ 0, but still far away from regimes with noticeable finite-size effects. Again, we study several truncations l max ¼ 1, 2, 3, 4 and bond dimensions D ∈ ½80; 240 to estimate the errors of the truncation as well as our numerical errors. Figure 4 shows an example of the total charge in the ground state as a function of the chemical potential. Indeed, we see the theoretically predicted discontinuous changes in the values of q by one unit as aμ is growing. In addition, the inset in Fig. 4 shows that the behavior of the ground-state energy is in excellent agreement with the theoretical prediction from Eq. (11), and we observe a linear scaling with a slope given by the total charge of the sector.
Similar to the case of vanishing chemical potential, we can extrapolate our results for the ground-state energies to the limit D → ∞, and, using Eqs. (11) and (12), we can estimate the exact location of the transitions (see Appendix B for details). Our results for the various truncations, system sizes and couplings are shown in Fig. 5. In general, we observe that our data do not show very strong truncation effects and there is hardly any difference between results with l max ≥ 2. Only for the simplest nontrivial truncation, l max ¼ 1, the locations of the transitions are significantly shifted towards higher values of the chemical potential. Comparing data with different values of the coupling, we see that the transitions between different charge sectors occur for smaller values of aμ c as we go closer to the continuum, independently of the truncation and the volume. Moreover, for a fixed value of β the widths of the plateaus in each charge sector shrink with increasing system size, consistent with the expectation that in the FIG. 4. Charge as a function of the chemical potential for N ¼ 80, β ¼ 1.2, l max ¼ 4 and D ¼ 200. Inset: Corresponding ground-state energies for data points in the sector for q ¼ 0 (red triangles), q ¼ 1 (green squares) and q ¼ 2 (purple diamonds). The solid lines represent the theoretical prediction according to Eq. (11), where the slope is fixed by the charge value and aE 0 ½aW aux j q has been determined from the numerical data. infinite-volume limit the plateaus for q > 0 merge, and there is only a single second-order quantum phase transition.
In particular, independent of the truncation the first transition should occur at aμ c ¼ aΔE. In Table II we compare the values for both quantities obtained in our simulations for N ¼ 80. Indeed, the critical aμ c at the first transition is in almost perfect agreement with the energy gap. We recapitulate that the latter has been obtained at the same parameters ðβ; l max ; NÞ but with aμ ¼ 0 by an independent computation, determining the lowest energy in the Hilbert space orthogonal to the ground state.

IV. CONCLUSION
We have studied the spectral properties, the entanglement entropy in the ground state and the zero-temperature phase structure at nonvanishing chemical potential for the O(3) nonlinear sigma model using MPS. To render the Hilbert space finite dimensional, we have restricted the angular momentum at each site to a certain maximum value l max .
Looking at the ground-state energy density at vanishing chemical potential, we observe that it is rather insensitive to the truncation, and does not show strong finite-size effects. We find that l max ≥ 2 is already sufficient to avoid noticeable effects even in the weak-coupling regime. In contrast, the energy gap is much more sensitive to the truncation.
For l max ¼ 1, 2 the gap does not show the expected asymptotic scaling, but rather seems to approach a constant value as we go towards the continuum limit, independently of the volume. Increasing l max to 3, 4 we observe that for large enough volumes the gap closes exponentially in the weak-coupling regime, as predicted by perturbation theory.
Investigating the scaling of the bipartite entanglement entropy in the ground state towards the continuum limit, we observe a similar picture as for the mass gap. For l max ¼ 1 we do not recover the expected divergence as we approach the continuum limit and the model does not become critical, consistent with the observation that the mass gap does not close according to the theoretical prediction. For l max ¼ 2 our data do not allow us to fully rule out that the energy gap closes and the model becomes critical. However, the central charge obtained in this case is significantly smaller than those for larger l max and the deviations become more pronounced the closer we go to the continuum limit. On the contrary, for l max ¼ 3, 4 our data for the entropy show the predicted divergence, and the central charges obtained in those cases are consistent with 2.
Our MPS approach does not suffer from the sign problem and readily allows us to investigate the zerotemperature phase structure at nonvanishing chemical potential. We clearly observe the theoretically predicted transitions between different charge sectors of the Hamiltonian, and we can determine the location of the transitions with great precision. Here truncation effects are only clearly visible for l max ¼ 1, and the transition points are significantly shifted towards larger values of the chemical potential compared to results with l max ≥ 2, which show hardly any difference among each other. With increasing volume the transition points between phases with charges q ≥ 1 come closer to each other, consistent with the expectation that in the infinite-volume limit only single second-order quantum phase transition survives. In particular, we have verified that the first transition occurs at chemical potential values equal to the gap, independently of the truncation and the coupling.
Our findings suggest that even moderate values of l max are sufficient to capture the relevant features of the model, at least for the ground state and the phase structure at nonvanishing chemical potential. Similar observations have been made in MPS studies for O (2) and O(4) rotor models [30] as well as for simple (1 þ 1)-dimensional gauge models [29,32]. This is especially encouraging for potential future quantum simulators for rotor models [78]. Although current quantum simulation experiments for lattice field theories [79][80][81][82] are rather limited, one might nevertheless be able to study relevant phenomena even with a modest amount of resources. Moreover, our results can serve as a test bench for these experiments.
Concerning the phase structure of the model, more features of the ground state at high chemical potentials can be investigated. As the system tends to become planar [15] it leans towards the O(2) model, but the conjectured Berezinskii-Kosterlitz-Thouless transition of the latter [54] has not been seen in numerics [15]. Secondly, the induced particles form a massive triplet (with known S-matrix [16]), which could modify entanglement properties and consequently the central charge.
On the technical side, it would also be very intriguing to find another basis for the weak-coupling regime such that truncation effects become much milder. This would allow a much improved investigation of the continuum limit. Finally, it would be interesting to probe the potential effects of spatially twisted boundary conditions which are crucial for resurgence conjectures of sigma models [8][9][10][11][12].

APPENDIX A: BASIS FORMULATION AND TRUNCATION TO A FINITE DIMENSION
Here we show how to calculate the matrix elements for the components of the operator n in the angular momentum eigenbasis. The calculation can be simplified by realizing that the operators n AE , n z in position representation are related to the spherical harmonics with l ¼ 1, which in turn form a spherical tensor operator of rank 1, as follows: 0 Thus the matrix elements can be computed with the Wigner-Eckart theorem and the expectation values become integrals over three spherical harmonics, where we have used that Y Ã l;m ¼ ð−1Þ m Y l;−m . The integral is related to Wigner-3j symbols (which are proportional to Clebsch-Gordan coefficients) by a well-known formula and we obtain The right-hand side is nonzero if and only if the conservation m ¼ m 0 þ M holds and jl 0 − lj ¼ 1 (as a consequence of the triangle relation, however, l 0 ≠ l, from the first symbol). The expectation values needed in this work are obtained by identifying n AE ⇌ ∓ X AE1 and n z ⇌X 0 . As discussed in the main text, Eq. (13), we truncate the local Hilbert space through limiting l by l max at every site. Consequently, the normalization of n and the commutator relations from Eq. (4) will be violated, as we demonstrate now. The truncation means that we represent all operators by finite dimensional matrices, where the n α -expectation values have been calculated above. For bilinears of n α we have to take the matrix product; e.g., the normalization becomes Obviously, if l 0 were summed up to infinity, the second line upon completeness would simply yield δ l 1 ;l 2 δ m 1 ;m 2 . From above we know that l 0 s up to maxðl 1 þ 1; l 2 þ 1Þ contribute to such sums of expectation values, but they are neglected in our truncation once l 1 or l 2 equal l max . Thus, the normalization is violated in the highest sectors of angular momentum. This violation does not affect lower l-sectors or operators made of just L's, since their expectation values are nonzero only within the same l-sector; so larger l 0 are irrelevant anyhow. That the truncated Hamiltonian still conserves the total charge Q follows from a similar argument. Looking at the derivation below Eq. (4), we have to check the commutator of n α with L z at some site (the neighboring site being a spectator) in the truncated representation. In the corresponding expectation values, as in the second line of Eq. (A6), the L z factor again limits the relevant l 0 to l 1;2 , which are kept such that the completeness is intact and ½L α ; n β ¼ iε αβγ n γ holds indeed.

APPENDIX B: NUMERICAL DETAILS
Here we discuss some details of the evaluation procedures for the numerical data and how we estimate our errors for the data shown in the main text.
1. Ground-state energy, mass gaps and ground-state entropy As we already mentioned in the main text, for every combination of system size N, coupling β and truncation l max we repeat the numerical simulation for a range of bond dimensions D ∈ ½80; 160. Subsequently we can estimate the exact value and the numerical error for the ground-state energy and the mass gap by extrapolating to limit D → ∞.
To this end, we extrapolate our data points for aE 0 and aΔE with the largest three bond dimensions linearly in 1=D (see Fig. 6 for an example). As an estimate for the central value we take the mean value of our data point with the largest bond dimension and the extrapolated value. The error is estimated as half of the difference between those two values.
For the half-chain entropies of the ground state we proceed in the same fashion (see Fig. 7 for an example). Similar to the ground-state energy, we observe that our results are well converged in bond dimension and the numerical errors due to limited values of D in our simulations are negligible.
The errors for the central charges in Table I are estimated by comparing the results of different fits. As explained in the main text, the central value is obtained by fitting our data to Eq. (8) in the interval 1.2 ≤ β ≤ 1.6. Subsequently we repeat the same fit in every possible interval ½β min ; β max with 1.2 ≤ β min < β max ≤ 1.8, which contains at least six data points. As an estimate for the error we take the maximum absolute value of the difference between the results from these fits and our central value.

Extracting the location of the transitions
The location of the transitions can be determined with the help of Eq. (12) from the main text. In our simulations we measure the ground-state energy and the charge. Hence, we can extract aE 0 ½aW aux j q using Eq. (11). Combining these two equations, we find for the location of two adjacent phases with q þ 1 and q aμ c ¼ aE 0;qþ1 ðμÞ þaμðq þ 1Þ − ðaE 0;q ðμÞþaμqÞ; ðB1Þ whereμ andμ refer to the values of the chemical potential at which we determined the constants aE 0 ½aW aux j qþ1 and aE 0 ½aW aux j q for the two phases. Notice that Eq. (B1) allows us to extract the location of the transition between two adjacent phases as long as we have a single data point in each of them available. Thus, this method is more efficient than extracting the transition points directly from the discontinuities of Q as we do not need a high resolution in μ to determine them precisely. The value of q can in principle be determined exactly, since the total charge is still an exact symmetry of the truncated Hamiltonian. In practice, finite bond dimension effects can break this symmetry and it is only restored for large enough values of D, as the example in Fig. 8 shows. Hence, in our simulations we start with D ¼ 80 and keep increasing the bond dimension until the deviations of q from integer values are negligible. For all the cases we study we find that a bond dimension of 240 is sufficient to stabilize the values of q and to allow us to unambiguously determine in which phase we are [see Fig. 9(b) for an example].  The exact values for the ground-state energy and the corresponding errors are again estimated analogous to the case of vanishing chemical potential by extrapolating linearly in the last two data points [see Fig. 9(a) for an example]. From our extrapolation we obtain an estimate for the exact ground-state energy as well as the systematic error δaE 0;q . Propagating this systematic error in Eq. (B1), we obtain the systematic error for the transition points, which is given by δaμ c ¼ jδaE 0;q ðμÞj þ jδaE 0;q ðμÞj: ðB2Þ In general, we observe that this systematic error due to the finite bond dimension in the simulations is rather negligible and for all data shown in the main text the error bars are smaller than the markers (see also Table II in the main text).