Approaching the Kosterlitz-Thouless transition for the classical XY model with tensor networks

TL;DR

This paper employs tensor-network methods, specifically uniform matrix product states with non-abelian symmetry, to accurately analyze the Kosterlitz-Thouless transition in the 2D XY model, capturing key universal and critical properties.

Contribution

It introduces a tensor-network approach with non-abelian symmetry to study the Kosterlitz-Thouless transition, providing high-precision estimates of critical behavior.

Findings

Universal drop in spin stiffness at critical point

Exponential divergence of correlation length in high-temperature phase

High-precision estimate of critical temperature

Abstract

We apply variational tensor-network methods for simulating the Kosterlitz-Thouless phase transition in the classical two-dimensional XY model. In particular, using uniform matrix product states (MPS) with non-abelian O(2) symmetry, we compute the universal drop in the spin stiffness at the critical point. In the critical low-temperature regime, we focus on the MPS entanglement spectrum to characterize the Luttinger-liquid phase. In the high-temperature phase, we confirm the exponential divergence of the correlation length and estimate the critical temperature with high precision. Our MPS approach can be used to study generic two-dimensional phase transitions with continuous symmetries.

Tables

Table 1. Table 1: The leading eigenvalue (per site) of the XY transfer matrix in the critical phase ( T = 0.8 𝑇 0.8 T=0.8 ) as obtained by imposing an MPS approximation for the fixed point, where we have imposed different symmetry properties on the virtual legs of the MPS.

	$D$	$\lambda$
no symmetries	19	2.5869206
integer charges	45	2.5869172
half-integer charges	46	2.5869184

Table 2. Table 2: Numerical estimates for the critical temperature

Monte Carlo (1979) [10]	0.89
Monte Carlo (2005) [11]	0.8929(1)
series expansion (2009) [13]	0.89286(8)
Monte Carlo (2012) [12]	0.89289(6)
Monte Carlo (2013) [14]	0.8935(1)
tensor-renormalization group (2014) [25]	0.8921(19)
uniform MPS (current work)	0.8930(1)

Equations

$$H=-\sum _{⟨ij⟩}\cos \left({\theta }_i-{\theta }_j\right)-h\sum _i\cos ({\theta }_i),$$

$$\underset{T\to T_c^{-}}{\lim }\rho (T)=\frac{2T_c}{\pi },$$

$$\xi (T)\propto \exp \left(\frac{b}{\sqrt{T-T_c}}\right),T\to T_c^{+},$$

$$Z=\prod _i\int \frac{\mathrm{d}{\theta }_i}{2\pi }\prod _{⟨ij⟩}{\mathrm{e}}^{\beta \cos ({\theta }_i-{\theta }_j)}\prod _i{\mathrm{e}}^{\beta h\cos ({\theta }_i)}.$$

$${\mathrm{e}}^{x\cos ({\theta }_i-{\theta }_j)}=\underset{N\to \infty }{\lim }\sum _{n=-N}^N{I}_n(x){\mathrm{e}}^{in({\theta }_i-{\theta }_j)},$$

$$Z=\underset{N\to \infty }{\lim }\prod _{l\in \mathcal{L}}\left(\sum _{n_l=-N}^N{I}_{n_l}(\beta )\right)\prod _s{F}_{n_{s,1},n_{s,2}}^{n_{s,3},n_{s,4}}$$

$$F_{n_1,n_2}^{n_3,n_4}=\int \frac{\mathrm{d}\theta }{2\pi }{\mathrm{e}}^{\beta h\cos \theta }{\mathrm{e}}^{i\theta (n_1+n_2-n_3-n_4)}.$$

$$Z=\;\vbox{\hbox{\includegraphics[scale={0.32},page=1]{./intro.pdf}}}\;,$$

$$\;\vbox{\hbox{\includegraphics[scale={0.32},page=2]{./intro.pdf}}}\;=\left(\prod_{i=1}^{4}I_{n_{i}}(\beta)\right)^{1/2}F_{n_{1},n_{2}}^{n_{3},n_{4}}$$

$$T(\beta,h)=\;\vbox{\hbox{\includegraphics[scale={0.32},page=3]{./intro.pdf}}}\;\;,$$

$$f(\beta ,h)$$

$$=\underset{N_x{N}_y\to \infty }{\lim }\left(-\frac{1}{\beta }\frac{1}{N_x{N}_y}\log \left(T(\beta ,h{)}^{N_y}\right)\right)$$

$$=-\frac{1}{\beta }\log \lambda (\beta ,h).$$

$$T(\beta ,h)|{\Psi }_{\beta ,h}⟩=\Lambda |{\Psi }_{\beta ,h}⟩,$$

$$\ket{\Psi(A)}=\;\vbox{\hbox{\includegraphics[scale={0.32},page=1]{./fixedpoint.pdf}}}\;.$$

$$\;\vbox{\hbox{\includegraphics[scale={0.32},page=2]{./fixedpoint.pdf}}}\;\propto\\ \;\vbox{\hbox{\includegraphics[scale={0.32},page=1]{./fixedpoint.pdf}}}\;.$$

$$\underset{A}{\max }\frac{⟨\Psi (\bar{A})|T(\beta ,h)|\Psi (A)⟩}{⟨\Psi (\bar{A})⟩\Psi (A)}.$$

$$\Lambda$$

$$\displaystyle=\;\vbox{\hbox{\includegraphics[scale={0.32},page=3]{./fixedpoint.pdf}}}\;,$$

$$\lambda=\rho\left(\;\vbox{\hbox{\includegraphics[scale={0.32},page=4]{./fixedpoint.pdf}}}\;\right).$$

$$\rho\left(\;\vbox{\hbox{\includegraphics[scale={0.32},page=5]{./fixedpoint.pdf}}}\;\right)=1.$$

$${F_{n_1,n_2}^{n_3,n_4}|}_{h=0}$$

$$={\delta }_{n_1+n_2}^{n_3+n_4}.$$

$$\;\vbox{\hbox{\includegraphics[scale={0.32},page=4]{./intro.pdf}}}\;=\delta_{n_{1}}^{n_{2}}n_{1},$$

$$U(\theta )={\mathrm{e}}^{i\theta {\sum }_j{Q}_j},\left[T(\beta ),U(\theta )\right]=0.$$

$$\;\vbox{\hbox{\includegraphics[scale={0.32},page=6]{./intro.pdf}}}\;=\delta_{n_{1}}^{n_{2}}\exp\left(i\theta n_{1}\right),$$

$$\;\vbox{\hbox{\includegraphics[scale={0.32},page=6]{./fixedpoint.pdf}}}\;\\ =\;\vbox{\hbox{\includegraphics[scale={0.32},page=1]{./fixedpoint.pdf}}}\;$$

$$\;\vbox{\hbox{\includegraphics[scale={0.32},page=7]{./fixedpoint.pdf}}}\;=\;\vbox{\hbox{\includegraphics[scale={0.32},page=8]{./fixedpoint.pdf}}}\;.$$

$$⟨h({\theta }_j)⟩=\frac{1}{Z}\prod _i\left(\int \frac{\mathrm{d}{\theta }_i}{2\pi }\right){\mathrm{e}}^{-\beta E(\{{\theta }_i\})}h({\theta }_j),$$

$$\braket{h(\theta_{i})}=\frac{1}{Z}\left(\;\vbox{\hbox{\includegraphics[scale={0.32},page=1]{./observ.pdf}}}\;\right)$$

$$\;\vbox{\hbox{\includegraphics[scale={0.32},page=2]{./observ.pdf}}}\;=\left(\prod_{i=1}^{4}I_{n_{i}}(\beta)\right)^{1/2}\\ \times\int\frac{\mathrm{d}\theta}{2\pi}\;\mathrm{e}^{\beta h\cos\theta}\mathrm{e}^{i\theta(n_{1}+n_{2}-n_{3}-n_{4})}h(\theta).$$

$$\braket{h(\theta)}=\frac{\left(\;\vbox{\hbox{\includegraphics[scale={0.32},page=3]{./observ.pdf}}}\;\right)}{\left(\;\vbox{\hbox{\includegraphics[scale={0.32},page=4]{./observ.pdf}}}\;\right)}\;.$$

Full text

I Introduction

In contemporary theoretical physics the interplay between symmetry and dimensionality is widely appreciated for giving rise to fascinating physical phenomena. Historically, one of the crucial results in establishing this viewpoint was the phase transition in the two-dimensional classical XY model. This model is introduced by placing continuous angles $\{\theta_{i}\}$ on a square lattice and characterizing their interactions by the hamiltonian

[TABLE]

where $\braket{ij}$ labels all nearest-neighbour pairs and we introduce an external magnetic field $h$ for further reference.

The first ingredient for understanding the phase diagram of the XY model (without magnetic field) is the Mermin-Wagner theorem [1, 2]: No conventional long-range order can exist at any finite temperature in this model, because of the proliferation of spin-wave excitations in two dimensions. Still, one expects a phase transition in this system. At low temperatures, the system can be described by a simple continuum field theory with algebraically decaying correlation functions. A high-temperature expansion, however, suggests an exponential decay of correlations in this system at sufficiently high temperature. As Berezinskii [3] and Kosterlitz and Thouless [4] (BKT) have showed, the phase transition between the low- and high-temperature regime is driven by the unbinding of vortices and is, therefore, of a topological nature. The transition is captured by a renormalization-group (RG) analysis [5, 6], where in the critical phase the long-wavelength properties are described by a free bosonic field theory with continuously varying critical exponents, until vortices become relevant and drive the system into a gapped phase at a critical temperature $T_{c}$. The quantity that characterizes the transition (in the absence of a local order parameter) is the spin stiffness $\rho(T)$, which discontinuously drops [7] from a finite value in the critical phase to zero in the gapped phase; from the renormalization analysis, the size of the drop is known to be given by

[TABLE]

but the value of $T_{c}$ is not known exactly. When approaching the critical point from the gapped side, the BKT transition is characterized by an exponential divergence of the correlation length [5, 8],

[TABLE]

with $b$ a non-universal parameter.

The BKT phase transition in the XY model has been a notoriously hard case for numerical simulations, because of the exponentially diverging correlation length and the ensuing logarithmic finite-size corrections around the phase transition. Nonetheless, shortly after the theoretical work, Monte-Carlo simulations [9, 10] provided considerable evidence for the correctness of the theory and rough estimates for the transition temperature. In recent years, these estimates were significantly improved [11, 12], obtaining a value around $T_{c}\approx 0.8929$, in agreement with high-temperature expansions [13]. Still, these results depend heavily on assumptions about the logarithmic finite-size corrections and an improved extrapolation [14] places the transition temperature at a higher value of $T_{c}\approx 0.8935$. Crucially, the best estimates use the universal value for the spin stiffness at the phase transition for pinpointing the critical point.

Tensor networks [15, 16] provide an original framework for capturing the symmetry-dimensionality interplay, both from the theoretical and the numerical side. Although orginally devised for capturing the entanglement in strongly-correlated quantum lattice systems, tensor networks are increasingly being applied to problems in statistical mechanics. Since the framework is entirely different from traditional approaches such as Monte-Carlo sampling, it can shed a new, entanglement-based, light on statistical-mechanics problems. Indeed, tensor networks encode all physical properties of a given system into local tensors, and they allow to understand the relation between physical symmetries of the local degrees of freedom and the global properties of the system in a transparent way. One particular example is the classification of symmetry-protected topological phases in one-dimensional quantum systems, which is brought back to the symmetry properties of the local tensors that make up a matrix product state (MPS) [17, 18]. This was extended to two dimensions, where the topological order of a wavefunction can be related to the symmetry properties of the local tensor in a projected entangled-pair state (PEPS) [19]. On the numerical side, physical symmetries can be explicitly incorporated in tensor-network algorithms and lead to an improved efficiency and performance [20, 21, 22, 23].

The encoding of symmetries in tensor networks is performed most elegantly when working directly in the thermodynamic limit, because one can purely focus on the symmetry properties of the bulk tensors without bothering about what happens at the boundaries of the system. Yet, applying uniform tensor networks with an explicit encoding of the physical symmetries seem to break down when considering critical phases. For example, in a critical phase a uniform MPS [24] typically favours an artificial breaking of a continuous symmetry, where the associated order parameter decreases very slowly as the bond dimension is increased. The apparent reason for this artificial symmetry breaking is that MPS have a built-in limitation for the amount of entanglement in the state, which makes it energetically favourable to break a continuous symmetry. This seems to imply that uniform tensor networks fail to capture the essential properties of critical phases with a continuous symmetry.

In this paper, we explore this question in more detail by investigating the precise sense in which uniform MPS capture critical phases with a continuous symmetry. As explained, the XY model serves as the paradigmatic example of a system where the absence of symmetry breaking leads to a critical phase, and, therefore, we take the XY model as our test case. In contrast to earlier tensor-network approaches for the XY model [25], we use uniform MPS methods for transfer matrices [26, 27, 28] as a means for characterizing the BKT phase transition. In the first two sections we explain the duality transformation [6] that allows us to define a row-to-row transfer matrix, approximate its fixed point as a uniform MPS and to compute local observables. In the next section, we focus on the spin stiffness as the characteristic quantity in the BKT phase transition. Afterwards, we use the Luttinger-liquid formalism to characterize the critical phase. Finally, in the last section, we focus on the gapped phase and locate the critical temperature with high precision.

II MPS for the XY transfer matrix

II.1 Partition function

We start by writing down the partition function for the XY model as a tensor network. The partition function at a given inverse temperature $\beta=1/T$ is given by

[TABLE]

In order to arrive at a tensor-network representation, we introduce a duality transformation [6] that maps the above partition function to a representation in terms of bosonic degrees of freedom on the links. Such a map is obtained by introducing the following decomposition on every link in the lattice

[TABLE]

where $I_{n}(x)$ are the modified Bessel functions of the first kind. Then, by integrating over all the $\theta$’s, the partition function is transformed into

[TABLE]

where $F$ is a four-index tensor

[TABLE]

The first product runs over all the links in the lattice, and $s$ labels all the sites in the lattice. We can now represent this partition function as a network of tensors,

[TABLE]

where every tensor $O$ is given by

[TABLE]

and the virtual legs $n_{i}$ have infinite dimension. In practice, however, it will be possible to truncate these indices without loss of accuracy. We have introduced arrows on the legs to indicate the signs in which the $n_{i}$’s appear in the $F$ tensor above.

The fundamental object in this representation of the partition function is the row-to-row transfer matrix $T(\beta,h)$

[TABLE]

which is an operator acting on an infinite chain of bosonic degrees of freedom. The value for the partition function and, therefore, the free energy of the model, is determined by the leading eigenvalue $\Lambda$ of the transfer matrix. Indeed, this leading eigenvalue is expected to scale as the number of sites per row, i.e. $\Lambda\sim\lambda^{N_{x}}$, such that the free energy per site is

[TABLE]

The eigenvector corresponding to the leading eigenvalue, refered to as the fixed point of the transfer matrix,

[TABLE]

will be of crucial importance in all computations. In a number of applications, it has been shown that fixed points of transfer matrices can be approximated accurately using the variational class of matrix product states (MPS) [26]. For translation-invariant transfer matrices, we can describe the fixed point as a uniform MPS described by a single tensor $A$,

[TABLE]

The tensor $A$ has virtual legs of dimension $D$, which we call the MPS bond dimension, and by repeating the tensor on every site and contracting over the virtual legs, we obtain a translation-invariant state. For simplicity, we take the virtual legs of the MPS to have no arrows. In terms of the fixed-point MPS, the fixed-point eigenvalue equation is rephrased as

[TABLE]

We aim at finding a tensor for which this eigenvalue problem is obeyed in an optimal way. Since the transfer matrix is hermitian, we can use the variational principle for establishing an optimization problem for the tensor $A$:

[TABLE]

This optimization problem can be efficiently solved using tangent-space methods for uniform MPS [26, 28]; in particular, we use the vumps algorithm [27] for finding the optimal MPS tensor. The eigenvalue, and therefore the free energy, is then obtained as the contraction of an infinite channel of $O$ tensors sandwiched between the fixed-point MPS and its conjugate,

[TABLE]

or we find $\lambda$ as the leading eigenvalue (spectral radius $\rho$) of the channel operator

[TABLE]

Here, we have assumed that the MPS itself is normalized as

[TABLE]

In Fig. 1 we have plotted the result for the free energy per site as obtained by a variational MPS simulation of the XY model in a certain temperature range and without magnetic field, where we use MPS with a bond dimension of $D=30$. It is clear that the free energy shows no signs of a phase transition as it is perfectly smooth everywhere. Note that, whereas the free energy cannot be directly computed using Monte-Carlo sampling, it appears as the fundamental quantity in the variational MPS setup.

II.2 Symmetries in the MPS representation

The case without magnetic field ($h=0$) is of particular importance, because the model is invariant under a global transformation $\theta_{i}\to\theta_{i}+\alpha$. This $\mathrm{U}(1)$ invariance is reflected in the tensor-network representation as a symmetry of the transfer matrix. Indeed, if $h=0$ we have

[TABLE]

such that the tensor $O$ has conservation of $\mathrm{U}(1)$ charges. We can introduce the operator $Q$

[TABLE]

which counts the charge on a given leg in the tensor network; this operator is the generator of the $\mathrm{U}(1)$ symmetry of the transfer matrix. Indeed, from the conservation property for the tensor $O$ the transfer matrix clearly commutes with the symmetry operation

[TABLE]

For future reference, we introduce the tensor $u(\theta)$

[TABLE]

such that the symmetry operation is $U(\theta)=\bigotimes_{i}u_{i}(\theta)$.

The Mermin-Wagner theorem now dictates that this $\mathrm{U}(1)$ symmetry cannot be broken at any finite temperature. On the level of the transfer matrix this implies that the leading eigenvector (fixed point) necessarily is invariant under the $\mathrm{U}(1)$ transformation. For the MPS approximation of the fixed point, this implies that

[TABLE]

for any $\theta$. The fundamental theorem of MPS now dicates that we can associate this symmetry of the MPS to a symmetry property of the local tensor $A$ [29]. Specifically, we have that the virtual legs of the MPS tensor transform themselves according to representations of the $\mathrm{U}(1)$ symmetry,

[TABLE]

In general, these representations on the virtual level can be projective. Imposing that the MPS is invariant under group transformations implies that the tensor $A$ has a certain block structure, where each block can be labeled by a choice of representations on each leg of the tensor.

Although the characteristic physics of the XY model is determined by the $\mathrm{U}(1)$ symmetry, we can also take into account the charge-conjugation symmetry $C$ that flips the sign of the charges in the transfer matrix. The total symmetry group that we consider is the (non-abelian) group $\mathrm{U}(1)\rtimes C$, which corresponds to the full $\mathrm{O}(2)$ symmetry of the XY model. There are two types of irreducible representations (irreps) of this $\mathrm{O}(2)$ symmetry. First, there are the irreps corresponding to integer charges $n=0,1,2,\dots$; for $n=0$ there are two one-dimensional irreps, whereas for $n>0$ they are all two-dimensional. Second, there are half-integer charges $n=\frac{1}{2},\frac{3}{2},\dots$, corresponding to projective representations; these irreps are all two-dimensional.

Since the transfer matrix only contains legs with integer representations, the ‘physical legs’ of the MPS all transform according to integer representations as well. This, however, still allows for the freedom on the virtual level of the MPS to impose either integer representations or half-integer representations, giving rise to two different classes of $\mathrm{O}(2)$ invariant MPS. In gapped phases, the occurence of either set of representations is known to characterize a trivial (integer) or symmetry-protected topological (half-integer) phase [17, 18]. In order to illustrate these different classes of MPS and to see which one is realized for the fixed point of the XY transfer matrix, we have performed simulations (i) without imposing symmetries on the legs, (ii) with integer irreps on the legs, and (iii) with half-integer irreps on the legs. In Table 1 we list the corresponding values for the transfer-matrix leading eigenvalue $\lambda$ at a point in the critical phase.

A first observation is that we find a higher value for $\lambda$ when no symmetries are imposed at significantly lower bond dimension. This shows that a variational MPS favours the breaking of the $\mathrm{U}(1)$ symmetry. As we have discussed in the introduction, this artificial symmetry breaking is also observed in other critical systems with a $\mathrm{U}(1)$ symmetry [24]. The explanation for this effect can be sought in the fact that MPS necessarily induces a finite correlation length. Indeed, simulating a critical system with MPS with a finite bond dimension can be thought of as slightly perturbing the system such that a gap is opened and a finite correlation length is introduced. From the perspective of an effective field theory describing the critical system, the MPS adds a relevant perturbation that opens up a gap in this field theory. In the $\mathrm{U}(1)$ phase, however, only symmetry-breaking terms are relevant, such that we expect that the MPS approximation induces an artificial symmetry breaking. This observation is confirmed when we compute the $\mathrm{U}(1)$ order parameter in the next section.

Secondly, upon imposing the $\mathrm{U}(1)$ or $\mathrm{O}(2)$ symmetry explicitly, we observe that the eigenvalue reaches a similar value at a comparable bond dimension in the two (normal and projective, resp.) sectors. If we again interpret the MPS approximation as introducing a gap, this result suggests that the critical $\mathrm{U}(1)$ phase can be perturbed into both an SPT phase and a regular phase. As we will see in Sec. V, in the gapped phase it can be determined without ambiguity what irreps should be chosen on the MPS virtual legs.

II.3 Computing local observables

Local observables such as the internal energy and the magnetization can be represented in the tensor-network language as follows. A generic one-angle observable at site $j$

[TABLE]

can, under the duality transform that we introduced above, be represented diagrammatically as

[TABLE]

with

[TABLE]

Using the MPS representation of the transfer-matrix fixed point, we can simplify this to

[TABLE]

The contractions of these infinite channels are evaluated by finding the leading eigenvectors of the channel operators.

Similarly, a generic nearest-neighbour two-angle observable at sites $j$ and $k$ is given by

[TABLE]

with

[TABLE]

Using the MPS fixed points that we have optimized earlier for computing the free energy – i.e. without using symmetries on the virtual level – we now compute the internal energy $e$, the entropy $s$ and the order parameter $o$

[TABLE]

as a function of temperature, and plot the results in Fig. 2. Again, the energy and entropy show no sign of the phase transition. Also, we observe that the value for the energy has already converged up to an error of $\epsilon=10^{-6}$ at a bond dimension of $D=30$; since the entropy is computed from the free energy and the internal energy, it has the same accuracy. As anticipated in the previous section, the order parameter shows a very large value in the critical region, which decreases very slowly as the bond dimension increases. This shows that the MPS breaks the continuous $\mathrm{U}(1)$ symmetry significantly in the critical phase, whereas in the gapped phase the symmetry is restored. The fact that the order parameter decays very slowly with increasing bond dimension shows that this is an essential property of MPS approximations for $\mathrm{U}(1)$ phases. On the other hand, from the convergence of the free and internal energy, we see that this does not prohibit an accurate evaluation of the system’s physical properties. Note that the order parameter is identically zero if we would impose the $\mathrm{U}(1)$ or $\mathrm{O}(2)$ symmetry on the virtual level of the MPS, but, as we have seen in Tab. 1, at a large variational cost in the free energy.

III The spin stiffness

The phase transition in the XY model can, in the absence of a local order parameter, be characterized by the so-called spin stiffness. This quantity is defined as the response to a twist field $\vec{v}$, which rotates the angles as

[TABLE]

If we take the twist field along the $y$ axis, this modifies the classical hamiltonian (without magnetic field) to

[TABLE]

On the level of the partition function, this introduces an extra phase factor on the vertical links,

[TABLE]

Incorporating these extra phase factors, we can represent the partition function as

[TABLE]

where the tensor $u_{v}$ was defined in Eq. (4). We find for the transfer matrix after the twist

[TABLE]

because $U(v)$ commutes with the transfer matrix. This implies that all eigenvectors of the transfer matrix remain unchanged after a twist and the eigenvalues are multiplied by a phase. In particular, the fixed point remains unchanged and, as the fixed point lives in the charge-zero sector, so does the leading eigenvalue. This, in turn, implies that the free energy density is a constant when applying a twist

[TABLE]

and that the spin stiffness, defined as

[TABLE]

is identically zero for all temperatures.

This is a surprising result, because $\rho$ is supposed to jump discontinuously at the transition. The reason for this discrepancy lies in the fact that we work directly in the thermodynamic limit. Indeed, the spin stiffness is typically defined in a system that is finite and has periodic boundary conditions in the direction of the twist. In that setting, applying the twist field is equivalent to imposing twisted boundary conditions – this can be observed from the above representation of the partition function, where the operator $U(v)$ can be pulled through the lattice – and is, therefore, a finite-size property. It is only after it has been defined on this periodic system, that the infinite-size limit can be taken. The above definition of the spin stiffness in terms of the infinite-size transfer matrix, on the other hand, assumes open boundary conditions. If the transfer matrix has no gap in the thermodynamic limit, both definitions are not equivalent.

We can work around this by computing the spin stiffness as a two-point function. First we note that we can rewrite the expression for the spin stiffness as

[TABLE]

where the twisted partition function can be written as

[TABLE]

with the twisted tensor

[TABLE]

In the tensor-network representation of $Z_{v}$ we can easily differentiate with respect to $v$. Indeed, the first derivative is given by

[TABLE]

with the tensor $R=\left.\frac{\partial O_{v}}{\partial v}\right|_{v=0}$, or

[TABLE]

For the second derivative we have to differentiate two different tensors, and twice the same tensor. The result is given by (the site $i$ is arbitrary)

[TABLE]

where we have introduced the tensor $S=\left.\frac{\partial^{2}O_{v}}{\partial v^{2}}\right|_{v=0}$, or

[TABLE]

This reduces the spin stiffness to a summation of two-point functions, and therefore has the form of a structure factor. We should note that bringing the factor $in_{3}$ down in the tensor is equivalent to introducing a factor $\sin(\theta_{i}-\theta_{j})$ in the partition function because of the identities

[TABLE]

The spin stiffness can therefore be brought into the form

[TABLE]

with

[TABLE]

In this form, it can be evaluated in Monte-Carlo simulations on a system with periodic boundary conditions without explicitly applying a twist [30, 14].

Again, the evaluation of this two-point function is identically zero on an infinite system with open boundary conditions in the presence of an unbroken $\mathrm{U}(1)$ symmetry. Indeed, if we represent the infinite upper plane of the above expression by the fixed point of the transfer matrix, we observe that

[TABLE]

which vanishes because the fixed point is $\mathrm{U}(1)$ symmetric. This is, of course, consistent with the vanishing of the spin stiffness as defined above.

For evaluating the spin stiffness, we can, however, introduce a magnetic field in the hamiltonian, which breaks the $\mathrm{U}(1)$ symmetry of the model and induces a gap. In the presence of this extra field, the spin stiffness as defined by the above structure factor [Eq. (6)] can be evaluated in an infinite system with open boundary conditions. Similarly to taking the infinite-size limit for a periodic system, we can then take the limit $h\rightarrow 0$ to obtain the result in the zero-field case. In Fig. 3 the results for the spin stiffness (as defined from the structure factor in Eq. (6)) as a function of temperature are given for a few values of the magnetic field. These results were obtained by using the channel-environment construction of Ref. 31 for evaluating structure factors in two-dimensional tensor networks. The figure shows that the drop in the spin stiffness becomes sharper as the magnetic field is decreased, but for very small values of $h$ the effects of a finite bond dimension become more pronounced.

IV Luttinger liquid mapping

In the previous section we have computed the spin stiffness by explicitly breaking the $\mathrm{U}(1)$ symmetry via an external magnetic field. In this section, we show that we can use the Luttinger-liquid formalism [32] to compute the same quantity without an explicit breaking of the symmetry, directly in the thermodynamic limit.

Early on, it was realized [6] that the proper effective field theory for the partition function is given by the sine-Gordon model [33], which is described by the one-dimensional quantum hamiltonian

[TABLE]

The first line is the hamiltonian of the Luttinger-liquid field theory [32] describing the critical spin-wave excitations, and the second line adds vortices to the picture. The microscopic $\mathrm{U}(1)$ symmetry of the model is reflected in the sine-Gordon field theory by the generator

[TABLE]

Since vortices are irrelevant in the critical phase, we expect that under an RG transformation the vortices will drop from the sine-Gordon hamiltonian and the parameters $u$ and $K$ will be renormalized. In other words, we expect that the low-energy properties of the XY model will be described by an effective Luttinger-liquid field theory

[TABLE]

with the effective parameters $\tilde{u}$ and $\tilde{K}$ determined by the inverse temperature $\beta$. As soon as the effective Luttinger parameter reaches $\tilde{K}=2$, however, vortices become relevant, leading to a gapped phase where the $\phi$ field is locked in the minima of the cosine term; the elementary excitations are kinks and anti-kinks between the different ground states.

We can lift this notion of an effective field theory to the level of a transfer matrix, where we formally introduce a Luttinger-liquid transfer matrix

[TABLE]

The idea is that this captures the low-energy behavior of the XY transfer matrix, where the parameters $\tilde{u}$ and $\tilde{K}$ depend on the inverse temperature $\beta$. We confirm the low-energy correspondence of the XY transfer matrix $T(\beta)$ with this effective Luttinger-liquid form in three different ways.

The first correspondence can be found by computing the low-lying spectrum of the XY transfer matrix, using the MPS quasiparticle ansatz for the low-lying excited states [26, 28]. We can label the excitations with a momentum $p$, yielding a dispersion relation $\omega(p)$ as the logarithm of its eigenvalues,

[TABLE]

and observe that we find a gapless spectrum with a linear dispersion relation $\omega(p)=\tilde{u}|p|$. Even stronger, we find the effective velocity $\tilde{u}=1$ for all values of the temperature within the critical phase (see Fig. 4(a)). It is, of course, expected that the effective field theory is isotropic, but this serves as an excellent test of our transfer-matrix approach.

A second test consists of measuring the central charge, which can be determined in an MPS simulation by comparing the scaling of the entanglement entropy and the correlation length as a function of the bond dimension. It is known that, for a system that is described by a conformal field theory the scaling is determined by the central charge $c$ as [34, 35]

[TABLE]

In Fig. 4(b) we show clear evidence for a central charge $c=1$, which is precisely the value for a Luttinger liquid.

An estimate for the effective Luttinger parameter $\tilde{K}$ is obtained by computing the response to a chemical potential. It is easily seen [32] that adding the generator of the $\mathrm{U}(1)$ symmetry,

[TABLE]

yields a compressilbity that is equal to

[TABLE]

We have an explicit microscopic form for the generator of the $\mathrm{U}(1)$, so we can explicity implement a chemical potential on the level of the microscopic transfer matrix as

[TABLE]

This shift in the transfer matrix does not affect its eigenvectors, but reshuffles the eigenvalues. Therefore the fixed point $\ket{\psi_{\mu}}$ will change as a function of $\mu$, and will give rise to a finite expectation value of the generator. The associated compressbility

[TABLE]

will therefore give a finite value. This, in turn, yields a direct estimate of the effective Luttinger parameter $\tilde{K}$. Since $Q$ commutes with the transfer matrix, the partition function will only be affected in second order, and the second derivative yields the same value for the compressibility,

[TABLE]

where $\lambda_{\mu}(\beta)$ is the usual eigenvalue per site of the transfer matrix

[TABLE]

We can compute this compressibility straightforwardly using the uniform MPS framework. Indeed, as explained above, the eigenvalue $\lambda_{\mu}(\beta)$ is the quantity that is variationally optimized in an MPS fixed-point simulation, whereas the charge density $\braket{Q_{i}}$ can be easily computed as an expectation value. In Fig. 4 it is shown that both quantities yield a consistent numerical value for the compressibility.

We should note, however, that the transfer matrix $T_{\mu}(\beta)$ is equivalent to the one of the twisted XY model but with an imaginary value for the twist field

[TABLE]

Therefore, the effective Luttinger parameter that we have defined here is related to the spin stiffness from before,

[TABLE]

This correspondence, which can be readily seen from the mapping of the XY model to the sine-Gordon field theory [6, 33], yields the famous [7] value for the spin stiffness at the critical temperature $\rho=2T_{c}/\pi$.

In the previous section, we had anticipated that computing the spin stiffness without introducing a symmetry-breaking term would not be possible in the thermodynamic limit directly. The reason that we are here able to compute the spin stiffness without breaking the $\mathrm{U}(1)$ symmetry consists in the fact that we have expressed it as a thermodynamic quantity (the compressibility) for which the extensivity properties of the uniform MPS simulations are ideally suited. This thermodyamic quantity, however, is necessarily formulated on the level of the transfer matrix after the duality transformation, since an imaginary twist does not translate to a realistic modification of the classical XY hamiltonian.

As a final signature of the Luttinger-liquid phase, we investigate the entanglement spectrum of the MPS fixed point. As was observed in Ref. 36, the low-lying part of the entanglement spectrum for a bipartition of the MPS, should resemble the energy spectrum of a boundary conformal field theory (CFT). In Fig. 5 we plot the entanglement spectrum of the MPS fixed point, where we have imposed $\mathrm{O}(2)$ symmetry on the MPS tensor such that we can label the spectrum with the appropriate quantum numbers. We observe that the spectrum has a quadratic envelope, and that we obtain an equidistant spectrum after rescaling the different sectors, in perfect correspondence with the spectrum of a free-boson boundary CFT. Moreover, from the rescaling parameter we can deduce an estimate of the Luttinger parameter $\tilde{K}$, which nicely converges to the same value as the one we find using the compressibility.

V The gapped phase

The characterization of the gapped phase of the XY transfer matrix using uniform MPS is a lot more straightforward. Indeed, we expect that the fixed point of a gapped transfer matrix can be approximated by an MPS with arbitrary precision. Therefore, we no longer expect that the MPS will spontaneously break the $\mathrm{U}(1)$ or $\mathrm{O}(2)$ symmetry, and we can safely use the fundamental theorem to realize that the virtual legs of the MPS should transform under (projective) representations as well (according to Eq. (5)).

As explained in Sec. II.2, we are ignorant on which representations should be chosen on the virtual legs of the MPS. For that reason, we first plot the entanglement spectrum of a fixed-point MPS around temperature $T=1.2$ without explicit symmetries on the MPS representation (see Fig. 6). We find both isolated and twofold degenerate Schmidt values, which point to integer representations on the virtual legs. Indeed, in Fig. 6 we plot the entanglement spectrum with integer representations on the virtual legs imposed, showing that the isolated values correspond to either of the two one-dimensional irreps with $n=0$ charge sector; the twofold degenerate ones correspond to $n=1,2,\dots$.

One of the hallmarks of the BKT transition is the exponential divergence of the correlation length when approaching the critical point from the gapped side [Eq. (2)]. Using MPS, we can confirm this behaviour and use this form to obtain an estimate for the critical point. The correlation length is a notoriously hard quantity to converge in MPS simulations, but using the second gap in the transfer matrix (typically denoted as $\delta$) it is possible to extrapolate its value in a reliable way [37]. In Fig. 7 we plot this extrapolation procedure for $T=0.93$, yielding an accurate value for a correlation length of more than a thousand sites. In Fig. 7 we then plot the extrapolated correlation lengths as a function of temperature, and fit this to

[TABLE]

where the extra terms are added to account for deviations away from the critical point. This fit yields an estimate for the correlation length of $T_{c}=0.8930(1)$, which agrees well with other numerical results (see Tab. 2). We should note, however, that estimates for the critical point can depend strongly on higher-order corrections to the scaling behavior – for the finite-size extrapolation of the spin stiffness this is clearly the case [14]. We leave the careful incorporation of higher-order contributions to the scaling of the correlation length, and a more accurate estimation of the critical point, for further study.

The gapped phase is further characterized by the low-lying spectrum of the transfer matrix, which we define as before in Eq. 7. In Fig. 8 we have plotted the spectrum at a temperature $T=1.3$, showing an isolated two-fold degenerate quasiparticle line; the excited states on this line carry $\mathrm{U}(1)$ quantum numbers $q=\pm 1$. Above this elementary one-particle excitation, we find the two-particle continuum and, interestingly, around momentum $p=\pi$ we find a slightly bound state below the continuum with charge $q=0$.

VI Outlook

In this paper we have investigated the classical two-dimensional XY model using uniform MPS methods. We have shown that an MPS approximation for the fixed point of the XY transfer matrix breaks the $\mathrm{U}(1)$ heavily in the critical phase, which is expected because the MPS always induces a finite correlation length in the system. In a similar vein, a uniform MPS calculation of the spin stiffness in the critical $\mathrm{U}(1)$ phase is always zero. The reason for the latter was sought in the fact that uniform MPS work in the thermodynamic limit directly, whereas the spin stiffness is a quantity that is necessarily defined in a finite periodic system; only for systems with a gap the two definitions intersect.

Nonetheless, we showed that uniform MPS are an ideal framework for characterizing the XY model and its phase transition. First of all, despite the fact that an MPS breaks $\mathrm{U}(1)$ symmetry, we can evaluate local observables with very high precision in the critical phase. Secondly, we have shown that the spin stiffness can be evaluated by introducing a small magnetic field $h$, and taking the limit $h\to 0$. The mapping to an effective Luttinger-liquid field theory can be made explicit using the MPS framework by computing the central charge, the dispersion relation and the compressibility; the latter, which is the response to a twist with imaginary magnitude, is used to find very accurate values for the effective Luttinger parameter $\tilde{K}$. In addition, we find that the entantglement spectrum of the MPS is in agreement with a boundary CFT spectrum. In the gapped phase, the MPS leaves the (non-abelian) $\mathrm{O}(2)$ symmetry unbroken, which we can use to find accurate values for the correlation length upon approaching the critical point; from fitting the exponential divergence of the correlation length we find $T_{c}\approx 0.8930$, in agreement with other numerical studies.

We expect that the current setup can be applied to other two-dimensional classical systems with continuous symmetries. Whereas the standard ferromagnetic Heisenberg model has no phase transition [38], the so-called $\mathrm{RP}^{2}$ models with a classical hamiltonian $H=-\sum_{\braket{ij}}\left(\vec{s}_{i}\cdot\vec{s}_{j}\right)^{2}$ (with $\vec{s}_{i}$ a three-dimensional unit vector) potentially hosts $\mathbb{Z}_{2}$ vortices that drive a phase transition [39]. A similar phase transition might be present in the frustrated antiferromagnetic Heisenberg model on the triangular lattice [40, 41]. Also, our methods can be readily applied for simulating KT transitions in one-dimensional quantum systems. The relation between these classical topological transitions and one-dimensional SPT phases [42] that we have investigated in this paper might prove very interesting in this context.

The results in this paper will prove instrumental in the program of simulating systems with unbroken continuous symmetries with uniform tensor networks and should, in particular, be useful in the study of two-dimensional quantum spin liquids with PEPS. Indeed, the norm of a PEPS can be naturally interpreted as a two-dimensional partition function and the question often poses itself in what phase the corresponding PEPS transfer matrix is. The paradigmatic example here is the resonating valence-bond (RVB) wavefunction on the square lattice, for which the transfer matrix is known to be in a $\mathrm{U}(1)$ phase [43, 44], and also other symmetric PEPS parametrization for (chiral) spin liquids seem to give rise to critical transfer matrices [45]. The relation between the critical properties of the transfer matrix, which are esssentially the properties of a two-dimensional classical system, the symmetries of the PEPS tensors, and the quantum properties of the PEPS wavefunction (physical correlation functions, entanglement spectra, etc.) remains, however, largely unexplored.

In a different direction, the two-dimensional partition functions that we have considered here, can be naturally lifted to the quantum level by promoting the charges on the bonds to quantum-mechanical degrees of freedom. Upon doing that, we find quantum-mechanical wavefunctions for $\mathrm{U}(1)$ gauge theories on the lattice. This PEPS construction can be generalized to a whole variational class of states that are ideally suited to study the phase diagram of two-dimensional lattice gauge theories, and, in order to understand the phase transitions, we will need the tools that were explored in this paper.

We acknowledge inspiring discussions with Nick Bultinck, Thierry Giamarchi, Jutho Haegeman, and Masaki Oshikawa. This work was supported by the Flemish Research Foundation, the Austrian Science Fund (ViCoM, FoQuS), and the European Commission (QUTE 647905).