The classical two-dimensional Heisenberg model revisited: An $SU(2)$-symmetric tensor network study

The classical Heisenberg model in two spatial dimensions constitutes one of the most paradigmatic spin models, taking an important role in statistical and condensed matter physics to understand magnetism. Still, despite its paradigmatic character and the widely accepted ban of a (continuous) spontaneous symmetry breaking, controversies remain whether the model exhibits a phase transition at finite temperature. Importantly, the model can be interpreted as a lattice discretization of the $O(3)$ non-linear sigma model in $1+1$ dimensions, one of the simplest quantum field theories encompassing crucial features of celebrated higher-dimensional ones (like quantum chromodynamics in $3+1$ dimensions), namely the phenomenon of asymptotic freedom. This should also exclude finite-temperature transitions, but lattice effects might play a significant role in correcting the mainstream picture. In this work, we make use of state-of-the-art tensor network approaches, representing the classical partition function in the thermodynamic limit over a large range of temperatures, to comprehensively explore the correlation structure for Gibbs states. By implementing an $SU(2)$ symmetry in our two-dimensional tensor network contraction scheme, we are able to handle very large effective bond dimensions of the environment up to $\chi_E^\text{eff} \sim 1500$, a feature that is crucial in detecting phase transitions. With decreasing temperatures, we find a rapidly diverging correlation length, whose behaviour is apparently compatible with the two main contradictory hypotheses known in the literature, namely a finite-$T$ transition and asymptotic freedom, though with a slight preference for the second.


Introduction
The two-dimensional classical Heisenberg model has moved into the focus of attention as early as in the 1960ies [1,2] as a paradigmatic model for understanding magnetism in spin systems. The presumably most striking result achieved at that time was the Mermin-Wagner theorem: In systems featuring continuous symmetries in one and two spatial dimensions, there cannot be an instance of spontaneous symmetry breaking at any non-zero temperature [3]. As such, this statement applies to the classical Heisenberg model, too. Interestingly, at the same time, work was performed that seemed instead to indicate the existence of a phase transition at finite temperature [4,5], building on high temperature expansions and spin-wave arguments. This apparent contradiction naturally sparked a significant interest in the community, with different authors employing different techniques to answer this question. This includes studies using Monte Carlo techniques [6][7][8][9][10][11], high temperature expansions [12,13], and other methods with contradictory results [14,15]. Meanwhile, two other important aspects of the puzzle emerged, both with profound implications on our modern understanding of field theories, including quantum ones. The first is the fact that the absence of spontaneous symmetry breaking strictly speaking does not imply the absence of quasi-long-ranged order, and so-called Berezinskii-Kosterlitz-Thouless (BKT) transitions [16,17] may occur. Some studies have indeed found signatures resembling those of a BKT phase transition at T c ∼ 0.6 [18]. More recent works, using accurate finite-size scaling, seemed rather to suggest a pseudo-critical region [19] and quasi-long-range order [20]. The second one relies on the mapping between D-dimensional classical models and (d + 1)-dimensional quantum ones, with D = d + 1 [21][22][23]. The classical two-dimensional Heisenberg model should then share properties with the quantum field theory known as the O(3) non-linear sigma model [21][22][23], with a bare coupling directly proportional to the classical temperature (see Section 2 for more details). The most striking feature of the latter is the non-perturbative phenomenon of asymptotic freedom [24,25]: Namely the effective coupling constant being vanishing at large energies (short distances), while tending to get stronger and stronger at low energies (large distances), similarly to what happens in the celebrated case of 3 + 1-dimensional quantum chromodynamics (QCD) [26][27][28]. According to this line of thinking, there should then be no doubt about the absence of whatsoever phase transition at finite temperature, with the essential physics being captured by an exponentially diverging correlation length while approaching T = 0 (again, see Section 2). However, such mappings to quantum field theories are strictly valid only in the continuum limit, with corrections originating from a finite lattice spacing possibly giving rise to additional effects, which in turn could be considered misleading or interesting features [29]. All in all, it is fair to say that the specifics of the phase transition of the classical Heisenberg model in two spatial dimensions has remained inconclusive despite all these attempts. This state of affairs is aggravated by the central role this model plays in quantum statistical mechanics. In this work, we provide a fresh perspective to this long-standing question, by making use of state-of-the-art tensor network (TN) [30][31][32] approaches to address the problem at hand. Specifically, we express the partition function of the classical Heisenberg model in two dimensions as an infinite tensor network consisting of O(3)-symmetric rank-four tensors for a given inverse temperature β > 0. We achieve such a formally exact rewriting via an expansion over the fundamental characters of the (non-Abelian) group, in a similar fashion to what was recently performed for the Abelian O(2) case in Refs. [33,34]. We then compute a number of relevant quantities such as spin-spin correlators, the thermodynamic entropy S and geometric entanglement , among others. This requires contracting the whole infinite two-dimensional tensor network. We employ a symmetry-preserving corner transfer matrix renormalization group (CTMRG) approach [35][36][37][38] in order to achieve this, similarly to what has been done for the quantum case of infinite projected entangled pair states (iPEPS) [39]. Our techniques have the added value of directly tackling the thermodynamic limit: Residual finite-size effects are encompassed by the so-called bond dimension of the ansatz and could be accounted for in a rather systematic way [40][41][42][43][44]. These features have made two-dimensional tensor networks a very suitable tool for studying intricate condensed matter problems, not only via their ground states [45][46][47][48] but even beyond [49][50][51][52][53][54], as well as finite temperature properties of both classical and quantum models in two spatial dimensions [33,34,[55][56][57][58][59][60][61][62]. The present work builds upon and develops this substantial technical machinery. Taking profit from an SU (2)-invariant framework for tensor manipulations [63], we are able to achieve very large bond dimensions of the environment and thus to extrapolate our predictions in a presumably unprecedented way. To the best of our knowledge, this is the first study of the classical Heisenberg model in two dimensions using tensor network approaches and thus helps to shed further light onto the debated questions. The stable and efficient numerical machinery that we provide here also opens up a lot of exciting new opportunities to explore related models and materials, including ones in the realm of classical statistical mechanics.
This work is organized as follows. In Section 2, we define the model and we revise its diverse interpretations sketched above. In Section 3 we introduce the formalism to write the partition function as a tensor network. Following this construction, Section 3.3 describes how the infinite two-dimensional network can be contracted using a corner transfer matrix (CTM) scheme. Hence we not only compute fundamental thermodynamic quantities as the free energy and the entropy per lattice site (Section 3.4), but we also glance at quantum-inspired quantities such as the corner entropy and the geometric entanglement (Section 3.5) which could possibly highlight phase transitions. Moreover, as explained in Section 4, we get efficient access to a number of observables, both local and at a distance. In particular, we discuss the bond energy as well as spin-spin correlation functions, that unveil the characteristic behaviour of the Heisenberg model at different temperatures. Finally, our work is concluded in Section 5 with some discussions and an outlook on possible future applications.

The model
The classical Heisenberg model represents the workhorse of statistical mechanics for capturing (anti-)ferromagnetism when no lattice distortion or any other anisotropy source plays a major role. Consequently, the constituent objects are three-dimensional vectorsŜ k of unit length, |Ŝ k | 2 = 1, living on the sites k ∈ V and interacting along the links i, j ∈ E of a graph G = (V, E) capturing the lattice structure of the system. The defining Hamiltonian then reads evidently invariant under actions of the O(3) symmetry group, i.e., rotations and reflections. The central object for the study of statistical mechanics models is the partition function Z, from which thermodynamic quantities, as well as local expectation values and correlation functions can be computed. At inverse temperature β = 1/T > 0 (setting k B = 1), it reads where dΩ k := dθ k dφ k sin θ k is the infinitesimal surface element of the unit sphere at site k. We focus here on the apparently simple, yet very rich, square lattice in two spatial dimensions for the reasons we illustrated in the introduction and we deepen below in this section. Moreover, we deal here with the ferromagnetic version (FM) of the model, though we could equivalently work with the anti-ferromagnetic version (AFM) instead, since in the classical case the two models can be mapped into each other by a suitable transformation on bipartite lattices. We stress, however, that the entire tensor network construction of Section 3 is generally valid for any graph G and interaction sign, a property which discloses the path for future studies in different scenarios, e.g., frustrated lattices.
Before dwelling on the details of our tensor network approach, let us recollect an important strategy to uncover the physics of the model at hand, namely through the investigation of an associated field theory. Condensed matter theorists are indeed usually interested in the universal properties arising in the neighbourhood of a phase transition. There, the correlation length ξ diverges to values much larger than the lattice spacing a, or otherwise said the typical energy scale of excitations (gap) ∆ vanishes with respect to any bare energy scale. Microscopic lattice details should not matter anymore: Correlation and response functions can be then extracted from a field theory defined in the continuum without intrinsic cutoffs (the lattice instead gives an ultraviolet momentum cutoff Λ 1/a). The exponent involved in the partition function is then interpreted as the action S of the theory, with one of the original dimensions playing the role of the (imaginary) time τ , along which (Wick-rotated) Feynman path integrals are performed [21][22][23]. The most standard of these procedures consists in replacing the originalŜ k unit vectors with some N -component coarse-grained field ( x, τ ) → φ( x, τ ) and then getting out a φ 4 Ginzburg-Landau-Wilson functional: This is the realm to which the machinery of symmetry-breaking mechanisms or its forbiddance via the Hohenberg-Mermin-Wagner theorem finds application. However, an even more suggestive possibility -making use of vector fields ( x, τ ) → n( x, τ ) which preserve the essential property of having unit length -is offered by the O(N ) non-linear sigma model (NLSM) [21][22][23] with the spin-wave velocity set to unity, c = 1, and the bare coupling g = N a d−1 /β reading simply g = 3/β for the N = 3, d := D − 1 = 1 case of our interest. This simply formulated model shares a wealth of peculiarities with more intricate higher-dimensional ones, prominently with quantum chromodynamics (QCD) in d = 3. There, particle physicists want to make use of the lattice only as a convenient tool to investigate the field theory per se, i.e., they want to send the UV cutoff Λ and the energy scale J to infinite, while keeping all other scales (∆, ξ, etc.) fixed, which is usually achieved via renormalization group (RG) procedures. Along the decades, results first obtained in the large-N limit via formally exact saddle-point expansions of the action have been refined to all N ≥ 3. 1 The arguably most striking one is the so-called RG beta-function [21,24,25] which unveils the cornerstone phenomenon of asymptotic freedom. The meaning thereof becomes apparent when we integrate the RG beta-function to predict what happens at a given scale of energy, or in this framework of (quantum) temperature T q > 0, where g(Λ = 1/a) is the bare physical coupling we defined above at lattice scale. For T q → ∞, i.e., at large energies and for short distances, the effective coupling vanishes and the theory is almost free, while for T q → 0, i.e., low energies and large distances, the effective coupling gets stronger and stronger. This means that the physics resembles more and more the one of lower and lower β values, i.e., the system is always magnetically disordered. Correspondingly, a sizeable gap for excitations develops, a behaviour often called "dynamical mass generation" and quite simplistically sketched as ∆ + = c Λ exp −2πN/((N − 2)g) + O(g 0 ) , which however holds only asymptotically and is therefore misleading in trying to fit numerical data. Summarising, the field theory approach strongly indicates a divergent correlation length ξ = c/∆ + at low classical temperatures (i.e., high β and low bare g values), without any signature of a phase transition at finite values of it. However, all the mappings conducted above are strictly speaking valid only under the assumption of being close enough to the continuum limit (i.e., to the critical regime in classical terms), and lattice effects might give rise to additional effects. The latter are either treated as spurious by particle physicists or as interesting features by condensed matter ones [29]. Some observations are still in order before concluding this section and moving to our tensor network study and its results. The first one is that, by introducing the angular momentum operators L k as the canonical conjugate to the position of a particle on the unit sphere, we could establish a mapping with the so-called O(3) quantum rotor model, too, 1 We recall here that the cases N = 1 (Ising) and N = 2 (XY) are substantially different in that they exhibit a phase transition of symmetry-breaking and topological (BKT) nature, respectively.
where the bare couplings readg = 1/β 2 and J = β/a. We also introduced a quantum temperature T q which corresponds to the finite size of the system along the (imaginary) time direction in the previous two models. The second one is that an analogue construction could be started from a purely quantum Heisenberg model in d = 1, where the S k are now quantum spin operators, leading at first sight to the same result of an O(3)-NLSM field theory. However, along the path we should be particularly careful in at least two respects. First, in noticing that at quantum level the ferromagnetic case is pretty different from the anti-ferromagnetic one, namely it has a different dynamical critical exponent, z = 2, thus impairing a straightforward D = d + 1 classical analogy. Second, we should account for geometric phases (a.k.a., Berry phases and their generalizations) which could arise through the rotating spins and contribute to the action with a topological term with θ = 2πS and Q the so-called Pontryagin index. This is a topological invariant characterizing the configurations of n( x, τ ) in spacetime, i.e., the homotopy classes via π 2 (S 2 ) = Z. Such a θterm (often also dubbed Weiss-Zumino action term) is the one responsible for the peculiarly different behaviour, gapped vs. gapless, of integer vs. half-integer spin Heisenberg chains in d = 1, at the heart of the celebrated Haldane conjecture [64]. Noticeably enough, such a topological term is also the one generated in the semi-classical context by certain lattice distortions and plays a central role in the description of Skyrmions (i.e., the defects with |Q| = 1) in real materials [65][66][67][68]. Vice versa, in the absence of an explicit coupling, there exist quite compelling arguments to rule out such Skyrmions from playing the role of vortices in O(2) theories, i.e., being responsible for topological transitionsà la BKT [69].

The partition function as a tensor network
Here, we aim to use tensor networks as the numerical technique to compute the classical partition function. Tensor networks have been developed in the context of strongly-correlated quantum systems to represent many-body wave functions based on their entanglement structure. However, their network setup resembling the underlying physical lattice makes them equally well-suited for the simulations of classical models, by means of their partition function. In the following derivation, we will introduce a graphical convention that will conveniently lead to the tensor network structure of the partition function, resembling the underlying graph G = (V, E) the model is defined on. We here substantially contribute to the study of classical spin models by means of tensor networks, but remark that the earliest historical studies, that can be identified as being tensor network studies for systems in two spatial dimensions, have been dedicated to the study of classical statistical models [70]. In some way, therefore, we build on this mindset and tradition here. By employing the character expansion of the rotational group O(3), also known as the Rayleigh equation or plane-wave expansion in the contexts of electromagnetism and scattering theory, we can rewrite the Boltzmann factors on the edges E of the graph in terms of a linear combination of products of terms living on the vertices V according to Here β → I +1/2 (β) are modified Bessel function of the first kind and (θ, φ) → Y ,m (θ, φ) are the spherical harmonics. In Eq. (8), we also introduce the graphical notation for the Boltzmann factors, from which the full partition function can be graphically derived. The Boltzmann factors on the links E are drawn as an ellipse, which can be decomposed into a discrete sum of (products of) coefficients living on the vertices V . This rewriting, hence, transforms the integration over all site-variables (θ k , φ k ) in Eq. (2) into a product of independent integrations of all spherical harmonics living there. The explicit sum over and m is denoted by the line connecting the two lattice sites, for which we can define the following two objects The arrow on the edges is introduced in preparation for the symmetric tensor network ansatz we are going to develop. Note that the spirit of the mapping in Eq. (8) has also been used in previous works on the O(2) XY model [33,34]. The same decomposition can also be used in the anti-ferromagnetic version of the Heisenberg model, in which additional factors of (i) 2 appear in the plane-wave expansion of Eq. (8). The missing parts for the construction of the partition function in Eq. (2) are the integrations of the two continuous parameters (θ k , φ k ) on all vertices V of the lattice. In our graphical notation they will be denoted as Once a concrete graph for the lattice underlying the Heisenberg model is considered, the two parts of the graphical notation in Eq. (9) and Eq. (10) will directly lead to the tensor network structure, as we will demonstrate in the following sections.

Application to the one-dimensional Heisenberg chain
Let us start with the construction of the partition function for the classical Heisenberg spin-chain. This allows us to employ the graphical notation introduced above in a simple setting. Moreover, it will build an intuition for the tensor network appearing and the benefit of exploiting SU (2) symmetry therein. Assuming a linear chain of three-dimensional spinsŜ k , the full partition function becomes Exploiting the decomposition of the Boltzmann factors in Eq. (8), the spherical harmonics can be assigned to the vertices of the corresponding spins they act upon. The integration over the continuous variables can be performed locally, i.e., on every vertex independent of all the others. The partition function in Eq. (11) is then simply the product of objects which live on the vertices of the lattice. This unveils a fundamental property of the construction, namely that the local objects can be decomposed into two parts: i) A so-called degeneracy part that includes numerical factors stemming from the Boltzmann factors on the edges E, and depends on the temperature and on the angular momenta principal numbers j only; ii) a structural part that ensures that the angular momenta on the edges attached to each vertex are compatible. In particular, using the orthogonality of the spherical harmonics, we find for the two-edges vertices of a linear chain. For lattices with higher connectivity (i.e., a larger number of edges per vertex), the structural part will be Clebsch-Gordan coefficients (and combinations thereof) that dictate the coupling of spins, as we will demonstrate in the next section. At this point we can introduce the SU (2)-symmetric tensor network notation, in which the local objects in Eq. (12) will be denoted by a two-index tensor according to [63] ( 1 , m 1 ) with the degeneracy part being and the structural part imposing the conservation of the angular momentum, 1 = 2 , see Eq. (13). The evaluation of the partition function in Eq. (11) is now straightforward, since the angular momenta have to be the same across all links of the one-dimensional chain. Therefore, it reduces to a sum of powers of Bessel functions as it is also known from its exact solution (specifically, the model is Bethe integrable) [71,72]. Here, N denotes the number of edges of the chain. The factor (2 +1) is due to the degeneracy of the angular momentum eigenvalues . Considering a chain with open boundary conditions in the thermodynamic limit, only the largest eigenvalue = 0 contributes to Eq. (16) and the partition function becomes which is equivalent to Fisher's expression in Ref. [71]. In Section 3.4.1, we will see what this expression does tell us about thermodynamic quantities. Before doing so, we first proceed with the tensor network notation for the targeted square lattice.

Application to the two-dimensional square lattice
Making use of the graphical notation, the partition function for the two-dimensional Heisenberg model on a square lattice is given by . (18) Similarly to the construction of the one-dimensional chain, the decomposition of the Boltzmann factors allows us to collect all terms into local objects around the vertices k ∈ V . This time though, the integral over the two continuous parameters involves four spherical harmonics corresponding to the four links attached. Again, it is non-vanishing if and only if the total angular momentum is zero: This could be ensured by fusing the two incoming angular momenta ( 1 and 2 ) to an intermediate value , which is then split again into the two outgoing ones ( 3 and 4 ). The coupling rules of spherical harmonics are in close connection to the spin coupling rules of SU (2), which makes an SU (2)-symmetric tensor network ansatz again very convenient. The precise relation reads where the numerical factors converting between the two representations are given by and the fusion tree represents the Clebsch-Gordan coefficients. Besides the usual constraint m 1 + m 2 = m 3 + m 4 , we have here an additional one dictated by the reflection symmetry of the O(3) group, namely Ultimately, the partition function for the Heisenberg model on the two-dimensional square lattice can be constructed by the contraction of local four-index tensors of the form .
The degeneracy tensors encompass this time not only the numerical Boltzmann factors as in the onedimensional case, but also the conversion factors of Eq. (22) and are hence given by The handling of SU (2)-symmetric tensors is implemented in an automated library, that deals with the symmetry constraints efficiently in terms of fusion trees [63]. Some further remarks are in order. The representation of the partition function as a tensor network will necessarily be an approximation. This is due to the fact that the plane wave expansion of the Boltzmann terms in Eq. (8) is only exact in the limit → ∞. In practice, the links in our tensor network ansatz will always deal with a finite number of angular momenta 0, 1, . . . , j max (we will use for angular momenta and j for SU (2) quantum numbers in our TN). This in turn means that the number of representations kept will determine the accuracy of our simulation and provides a systematic refinement parameter, as we will demonstrate in the next sections.

Partition function with corner transfer matrices
Let us now demonstrate how the tensor network representation is used in practice to compute the partition function. The central object for the Heisenberg model on the square lattice is the four-index tensor of Eq. (23). The full partition function is then the contraction of the infinite square lattice with a four-index tensor on every vertex. Although the construction is exact in principle, practical reasons of the tensor network ansatz force us to approximate the calculations of Z. These approximations are controlled by the number of angular momenta used in the plane-wave expansion of the Boltzmann factors, see Eq. (8). The number of angular momenta kept is then directly related to the bond dimension of the bulk tensors. All the links of the tensor network are described in terms of SU (2) quantum numbers denoted by j t j , where j is the spin in the language of Clebsch-Gordan coefficients (or angular momentum in the language of spherical harmonics, with the above mentioned further restrictions to fusion rules) and t j its degeneracy (see also Ref. [63] for details about SU (2) symmetry in TNs). Keeping only the lowest value of j B = 0 1 amounts to a bond dimension of χ B = 1, i.e., a product state. Naturally, as we increase the bond dimension by keeping more spin quantum numbers, Figure 1: Decay of the Bessel functions I +1/2 (β) normalized to I 1/2 (β) as a function of the angular momentum , for several values of the inverse temperature β = 1/T . the accuracy of the approximation of Z increases as well. Here, χ B counts the number of terms in the plane-wave expansion of the Boltzmann factors, such that a tensor network with bond dimension χ B has spins j B = [0 1 , 1 1 , . . . , (χ B − 1) 1 ] on the links. Notice also that χ B is the symmetric bond dimension, while the effective bond dimension is larger due to the intrinsic degeneracy of SU (2) quantum numbers (a spin j is (2j + 1)-fold degenerate, see also Table 1). In our study we consider the temperature range of T ∈ [0.01, 100]. For each individual temperature we initialize an infinite tensor network with a single-site unit cell tensor according to Eq. (23) for bond dimensions up to χ B = 6, i.e., with a maximal spin of j max B = 5 on the links. It is expected, that this bond dimension is more than sufficient for large temperatures. However, at low temperatures the approximations due to a finite χ B will become strong, as forecasted by the behaviour of the Bessel functions in Fig. 1. While the Bessel functions drop off very quickly for large T and χ B = 4 already yields a relative accuracy of ≈ 10 −5 , the bond dimension would have to be χ B ∼ 50 to achieve the same level of accuracy at low T . The rather small chosen value of χ B = 6 is due to limited computational power, especially when actually contracting the infinite TN, as we will elaborate next. The contraction of the infinite two-dimensional square lattice can only be computed approximately, for instance using a standard corner transfer matrix (CTM) procedure. This is basically a power method in which the tensor for the partition function is absorbed into so-called environment tensors, until they converge to an environment fixed-point. The contraction of the infinite lattice is then approximated by eight fixed-point tensors for a translationally invariant system with a single-site unit cell. Arrows are omitted here for better clarity. Since the indices of the bulk tensors (the white tensors in Eq. (25)) carry only integer spin representations, the CTM environment bond indices can generally carry integer, half-integer or mixed, i.e., integer and half-integer spin representations. The natural choice of only integer spin irreducible representations worked out to be the best one, with smooth convergence of the environment tensors. We notice here that this could be interpreted as the full O(3) symmetry being preserved, without room for spontaneous breaking of the Z 2 subgroup: The latter would instead result in half-integer, SU (2)like irreducible representations in the environment tensors. The CTM procedure necessarily introduces a second approximation in the calculation of the partition function. It stems from the contraction of the two-dimension tensor network, which can not be performed exactly without an exponential growth of the simulation time. The approximations are controlled by the environment bond dimension χ eff E , that is the bond dimension along the perimeter of the CTM tensors in Eq. (25). For each fixed χ B , the environment bond dimension should be increased up to the point, where the environment tensors are converged in χ eff E . In this case the approximations of the contraction of the infinite two-dimensional lattice are insignificant and the only approximation is the one of the bulk tensors, χ B . However, convergence in χ eff E is only achieved for rather small χ B , while approximations in χ eff E remain for larger χ B due to limited computational power. In our simulations we work with different environment bond dimensions between χ eff E = 100 and χ eff E = 1500. In contrast to the bulk bond dimension χ B , these are now effective bond dimensions. The spin sectors on the CTM environment tensors are fixed to be integer but the truncation procedure automatically adapts to the most relevant sectors. As an example, we lists the actual quantum numbers for the SU (2)-symmetric CTM tensors at T = 0.01 in Table 1, where the (effective) bond dimension of the environment tensors χ eff E is broken down into its symmetric bond dimension χ E and quantum numbers 2 . As expected, for larger χ B and hence higher j max B , the algorithm also keeps more spin sectors on the CTM bond indices. 2 Notice that χ eff E refers to the theoretical maximal bond dimension. Due to the preservation of SU (2) multiplets in the tensor network truncation procedures, the actual bond dimension can be slightly less in practice.

Thermodynamic quantities from the partition function
The partition function can be used to compute all the key properties of a thermodynamic system. In particular, it gives access to the free energy F , the internal energy U and the thermodynamic entropy S. In this section we will compute those quantities directly from the tensor network representations described in Section 3.1 and Section 3.2. Notice that we employ a notation where the quantities are always defined per lattice site, so that we avoid obvious divergences in the thermodynamic limit.

One-dimensional linear chain
In one spatial dimension, the partition function for the Heisenberg model has been shown to be where the SU (2)-symmetric matrices on each lattice site are constructed according to Eq. (14). Using the tensor network construction we can directly compute the free energy per site where λ(β) = sinh(β)/β is the ( = 0) leading eigenvalue of the transfer matrix. This matches the exact solution for the free energy, since there are no approximations present in the one-dimensional construction (see Section 3.1). One quantity we can directly compute from F (β) is the thermodynamic entropy per site, given by S(β) = −∂F/∂T . The free energy F and the thermodynamic entropy S are shown in Fig. 2. It is known that thermodynamic potentials are concave functions of their intensive variables and convex functions of their extensive variables [73]. This is indeed found for the free energy F , for which the logarithmic temperature scale skews the picture. However, due to the permanent positive slope, the thermodynamic entropy becomes negative. This is somewhat surprising and quite counter-intuitive, and the entropy even diverges as T approaches zero like S(T → 0) → −∞. This behaviour is however typical for classical spin models [71,73] with continuous degrees of freedom 3 , and not a misfeature of the tensor network construction. Now that we have validated our tensor network construction against an exactly solvable model, we will move our focus to the more interesting two-dimensional square lattice scenario.

Two-dimensional square lattice
Similarly to the one-dimensional case, we define the partition function per site as We can compute it again from the dominant eigenvalue of a transfer operator (an infinite row or column of four-index tensors), and then express the latter itself in terms of the dominant eigenvalue of its constituent transfer matrix in the other direction (column or row, respectively). Without going into the details, z(β) can thus be expressed in terms of the CTM tensors in such a way that their normalization fully cancels, see for instance Ref. [74]: We can now proceed to extract the free energy F and the entropy S per site. In both panels of Fig. 3 we show the data for all available bulk bond dimensions at the maximum environment bond dimension χ eff E = 1500. The overall behaviour of F and S resembles the one of the one-dimensional Heisenberg chain, and the negative entropy is again not a misfeature of our tensor network method. However, in the two-dimensional setting we are dealing with approximations in z(β) and hence in the thermodynamic quantities. Based on the truncation of the bulk bond dimension χ B and the behaviour of the Bessel functions for large β (refer to Section 3.3), we identify low temperatures, i.e., T ≤ 0.2, as the regime in which our approximations become strong, and we highlight this with a grey background. In this regime a much larger bond dimension χ B would be desirable for all links in the TN. It might well be that the approximations lead to larger error bars, while the data points could still exhibit a good trend towards their actual values. However, the concrete impact of the χ B -truncation is difficult to assess. The choice of the boundary value of T = 0.2 is supported by the spread of data curves in Fig. 3, and even clearer confirmed by further analysis in subsequent sections. A lighter gray highlighted area for 0.2 < T ≤ 0.5 is introduced to indicate the region where our tensor network data seem to be not dramatically spread with respect to χ B , but an extrapolation in the environment bond dimension χ E looks crucial to fully capture the physical behaviour (see Section 4).

Entanglement analysis
One of the many advantages of using tensor network-based approaches is that we can comparably easily get access to entanglement properties almost directly by construction. These quantities can then be used to describe and detect phases as well as the phase transitions between them. For instance, in one spatial dimension, the entanglement spectrum, i.e., (minus log of) the eigenvalue spectrum of a reduced density matrix, can be extracted from a bipartition of a matrix product state (MPS). This can then be analysed to spot signatures of distinct topological and trivial phases of matter protected by symmetries [75][76][77][78][79][80]. It can also be used to compute the entanglement entropy, as well as Rényi entropies, whose scaling behaviour can then be used to locate phase transitions. Similarly, in two dimensions, the environment tensors obtained from the CTM procedure after contracting the infinite two-dimensional lattice can also be used to compute different quantities that are associated to entanglement properties of some underlying quantum model. For the case at hand in this paper, this allows us to determine some characteristics of the model in different temperature regimes. To be more specific, the model that we study here is two-dimensional classical, but the associated "entanglement" computed from the CTM tensors is that of the associated dual one-dimensional quantum model in the same universality class [81]. Therefore, even though the model we study here is purely classical, these quantum-inspired quantities, nonetheless, give us valuable information also about the classical correlations present in our system, which can then be used to characterize phases and their transitions. We describe them in more detail below.

Corner entropy
The environment fixed-point tensors can be directly used to compute the corner entropy S C , a quantity that is known to contain information about the universal properties of the model at hand [81]. As such it can help to pinpoint phase transitions. The corner entropy is related to the entanglement entropy of a bipartition for some infinite one-dimensional quantum Heisenberg model, due to the general mapping between D-dimensional classical statistical mechanics models and (d + 1)-dimensional quantum models [21]. It can be computed from an eigenvalue decomposition of the matrix C := C 1 ·C 2 ·C 3 ·C 4 , whose graphical notation in the tensor network language is The eigenvalues λ α appear with different degeneracies due to the block-diagonal structure of the SU (2)-symmetric matrix C. The corner entropy can then be computed using where the degeneracies of the eigenvalues λ α in each spin sector need to be taken into account. Results for all available bond dimensions up to χ B = 6 are shown in Fig. 4. For temperatures T ≥ 0.2 the corner entropy is essentially converged in the virtual spin irreducible representations, however, in the low-temperature regime for T < 0.2 more quantum numbers would be needed to achieve convergence, as indicated by the separation of the data points and expected from the properties of the Bessel functions (cf. Fig. 1). The behaviour of the corner entropy is a first indication of different physical properties in the low-and high-temperature regime. While it hints at an uncorrelated phase at high temperature, the increasing correlations towards low temperature already indicate the emergence of a highly-correlated low-temperature phase. To be precise, in the figure we can distinguish three different regimes: above T ∼ 0.6 (vanishing entropy), between T ∼ 0.6 and T ∼ 0.2 (linear behaviour in T and converged in χ B ), and below T ∼ 0.2 (weak dependence on T and S C increasing with χ B ).

Geometric entanglement
As we have mentioned in the introduction (Section 1), some past works have hinted at the possible presence of a BKT-like transition at low temperatures -incidentally compatible with the change in the slope of the corner entropy in Fig. 4. Such kind of transitions are however challenging to detect numerically: Most of the usual figures of merit fail to capture them in a clear-cut way [82]. In order to confirm or disqualify such a transition, we resort here to the geometric entanglement (GE) for the boundary MPS obtained from the CTM procedure. It is indeed known that the global geometric entanglement is able to capture elusive phase transitions [82] even in the case where no indications are found in other entanglement measures such as the entanglement entropy [83]. The geometric entanglement describes the proximity of an entangled quantum state vector |ψ to the closest product state vector |φ via the overlap The geometric entanglement per site is then defined as We computed the geometric entanglement per site for the boundary MPS obtained after contracting a half-infinite two-dimensional lattice of the partition function tensors. As such, this MPS corresponds to the ground state of the associated dual one-dimensional quantum model (with entanglement entropies such as the corner entropy computed in the previous section). This MPS can be easily built from our CTM technique, just by considering an infinite MPS with tensor T 1 at every site. In Fig. 5 we We observe a smooth crossover, which is incompatible with a phase transition, even of the BKT type, as explained in Ref. [82].
show the geometric entanglement per site computed for such boundary MPSs. For large temperatures the proximity to a product state is almost perfect, which yields a vanishing GE. For lower temperatures it becomes increasingly difficult to approximate the boundary MPS with a product state, as already expected from the behaviour of the corner entropy. The GE shows however a crossover but not any sign of a phase transition. As shown in Ref. [82], the geometric entanglement per site shows a clear peak when crossing a BKT transition, being in fact very sharp and similar to a level crossing. This is due to the fact that it is a multipartite measure of entanglement, and therefore can capture the subtleties of global changes in the structure of the physical system, as opposed to other entanglement and correlation measures that are essentially more local. In our case, as we have discussed above, we see no peak whatsoever, and just a smooth crossover.

Computation of expectation values
In this section, we will show how to compute expectation values of observables, both local and at a distance, and discuss how to extract the model's correlation length, both from fitting correlations as well as directly from the tensor network formalism. The aim is to explore the full temperature regime and try to shed light onto one of the most important core issues of the model, namely that of a finitetemperature phase transition. In particular, we will apply a recently introduced scaling hypothesis for finite bond dimension data in order to better estimate the true correlation length: This will help avoiding to over-interpret the bare data at low temperatures, where numerical approximations intrinsic to the method might cap the (very large) correlation length. Nonetheless, at first sight, the data seem to fall well onto the typical fitting form of BKT-like transitions at the before mentioned finite temperature T c 0.5, although no other signatures thereof are found (see also Section 3.5.2). Conversely, we show that our data are overlapping with reference ones in the lattice gauge theory community (in the region where both are available): By performing the same kind of scaling to infinite correlation length customary in that context, we assess the consistency of our data with the asymptotic freedom hypothesis, i.e., with no finite-T transition. The puzzle thus still remains open, but we foresee that massive tensor network simulations in the regime T ∈ [0.1, 0.5] would considerably help to lift the mystery in the near future. We plan to take up this formidable challenge in upcoming works.

Construction of observables
The computation of expectation values is relatively straightforward in our present formulation of the classical partition function. An n-site observable results in the contraction of a tensor network with n local tensors modified in a non-trivial way, dictated by the decomposition of the observable in terms of the group characters, plus possibly a few more trivial tensors mediating the correlation without changing their own character. Let us see in practice what this means for certain relevant cases up to n = 3. In general, an n-site observable is given by whereŜ 1 , . . . ,Ŝ n act on n different lattice sites and not necessarily on site 1, . . . , n. Naturally, the infinite contraction of the tensor network surrounding the n-sites with non-trivial modifications can again be approximated by the fixed-point environment tensors. The only non-trivial one-site observable is in principle O(Ŝ k ) =Ŝ k , i.e., the magnetization per site. However, since this operator is not SU (2)-(resp. O(3)-) invariant, the local magnetization has to vanish. This is not only expected for a manifestly invariant Ansatz like ours, but it is also anyway forbidden by the Mermin-Wagner theorem, which excludes spontaneous breaking of the continuous symmetry of the Heisenberg model.
As a non-trivial two-site observable we will consider the spin scalar product O(Ŝ i ,Ŝ j ) =Ŝ i ·Ŝ j . When the operator acts on nearest-neighbours, it measures the (negative) bond energy, otherwise it measures generic spin-spin correlations at a distance |i − j|. The spin scalar product can be written in terms of spherical harmonics acting on sites i and j aŝ This additional spherical harmonic introduces a fifth tensor index carrying a spin-1, which then calls for straightforward modifications in the underlying fusion trees and, consequently, in the conversion factors to Clebsch-Gordan coefficients. Assuming, without loss of generality, to measure the spin-spin correlations along horizontal bonds, the whole expression reads ( 1 ) The two operators acting on vertical bonds are defined in complete analogy, and such constructions can be straightforwardly extended to three-sites observables, like for instance The central tensor (the bottom right one in this example) has two additional spin-1 channels, which connect to the scalar product operators on both sites respectively. The underlying fusion tree determines which expectation value gets measured, Setting the eminent internal spin to j = 0 yields an identity on the central site, and therefore nextto-nearest neighbour spin correlations Ŝ i · I j ·Ŝ k ; replicating this an arbitrary number of times can therefore be used to generate spin scalar product at larger distances. Choosing j = 1 instead generates the interesting triple product O(Ŝ i ,Ŝ j ,Ŝ k ) = Ŝ i · (Ŝ j ×Ŝ k ) . Measuring this observable on all triangles on a square lattice plaquette could be useful in uncovering twisted spin configurations, similar to those that appear in magnetic Skyrmion systems. Noticeably, this quantity is SU (2)-invariant but breaks the additional discrete Z 2 reflection symmetry of the O(3) group. Therefore, it could in principle exhibit finite expectation values without violating the Mermin-Wagner theorem.

Bond energy
We now focus on the bond energy, which is computed using the scalar product operators introduced previously in Eqs. (36)- (38) and setting |i − j| = 1 as Here the square tensors are the ones generating the scalar product, while all the other ones are the analytic tensor for the partition function as in Eq. (23). Notice the additional negative sign that appears due to the ferromagnetic coupling in the Hamiltonian. Using the fixed-point environment tensors computed with the CTM procedure, the evaluation of the expectation value simplifies to the evaluation of the two rightmost tensor network diagrams. The results are shown in Fig. 6 for the accessible bond dimensions. Due to the spatial isotropy of the Hamiltonian and the tensor network construction, the energy along horizontal and vertical links is identical. The behaviour of Ŝ i ·Ŝ j for a fixed χ B hardly Figure 6: Bond energy − Ŝ i ·Ŝ j on horizontal and vertical nearest-neighbour links. The inset shows the convergence of the bond energy vs. the inverse bulk bond dimension at T = 0.01. As T → 0, the relative weights of all angular momenta in Eq. (24) tend to become equal, so that the only sensible scaling indeed seems to be against the (inverse) number of the allowed ones.
depends on the environment bond dimension, so that the results are converged in χ eff E . This is due to the very local support of the operators and not the case for general operators or global properties.
As expected for large temperatures, the spins are fully disordered and the bond energy approaches 0 for T → ∞. In contrast, for low temperatures a tendency towards ferromagnetic alignment develops and results in a bond energy tending to −1 for T → 0. The data curves for different χ B nicely overlap down to temperatures of T ≈ 0.2. As we have argued before, at this point our approximations due to limited bond dimension become sizeable. However, the expected behaviour is recovered by increasingly larger values of the virtual spin j max , i.e., system bond dimension χ B .

Spin-spin correlations: The tensor network data and a finite-T c hypothesis
The bond energy in Section 4.2 can be seen as being a specific case of general spin-spin correlations Ŝ i ·Ŝ j at a distance |i − j| for lattice sites i and j Here the leftmost and rightmost tensors are the ones for the scalar product from Eqs. (36)- (38), while the ones in between act as identities and simply loop through the interaction (this is achieved by setting j = 0 in the fusion tree, see Eq. (40) and surrounding text). It is also important to stress that the denominator must be computed via the full contraction over the same finite distance |i − j| as the numerator, in order for all eigenvalues of the transfer matrix to be accounted in the same way. As expected for a ferromagnetic model, the correlations have matching signs. For large temperatures, the results nicely match an exponential decay with a short correlation length ξ, which can therefore be quite comfortably extracted by standard fitting procedures. For low temperatures, instead, it becomes increasingly difficult to distinguish the correlation behaviour from a quasi-algebraic decay. Very large separations between the spins would be needed due to a rapid increase of ξ, but these are computationally expensive to be performed according to Eq. (42). The tensor network representation of the partition function provides however a direct and much cleaner access to compute the correlation length, i.e. via the eigenvalue spectrum λ α of the transfer matrix T of the boundary MPS. The transfer matrix is constructed out of the infinite half-column CTM tensors T 1 and T 3 (or, equivalently, from T 2 and T 4 in the other direction) according to The correlation length is then computed as [30, 31] where λ 1 and λ 2 are the two largest eigenvalues, and the rightmost expression corresponds to a normalized spectrum such that λ 1 = 1. The additional advantage of this method is that the extracted ξ is automatically the largest one in the system, irrespective of which correlation function gives rise to it. In Fig. 7 we show the correlation length ξ for different values of χ B at the largest available environment bond dimension of χ eff E = 1500. The inset in Fig. 7 shows a good agreement of our data with Monte Carlo simulations [84], which are among the reference ones for the lattice gauge theory community. We will see later how to reconcile them with the predictions by asymptotic freedom, as it is customary in that context. Here, instead, we first point our attention to the apparently striking consistency of the data with a fit of the kind with T c = 0.509 ± 0.003 and c 1, which is commonly taken as a pristine signature of a Berezinskii-Kosterlitz-Thouless transition [16,17]. We try to heal finite bond dimension effects by adopting a scaling procedure recently introduced to cope with finite-entanglement in quantum tensor network simulations [40,41,85], as described in the following paragraph and displayed in Fig. 8. Thereby, the correlation length can be confidently extrapolated down to a bit below T = 0.5, which is shown in the main panel of Fig. 7 and the actual values are reported in Table 2. This indicates, most probably, that the BKT-like fit of Eq. (45) is not really reliable, though intriguing at first sight away from the putative T c . It is however not excluded, from the numerical data at disposal, that a divergence sets in at lower temperatures.
Our tensor network includes two approximation parameters, χ B and χ eff E . However, for every fixed χ B the second approximation in the CTM environment tensors could in principle be eliminated for Figure 7: The largest correlation length in the system as extracted from Eq. (44) for χ B = 2 to χ B = 6 at the largest available environment bond dimension χ eff E = 1500 (usual color legend). Additionally, we show the correlation length extrapolated to χ eff E → ∞. Importantly, there is no drift of the temperature dependence in the bond dimension χ B . The excellent agreement with the exponential fit of Eq. (45), down to T = 0.6 at least, seems at first sight to indicate a BKT-like transition at a finite temperature T c ≈ 0.51. The choice of the exponential fit is supported by the topright inset, which shows the correlation length rescaled to √ T − T c ln ξ(t) vs. √ T − T c . However, even the extrapolated correlation length starts to deviate from a BKT-like divergence for T ∼ T c . The bottom-left inset shows excellent overlap with the Monte Carlo data presented in Ref. [84], which are usually employed to claim agreement with the asymptotic freedom description of the underlying field theory. For a thorough discussion, see the main text.
χ eff E → ∞, which provides us with a reliable scaling. As has been demonstrated in Refs. [41,85], the eigenvalue spectrum of the transfer matrix T can be used to define a scaling parameter δ = δ(χ eff E ). The true transfer matrix eigenvalue spectrum would indeed be continuous above a possible initial gap, which in turn determines the correlation length of the system, as described above in Eq. (44). However, the necessarily finite bond dimension χ eff E naturally introduces a cutoff and therefore a maximal length scale. The parameter δ then describes the discreteness of the spectrum, i.e., its distance to a continuous spectrum. While the first gap is related to the correlation length, the second gap of T is already suitable to define a scaling using A more sophisticated scaling parameter can be computed with a linear combination of the first n eigenvalues, with coefficients that can be determined manually or even in an optimization procedure [41]. For our purposes, the basic definition in Eq. (46) works equivalently well. In order to perform an infinite-bond dimension extrapolation for the correlation length, we utilize a relation be-   Fig. 8. This procedure becomes unreliable for the values highlighted in gray, as discussed in the main text. Monte Carlo data from Ref. [84] is shown for comparison. tween = 1/ξ = − log |λ 2 | and δ according to presented in Ref. [85]. In Fig. 8 we depict the finite-entanglement scaling for χ B = 6 at different inverse temperatures β, with δ computed for all available χ eff E . For small values of β the scaling procedure yields a reliable extrapolation of 1/ξ, as demonstrated in the lower panel for β = 1.0 and β = 2.0. Notice that a change of slope sets in at values of δ < 0.02 in the central lower panel of Fig. 8. We cannot exclude that this will happen at even lower values of δ for lower temperatures (larger β), i.e., that the estimate of ξ might be strongly affected by data with even larger bond dimensions χ eff E . Therefore, we highlight extrapolated values for β ≥ 2.2 in Table 2. Nevertheless, the scaling in the confident temperature range demonstrates a rapidly increasing correlation length, which however does not automatically imply a true divergence. It might as well be that a saturation sets in for large χ eff E that are currently inaccessible and therefore point at no finite-T transition, but rather to a quasi-critical region at low temperatures. The absence of a finite-temperature phase transition would be consistent with the geometric entanglement analysis performed in Section 3.5.2, that seemed not to show any sign of a BKT-like transition despite its sensitivity to it.

Spin-spin correlations: A consistency check with asymptotic freedom
As we have recalled in Section 2, the plausibly most recognised theoretical framework to describe the classical Heisenberg model in two spatial dimensions is the one relating it to the O(3) non-linear sigma model in 1 + 1 dimensions and its asymptotic freedom. It would however be pointless to try demonstrating such mechanism from the behaviour of the correlation length itself, since no divergence of the kind ξ(T ) exp(−1/T ) will ever be visible in Fig. 7 within the achievable regimes. The lattice gauge theory community has instead developed a quite sophisticated consistency check with such a picture in terms of the measured spin-spin correlations [86][87][88][89], which we now illustrate and perform on our data for a maximal bond dimension of χ max E = 1000. For this sake, we should first extrapolate physical quantities to the limit of zero lattice spacing, a → 0, or, equivalently, to the limit of a diverging correlation length, ξ → ∞, while keeping the lattice spacing fixed as in the typical condensed matter context adopted here. This means that we should look at correlation functions at a fixed distance in terms of ξ, i.e. C(m, ξ) := Z(ξ) S(x = mξ) · S(0) .
We omit the vector for the position x here, since we will measure correlations along horizontal bonds in the square lattice tensor network. Z(ξ) is a constant with respect to m, to be chosen such that the limit Figure 9: Collective analysis for the quantity C(m, n = 1) at χ B = 2 and χ B = 6, for different values of m. A light to dark shading has been used for increasing environment bond dimensions χ eff E . The fit to extract the infinite correlation length limit has been performed for χ eff E = 1000 using Eq. (51). The legend applies to both panels.
is finite. Such a (renormalization) constant is guaranteed to exist, but it is quite intricate to compute it carefully, since Z diverges itself in the proximity of critical points. While this is usually happening in a power-law fashion, in strict relation to what we usually call critical exponents (or anomalous dimensions), in the context of asymptotic freedom the divergence of Z is itself logarithmically slow (a.k.a., the associated operator is RG-marginal). Therefore, a much cleaner approach is to define instead the ratio of such correlation values C(m, n) := lim ξ→∞ C(m, ξ)/C(n, ξ) = C(m)/C(n) , thus circumventing the need to determine Z(ξ) explicitly for different values of ξ. We will set n = 1 as the reference distance in the following.
In principle, one could first extrapolate the raw correlations to the limit of infinite bond dimension (χ B , χ eff E → ∞) for each value of the coupling β, as we did in Section 4.3, and then define C(m, n, ξ) to perform the ξ → ∞ limit of Eq. (50). Alternatively, we could also consider first the raw correlation data C(m, n, ξ(χ B , χ eff E )) as a function of ξ(χ B , χ eff E ), and then extrapolate to the limit of ξ → ∞ -most practically, first for each χ B separately. Both approaches lead to fully compatible results, though the second one turns out to be numerically more stable. In Fig. 9 we display the second approach for β ≤ 5, which extends beyond the window of overlap with the mentioned Monte Carlo data of Ref. [84]. While the behaviour fluctuates strongly for small β (which correspond to small ξ), it smoothens out for larger values, and an extrapolation is possible. Here we employ a simple fit of the form Additional logarithmic corrections of the form poly(1/ log(1/ξ)) predicted by the Symanzik effective theory of lattice artefacts [88,89] should also be included, if one is interested in the precise quantitative value of the limit as the QCD community is [86,87]. This however goes beyond the scopes of this first tensor network analysis, considering the currently accessible bond dimensions. In Fig. 10 and Table 3 we show the extrapolated physical values C(m, n = 1, 1/ξ → 0) for different bond dimensions χ B of the partition function approximation. While χ B = 2 is a too rough approximation, allowing for a limited range of 1/ξ only, the differences for χ B = 3 up to χ B = 6 are hardly visible at this stage. This is however expected to change as soon as larger β values are introduced in the analysis, and quantitative corrections to the extrapolated limits could be foreseen. They are however not forecasted to substantially change the behaviour in the logarithmic scale of such plot, as much as the error bars also do not (see Table 3), so we can safely leave them for future investigations. Most importantly, the behaviour of the (physical) correlation function at large distances is predicted by the so-called spectral representation with E(p) = p 2 + 1, which roughly gives C(m, 1) exp(−m)/ √ m. In Fig. 10, we display the result of a single-parameter fit to adjust the multiplicative constant A to our data points for m ≥ 10, which results in A = 63.19 ± 0.06. The agreement is evidently already very good, and could possibly improve with future data sets for larger (χ B , χ eff E ), though this could be quite a challenging task.
Summarising, we seem to obtain consistency of the data computed via our tensor network approach with the predictions of asymptotic freedom, though strictly speaking we are not providing evidence that this limit ξ → ∞ does take place indeed at T → 0 and not earlier. This leaves the door open for alternative scenarios as, e.g., the lattice model having its critical coupling at a value larger than zero, as suggested by some mathematical physicists [90,91].

Conclusions and outlook
In this work, we have presented a fresh approach for the study of the two-dimensional classical Heisenberg model based on tensor networks. The partition function for the statistical mechanics model is obtained by the contraction of a network of interconnected tensors residing on the lattice sites. The local tensors are constructed from an analytic rewriting of the Boltzmann factors, whose precision is controlled by a cutoff in the infinite sum over the rotational group characters (i.e., angular momenta). The correct fusion of angular momenta at each lattice vertex is automatically ensured by directly imposing the O(3) symmetry at the level of the constituent tensors.
Our technique gives direct access to thermodynamic quantities and physical observables for a formally infinite system over a wide (classical) temperature range. Noticeably, we can directly access a wealth of information about correlations in the system via the transfer matrix spectrum in our formalism, without the need for direct measures and/or for exponential fittings at prohibitively large distances. The finite bond dimensions employed throughout the algorithm play the effective role of spatial and thermal cutoffs. Therefore, data have to be extrapolated accordingly, which we did following protocols recently introduced in the (quantum) tensor network context. Moreover, we explored the behaviour of different entanglement measures, leveraging on the classical to quantum correspondence in different dimensionality.
Our findings bring new clues in the ever-lasting riddle about this fundamental model: On the one hand, the spin-spin correlation functions, once processed as commonly done in the lattice gauge theory community, display a satisfying agreement with the vastly endorsed mechanism of asymptotic freedom. This predicts a very rapidly diverging correlation length at low-temperatures but no finite temperature transition of any kind, i.e., the spin system remaining disordered down to zero temperature and the corresponding O(3) non-linear sigma model remaining gapped down to zero coupling. The so-called physical value of the correlations, i.e., the one extrapolated in the zero lattice spacingor, equivalently, infinite correlation length -limit, indeed matches the behaviour forecasted via the spectral representation.
On the other hand, the (largest) correlation length itself shows at first a striking agreement with a fitting form usually regarded as a bona fide signature of a Berezinskii-Kosterlitz-Thouless transition at a finite temperature T c 0.5J. We incidentally stress here that the above mentioned extrapolation to the physical limit -as it is performed -does not strictly imply that this is reached at zero temperature. Our numerical results, however, suffer from severe approximations in the very-low temperature limit and important deviations appear close to such a putative T c . We also cannot exclude the onset of some quasi-critical behaviour at very low temperatures. This seems in turn to be further supported by the entanglement-based quantities we have been looking at, for which we observe a strong entanglementscaling dependence in the low-temperature regime.
We could envision future (larger-effort) investigations to compute further quantities like Noethercurrents correlations, form factors as well as others and compare their behaviour with the exactly computable ones for the integrable O(3)-NLSM, in order to let the asymptotic freedom mechanism undergo more stringent tests. Conversely, the availability of reliable data for lower and lower temperatures shall also contribute to lift the mystery about the apparent finite-T BKT-transition and/or the quasi-critical behaviour. Furthermore, even more ingenious extrapolation schemes from the achievable finite bond dimension data would be of the utmost help in elucidating the subtle details of this intriguing problem. An example could consist in quantifying the intrinsic cutoff of the characters' expansion with a single parameter, and then searching for some sort of regression plane to achieve the simultaneous limit of both truncations of our ansatz.
Put into a broader perspective, the method that we have employed here offers a strong potential for the detailed and systematic exploration of a wealth of classical and quantum spin models. It will be interesting to flesh out its applications for physical systems defined on frustrated geometries and for real materials where tensor networks have shown great promises recently [47,[92][93][94][95][96][97]. These efforts will help us in gaining further insights into exotic phases of matter such as quantum spin liquids including their finite temperature properties and quantum to classical crossovers, a topic of intense Funding information. We acknowledge DFG funding through CRC 183 (B01), for which this is an inter-node Berlin-Cologne project of the CRC, EI 519/15-1, GZ OR 381/3-1 as well as GZ SCHM 2511/10-1. This work has also received funding from the Cluster of Excellence MATH+ and the European Union's Horizon 2020 research and innovation programme under grant agreement No. 817482 (PASQuanS). R. O. also acknowledges support from Ikerbasque and the DIPC.

A Numerical data for the consistency check with asymptotic freedom
For completeness and reproducibility we show the final infinite correlation length limit C(m, n) = lim ξ→∞ C(m, ξ)/C(n, ξ) used in the analysis of Section 4.4 in Table 3. Table 3: Numerical data for the infinite correlation length limit of C(m, n) = lim ξ→∞ C(m, ξ)/C(n, ξ), for all available bulk bond dimensions χ B at different values of m. As in the corresponding Fig. 10, the reference is set to n = 1.