Resolving the Berezinskii-Kosterlitz-Thouless transition in the 2D XY model with tensor-network based level spectroscopy

Berezinskii-Kosterlitz-Thouless transition of the classical XY model is re-investigated, combining the Tensor Network Renormalization (TNR) and the Level Spectroscopy method based on the finite-size scaling of the Conformal Field Theory. By systematically analyzing the spectrum of the transfer matrix of the systems of various moderate sizes which can be accurately handled with a finite bond dimension, we determine the critical point removing the logarithmic corrections. This improves the accuracy by an order of magnitude over previous studies including those utilizing TNR. Our analysis also gives a visualization of the celebrated Kosterlitz Renormalization Group flow based on the numerical data.


I. INTRODUCTION
The Berezinskii-Kosterlitz-Thouless (BKT) transition was historically the first example of topological phase transitions, which is now an essential concept in physics [1]. The canonical model exhibiting the BKT transition is the classical XY model in two dimensions, which is defined by the total energy (classical Hamiltonian) where θ takes the angular value (0 ≤ θ < 2π) at each site on the square lattice and only the nearest-neighbor interactions are considered. The topology of the configuration space allows topological point defects (vortices and antivortices), whose dissociation drives the BKT transition. The BKT transition can be described in terms of two running coupling constants: y K controlling the thermal fluctuation of spin waves, and y V representing the vortex fugacity (see Eq. (A3) for the precise definitions). The nontrivial interplay between these two couplings is described by the celebrated Kosterlitz renormalization group (RG) equation [2] where l ∼ log L is the logarithm of the length scale L [3]. The phase diagram is Fig.1 and the transition line becomes y K = y V . We can introduce new variables g and t by y V = g + t and y K = g − t, so that the phase boundary corresponds to t = 0. At the BKT transition, t = 0, the RG equation for g is reduced to dg/dl = −g 2 , which implies g ∼ 1/l ∼ 1/ log L. This slow decay is the source of the notorious logarithmic corrections.
The BKT transition is conceptually well understood in terms of the Kosterlitz RG equation. However, the FIG. 1. Phase diagram of the BKT transition of the XY model. The behavior of the low-temperature phase is controlled by the c = 1 critical line originated from U(1) symmetry. The RG flow separates the phase with yV → 0 from yV → ∞ phase when l increases famous "Kosterlitz RG flow" has remained a rather abstract concept which cannot be seen directly. Moreover, because of the logarithmic corrections, significant finitesize effects persist even in a large system, making conventional finite-size scaling of Monte-Carlo methods such as Binder plot [4,5] powerless. Even with considerable efforts over decades, an accurate determination of the critical temperature remains difficult even with a huge computational power.
On the other hand, many 1D quantum systems can be also described by the same effective theory and thus also exhibit the BKT transition. Interestingly, a powerful numerical finite-size scaling method called "Level Spectrscopy" was developed specifically for those 1D quantum systems by Okamoto and Nomura [6][7][8][9][10]. Based on the Conformal Field Theory (CFT) results on the finitesize energy spectrum [11,12], they found that the BKT transition can be identified with a level crossing between a certain pair of the energy levels, canceling the loga-rithmic corrections. This allows a surprisingly accurate determination of the BKT transition point with exact numerical diagonalization of rather small systems.
However, the applications of Level Spectroscopy have been limited to 1D quantum BKT system such as quantum spin chains so far. While it should be applicable to the spectrum of the transfer matrix for the classical 2D XY model, the Level Spectroscopy has not been applied there, because of the difficulty in calculating the spectrum for the system with continuous variables θ j . In this paper, we demonstrate a successful implementation of Level Spectroscopy on the classical 2D XY model, based on the Tensor Network Renormalization (TNR) scheme [13][14][15][16][17][18]. The TNR enables to obtain a precise spectrum of the transfer matrix, to which the Level Spectroscopy can be applied. As in the case of the 1D quantum system, it allows a very accurate determination of the critical point by removing the logarithmic corrections from systems of moderate sizes which can be described by a tensor network with a finite bond dimension. On the other hand, our TNR study covers larger system sizes than those in the existing Level Spectroscopy studies on 1D quantum systems. We find a new feature of the finitesize scaling, that leads to a further improvement of Level Spectroscopy. Moreover, we can also visualize the celebrated Kosterlitz RG flow of the BKT transition from numerical data, for the first time to our knowledge [19] II. SU(2) SYMMETRY ON THE BKT LINE The low-temperature critical phase of the BKT transition is described by the free boson field theory in 1 + 1 dimensions, also known as Tomonaga-Luttinger Liquid: [20][21][22].
where K is called the Luttinger parameter and φ is a dual field of θ with compactifications φ ∼ φ + π and θ ∼ θ + 2π. Π = 1 πKφ is the canonical conjugate field of φ. Note that we deal with a continuous field θ(x) instead of the lattice variables θ i . (We follow the convention of Ref. [23]. See also the Appendix.) The compactification of the fields entails the existence of vertex operators V m,n = : e imθ+i2nφ :, which have the conformal weights It is helpful to first consider the BKT transition in the quantum spin-1 2 XXZ chain. There, the single vertex operators V 0,±1 are forbidden [24,25] and the BKT transition is driven by the double vertex operators V 0,2 +V 0,−2 . The Hamiltonian is then modified to the sine-Gordon model as, where K = −2∂ µ φ∂ µ φ shifts the Luttinger parameter from K = 1 2 and H represents the RG-irrelevant perturbations. The RG equation is again Eq. (2) and g = y K = y V is the BKT transition line. This corresponds to the isotropic XXX chain and the system has SU(2) symmetry. More precisely, H TLL K=1/2 turns into SU(2) 1 WZW model and Eq. (4) can be written with SU(2) currents as [26,27] Ignoring the RG-irrelevant perturbation H for the moment, exactly at the BKT transition, t = 0 and the effective Hamiltonian has a SU(2) symmetry. Consequently, all the finite-size energy levels can be classified in terms of the representation of the SU(2). Therefore, the BKT transition point can be identified by the degeneracy (level crossing) of the finite-size energy levels. This is nothing but the Level Spectroscopy method proposed and developed in Refs. [6][7][8][9][10][28][29][30]. Thanks to the emergent SU(2) symmetry, the transition point determined by the Level Spectroscopy is exact in all orders of the marginal coupling g, as long as the irrelevant perturbation H is negligible. In CFT, there is a correspondence between the finite-size energy levels and local fields. The lowest excited states of H WZW k=1 correspond to the "spin-wave" operators W ±1 = e ±iθ and the single vortex operators V ±1 = e ±2iφ , all with the scaling dimension 1/2. These 4 states are split into a singlet and a triplet by the SU(2) symmetric marginal perturbation g.
In the 2D classical XY model, the single vortex operator V ±1 = e ±2iφ is not forbidden in the Hamiltonian and indeed it is what drives the BKT transition, rather than the double vortex operator V ±2 . Nevertheless, the effective field theory in terms of boson field is equivalent to Eq. (4), with the replacement 2φ → φ and θ → 2θ. This implies that the fixed point Hamiltonian for the BKT transition point has the Luttinger parameter K = 2 instead of K = 1/2. It appears that the effective field theory in this case no longer has the SU(2) symmetry and the Level Spectroscopy may not apply. However, in Ref. [9], for Level Spectroscopy of the corresponding class of quantum spin chains, "half-vortex operators" V ±1/2 were introduced by twisting the boundary condition. At the BKT transition, one of the half-vortex states V s 1/2 = √ 2 sin φ becomes degenerate with the spin-wave states W ±2 = e ±2iθ to form an SU(2) triplet, thereby enabling the application of the Level Spectroscopy. Thus, in order to apply the Level Spectroscopy to the 2D classical XY model, we need to calculate the spectrum of the transfer matrix (which corresponds to the energy levels of 1D quantum Hamiltonian), under the periodic and twisted boundary conditions.

III. LEVEL SPECTROSCOPY WITH TENSOR NETWORK RENORMALIZATION
Now we demonstrate a successful implementation of Level Spectroscopy based on TNR for the XY model in Eq. (1). First, the partition function of the model is represented in terms of a tensor network, using a series expansion [31,32] where I n is the modified Bessel functions of the first kind.
In the practical calculations, we cut off the sum at |n| = 15 which is verified as sufficient.
In the tensor network representation of the XY model on the square lattice, the tensor on each site has four "legs" corresponding to the interactions with the four nearest neighbors. If we contract the legs of L horizontally aligned tensors with the periodic boundary condition, the remaining tensor is nothing but the transfer matrix in the vertical direction for the system of width L. As the transfer matrix corresponds to the Hamiltonian of the 1D quantum system, the Level Spectroscopy can be applied to its spectrum. In practice, the transfer matrix is obtained by contracting horizontal legs of a single tensor obtained after N steps of TNR. Since each step of TNR corresponds to rescaling of the lattice spacing by √ 2, this procedure gives the transfer matrix for the system of width L = √ 2 N . The eigenvalues λ n (L) thus obtained are related to the energy spectrum E n (L) of the corresponding 1D quantum Hamiltonian as where we define the rescaled energy levels x n by E n (L) − E 0 (L) = 2πx n (L)/L [33]. In this way, we can read off x W±2 . On the other hand, x V s 1/2 appears in the spectrum of the transfer matrix with the twisted boundary condition. By introducing a twist of angle 2π/L in the original model as we can also obtain the transfer matrix spectrum for width L with the twisted boundary condition. In particular, x W±2 and x V s 1/2 appear as the fourth/fifth(degenerate) and the leading eigenvalues in the transfer matrix in the periodic and twisted boundary conditions, respectively in the vicinity of T c . The scaling dimensions of the twisted boundary condition is measured against the leading eigenvalue of the periodic boundary condition.
Let x W2 and x V s 1/2 be the rescaled energy levels corresponding to the SU(2) triplet. While x O is given by the scaling dimension of the operator O in a pure CFT without any perturbation, it receives corrections from the irrelevant perturbations (y K , y V and H ) and depends on the system size L. Up to the first order in t and H (but to all orders in g), and using t ∝ T − T c , we find on the BKT transition line t = 0. δ W and δ V represent the first order corrections due to the irrelevant perturbation H . In the Level Spectroscopy, the transition temperature is first estimated by the crossing point T * between the two levels x W2 and x V s 1/2 . As far as the irrelevant perturbation H is ignored, it gives the exact transition temperature which corresponds to t = 0, removing the notorious logarithmic corrections to the all order. However, because of the irrelevant perturbation H , the crossing point α W (L) and α V (L) can be extracted as the slope of the corresponding levels in Fig. 2 for each size L. Since the leading irrelevant perturbations (other than the  Table I).
marginally irrelevant K and V ) to the fixed point Hamiltonian are the square of the stress-energy tensor with the scaling dimension 4, the standard CFT analysis [11] implies δ W ∝ δ V ∝ 1/L 2 . Thus, the level-crossing point T * (L) between x W2 and x s V 1/2 obtained for the system size L is related to the true transition temperature T c as where α(L) = α W (L) − α V (L). Therefore, the transition point can be extrapolated linearly by plotting T * against 1 α(L)L 2 . The original Level-spectroscopy implicitly assumed that α(L) is independent of L. This was a reasonable assumption because the system size employed there only ranged from L = 4 to 12. However, in our study, the system size is extended up to L = 32 and the size dependence of α(L) is not negligible. Indeed, in our numerical estimates of α(L), we found a weak dependence on L, which is consistent with the asymptotic behavior α(L) ∼ log L determined by the RG equation (2). The left panel of Fig. 3 exhibits the size dependence of the level-crossing temperature T * . The data show a good agreement with the theoretical prediction. The data for L = 4, which appears off from Eq. (10), are presumably affected by more irrelevant perturbations with the scaling dimension 6 and higher. Thus we carry out the linear extrapolation from the data for L = 8, 16, 32 where Eq. (10) holds almost perfectly to obtain T c . On the other hand, any practical calculation based on a tensor network is inherently limited by a finite bond dimension. The TNR is also accurate for the system sizes only up to the maximum correlation length imposed by the finite bond di-mension. In the present Level Spectroscopy approach, however, we can obtain very accurate results from only moderately large systems up to L = 32, where the TNR with the bond dimension D = 48 is sufficient [34]. This is manifest in the right panel of Fig. 3, which shows the extracted T c as a function of the bond dimension D: the dependence on D saturates at D ∼ 40, and the BKT transition temperature is identified as T c = 0.892943(2) at D = 48. As shown on Table I, our result has a higher precision than previous studies by an order of magnitude,

IV. RG FLOW OF THE XY MODEL
Finally, one can extract the coupling constants of the sine-Gordon model to visualize the RG flow. Ignoring H here, up to the second order in y K and y V , the lowest Monte Carlo(1979) [35] 0.89 Monte Carlo(2005) [36] 0.8929(1) Monte Carlo(2012) [37] 0.89289(6) Monte Carlo(2013) [38] 0.8935(1) Series expansion(2009) [39] 0.89286(8) HOTRG(2014) [40] 0.8921(19) VUMPS(2019) [32] 0.8930(1) HOTRG(2020) [41] 0.89290(5) present work 0.892943(2) rescaled energy levels are given as (13) where V c 1/2 = √ 2 cos φ is the remaining singlet. We read off x W±2 from Ref [42] and determine x V s 1/2 by imposing a restriction by the symmetry. (See the Appendix) Then, the running coupling constants are extracted from the numerically obtained energy levels as which are valid up to O(y 3 ). Fig. 4 shows the obtained RG flow of the XY model in the vicinity of the transition temperature. The theoretical trajectories of Eq. (2) are hyperbolas y 2 V − y 2 K =const. and our result matches perfectly with it: The low-temperature phase shown with blue colors flows to the y V = 0 critical line, whereas the high-temperature region does not. In the middle of two phases are the RG flow of the critical temperature T c = 0.893 marked with yellow dots, which is on top of the blue BKT line. This estimate of T c is consistent with our earlier estimate shown in Table I [43].

V. CONCLUSION
We analyzed the eigenspectrum of the renormalized tensors for the classical 2D XY model at finite renormalization steps, in terms of finite-size scaling of CFT. The BKT transition is described in terms of two marginal couplings y K (spin wave stiffness) and y V (vortex fugacity), and the transition point can be identified with y K = y V where a hidden SU(2) symmetry emerges. By exploiting the SU(2) symmetry, we determine the transition temperature with a record precision from the spectrum, improving the Level Spectroscopy method developed for 1D quantum systems. Furthermore, we note that SU(2) symmetry is no longer a requisite by regarding Level Spectroscopy as determination of transition temperatures based on the coupling constants of the underlying field theory extracted from the finite-size spectrum. One can extract the running coupling constants from a finite-size spectrum by following Eqs. (B12, B13, B18) in the Appendix, and tracking the RG flow exhibits the phase separation visually similar to Fig. 4. Our method is advantageous because finite-size effects are substantially reduced by the removal of the logarithmic corrections and the effects of the finite bond-dimension are suppressed at the moderate system sizes needed for our analysis. We also numerically tracked the evolution of the two coupling constants y K and y V as the system size is increased The original Level Spectroscopy was first developed to determine the transition temperature of the BKT transition for quantum spin chains. In particular, they investigated the effective Hamiltonian for the spin 1/2 XXZ spin chain. It has a phase transition (dimer to CDW) and the gapless phase is described by the Tomonaga-Luttinger liquid (TLL). The Lagrangian of the gapless phase is essentially where K is called Luttinger parameter. The bosonic field φ is compactified in the bulk as The canonical quantization of the field theory Eq. (A1) is done by requiring the canonical conjugate field The resulting Hamiltonian is Following the equation of motion of the Lagrangian Eq. (A1) we can decompose the field φ into the right-and leftmover as The dual field θ is then defined as θ is also compactified as θ ∼ θ + 2π.
In physical applications, often θ represents the angular variable corresponding to the microscopic U (1) symmetry, and it is indeed the continuous counterpart of the lattice variable θ i of the classical XY model in Eq. (1) of the main text. (In the next section, we apply Wick rotation as τ = it, z = x + iτ, andz = x − iτ so as to consider the CFT on the plane.) Since the filed theory Eq. (A1) is scale-invariant, it describes the critical and gapless phases of 2D classical systems and 1D quantum phases, respectively. However, the lattice models, generally speaking, contain perturbations at finite sizes and the XXZ spin chain is no exception. The field theoretical description of the system is the sine-Gordon model in the long wave-length limit near the critical parameters as, Table II. The first and second terms in the perturbations represent the renormalization of K and the creation of vortex-antivortex pairs, respectively. Kosterlitz derived the RG equations, and the transition occurs when g = y K = y V . On the BKT line y K = y V , the system restores SU(2) symmetry [44,45] and the SU(2) triplet V s 1 = sin(2φ) and W ±1 = e ±iθ becomes degenerate to the higher order loops of of g.
The TLL is an example of Conformal Field Theory (CFT) in 1 + 1 dimensions. Cardy [11] showed that, in an unperturbed CFT, the energy eigenvalues E n of the primary states in a quantum Hamiltonian of length L is given as where E 0 is the ground state energy and x n is the scaling dimension of the corresponding primary operator. Furthermore, perturbations to the CFT Hamiltonian give corrections to this relation [12]. Conversely, by analyzing the finite-size energy spectrum, one can estimate the perturbations to the CFT. This is the foundation of the Level Spectroscopy method.

Derivation of Level Spectroscopy
In the original papers [6][7][8][9][10] on Level Spectroscopy, the degeneracy of V s 1 and W ±1 was discussed in terms of the correspondence between the sine-Gordon model and the SU(2) Thirring model [44,45]. Here, we present the derivation of Level Spectroscopy in more intuitive manner.
The comformal weights of : e imθ+i2nφ : are h m,n = The scaling dimension is defined as a sum of conformal weights as x m,n = m 2 4K + n 2 K as in the main text. At K = 1/2, (m, n) = (±1, ±1) and (m, n) = (±1, ∓1) become currents respectively and the Hamiltonian gets equivalent to SU(2) 1 WZW model. The currents respect the SU(2) current algebra as where k = 1 and i = 0, 1, 2(or equivalently i = z, x, y) in the current case. The notations are listed in Table II and one can easily check that they follow Eq. (A5). Since the algebra is a representation of SU(2), it can be interpreted as a spin system. For instance, |↑ L is e i2ϕ |0 because its SU(2) charge is whereas J + |↑ = 0 because it has no pole as J + (z)e i2ϕ(w) ∼ (z − w)e i(4ϕ(z)+2ϕ(w)) . Using these notations, Eq. (A3) at g = y K = y V can be re-written as where H W ZW k=1 is a sum of the Sugawara energymomentum tensors as [46]. The lowest lying eigenstates of dx 2π (J L · J R ) are the spin singlet and triplets with their eigenvalues − 3 4 and 1 4 . Using Table II, the primary operators corresponding these states are As J + commutes with J 2 , V s 1 and W ±1 are degenerate. In particular, the energy levels of the singlet and triplets to the first order can be evaluated as Thus, the scaling dimensions of the triplets are 1 2 − 1 4 g + O(g 2 ).

Implementation of Level Spectroscopy
For small system sizes, there are also irrelevant operators that do not respect SU(2) symmetry due to the lattice anisotropy. The leading irrelevant ones are T 2 ,T 2 , TT , where T andT are holomorphic and antiholomorphic components of the energy-momentum tensor (not to be confused with the temperature T ). Since they have scaling dimension 4, it splits the energy levels of the triplets by ∼ 1/L 2 . While W ±1 remain degenerate, V s 1 is no longer at the same energy level. Nonetheless, the OPE coefficients with these irrelevant operators must be almost the same (conformal weights are almost the same). Thus, the level-crossing point is closed to the true transition point even at small L.
In the case of the 2D classical XY model, as discussed in the main text, the BKT transition is described by making the replacements on that for the S = 1/2 XXZ chain: 2φ → φ, θ → 2θ, and K = 1/2 → K = 2. In this context, W ±2 and V s 1/2 should be degenerate and form a SU(2) triplet, where V s 1/2 is one of the linear combinations of V ±1/2 corresponding to the insertion of the π-twist operator (antiperiodic boundary condition).
Appendix B: Calculation of yK and yV

Lukyanov's result
Combining the exact Bethe ansatz solution and the low-energy effective field theory, Lukyanov [42] studied the"vacuum" (ground-state) energy of the XXZ chain for a given total magnetization under the twisted boundary condition with the twist parameter 2πθ. (Here we introduce θ as the twist parametrizing the boundary condition as in Ref. [42]. This should not be confused with the use of θ as a field variable elsewhere in this paper.) His remarkable results have been verified by several numerical studies. While the results are limited to θ < 1 in Ref. [42], we may access the excited states corresponding to V ±1 under the untwisted (periodic) boundary condition. However, in this limit θ = 1, the mixing between V ±1 states by the allowed vortex perturbation V c 2 must be taken into account, as it was indeed suggested by Lukyanov [42]. Below we will demonstrate that the mixing is necessary to reproduce the split between triplet and singlet levels at the BKT transition, which is required by the emergent SU(2) symmetry as discussed in the main text.
Lukyanov's results can be also translated into our problem of the classical 2D XY model by the simple replacement 2φ → φ as discussed above. In the following, we present an analysis in the context of the classical 2D XY model. According to Lukyanov, the scaling dimension to the third order is where x m,n is the scaling dimension of : e imθ+i2nφ :. In particular, x 0,1/2 = 1 2 Boson field Definition Currents (k = 1 WZW) Corresponding operator The highest weight Corresponding state φ ∼ φ + π ϕ(z) +φ(z) In Eq. (B4), however, the energy repulsion between x s 0,1/2 and x c 0,1/2 , due to the cosine term, is not considered here. Hence, we shall determine it in order to calculate y K and y V to the second-order.

Calculation of
We deduce the energy level of V s 1/2 using two facts: · W ±2 and V s 1/2 are exactly degenerate on the BKT line. · x V s 1/2 + x V c 1/2 = 2x 0,1/2 , which is Eq. (B4). First, from Eq. (B3) we find the energy level of the triplets on the BKT line g = y K = y V is Then, using the second fact, we find that Eq. (B5) is larger than Eq. (B4) by b. Simple calculation leads to the power expansions of b with g as Given that b should be odd under y V → −y V the explicit form of b is deduced as The final forms of the singlet and triplets are 3. Perturbative calculation to the first order based on Conformal Field Theory Generally speaking, x n (L) depends on the system size due to the relevant/irrelevant perturbations to the fixedpoint Hamiltonian. Let the Hamiltonian with the system size L beĤ whereĤ * is the Hamiltonian of the fixed point,Φ j (x) is an operator with the scaling dimension x j and g j is a corresponding coupling constant. The scaling dimensions of a finite-size system is then described by c ijk is the operator product expansion(OPE) coefficient [11,12]. Following a standard calculation from Tomonaga-Luttinger liquid, we find that the OPE coefficients as c KW±mW±m = − m 2 4K , c KV±nV±n = n 2 K, c V c Comparing Eq. (B14) with Eq. (B12), we identify g K = y K 4π and g V = y V 2 √ 2π . Substituting them into Eq. (B13), we find This is in agreement with Eqs. (B8-B11) to the first order. The scaling dimension x n (L) = 1 2π ln(λ 0 /λ n ) in the main text changes as we renormalize the tensor(as we change the system sizes L), and so are y K (L) and y V (L). Therefore, we can compute y K at each scale by measuring the deviation of x W±2 from 1 2 for example. The identification of the relevant energy levels can be done by comparing the exact values of the scaling dimensions in the UV/IR fixed point. This approach is applicable as long as the Hamiltonian is in the vicinity of the fixed-point. Combined with TNR, our approach would potentially become a powerful method to quantitatively compute running coupling constants of the field theory from lattice models near the fixed points.