Berezinskii-Kosterlitz-Thouless phase transition from lattice sine-Gordon model

We obtain nonperturbative results on the sine-Gordon model using the lattice field technique. In particular, we employ the Fourier accelerated hybrid Monte Carlo algorithm for our studies. We find the critical temperature of the theory based autocorrelation time, as well as the finite size scaling of the"thickness"observable used in an earlier lattice study by Hasenbusch et al.


Introduction
Many revolutions in theoretical physics occured in the early 1970s; one of these was the discovery that phase transitions did not always associate themselves with spontaneous symmetry breaking and long-range order. The physics of vortices, and the corresponding topological phase transitions, allowed one to avoid the theorems that opposed spontaneous symmetry breaking in two dimensions. Besides illustrating a defect driven topological phase transition, the related systems, the XY model and the two-dimensional (2d) Coulomb gas, became proving grounds for ideas related to the Wilsonian renormalization group. We illustrate the binding and unbinding of vortex -anti-vortex pairs on either side of the transition in Fig. 1.
The Berezinskii-Kosterlitz-Thouless (BKT) transition [4,5] was originally connected to the superfluid transition in two dimensions, such as the 4 He thin film. Later it was realized that it could also occur in superconducting thin films, in spite of being a charged superfluid. In this case the supercurrents screen the fluctuations. However, the screening length is given by Λ = λ 2 /d, where λ is the magnetic penetration depth and d is the film thickness. Thus for large disorder (hence large λ) and for sufficiently thin films, Λ can still be large enough for the algebraic order of the XY universality class to emerge. 1 The BKT transition is remarkably different from the usual Ginzburg-Landau transition. It is a transition without an order parameter. There is a line of conformal fixed points, rather than a single temperature at which the theory is critical. This entire low temperature region displays algebraic ordering, i.e., correlation functions that are power laws of the separation between operators. Indeed, one can regard it as a family of conformal field theories (CFTs), since the anomalous dimensions (critical indices) are continuously varying with temperatures. In this sense the XY universality class can be regarded as a two-dimensional (2d) toy model for N = 4 super-Yang-Mills (SYM), which also has continuously varying anomalous dimensions for (composite) operators depending on the gauge coupling g (or more generally, the complexified coupling which incorporates the θ angle).
Another view of the sine-Gordon model is in terms of a description of vortices that form in a 2d superfluid. Here the circulation is quantized in units of h/M A , where h is Planck's constant and M A is the atomic mass. The XY model angular variable θ(x) corresponds to this phase angle of the superfluid wave function.
With the conventions that we use in this report, the Euclidean action for the theory that we study is given by: Here, t is the "stiffness," which is identified with the inverse temperature in the corresponding XY model, while y ≡ g/2t is the fugacity of vortices. The small g behavior of this theory has been known for a long time. For t > 8π (the low temperature regime), the renormalization group (RG) flow (toward the infrared) of g is g → 0. We thus recover the theory with long range correlations (algebraic ordering) in this part of the phase diagram. Thus we see that t > 8π corresponds to the low temperature phase of the XY model, in terms of its long distance behavior. On the other hand, for t < 8π (the high temperature regime), the flow of g is g → ∞, which leads to screening and an absence of criticality. One should not confuse t here with the reduced temperature t red = (T −T BKT )/T BKT . Rather t maps to an inverse temperature in the XY model. Nevertheless the BKT behavior ξ ∼ e a/ √ t red translates in the sine-Gordon theory into ξ ∼ e b/ √ 8π−t for t < 8π. Thus we retain the essential singularity at the critical temperature/stiffness.
This sine-Gordon model falls into a class of solid-on-solid (SOS) models, which display the BKT transition. In the context of the SOS models, this is a roughening transition. The relationship between roughening transitions in crystal facets and the XY universality class was first described in [7,8]. The variable φ(x) is viewed as a height variable above a two-dimensional (2d) surface-the facet of a crystal. Hence above the critical t, where φ(x) becomes highly nonuniform, it is described as the rough surface that grows at high temperatures. The critical line is described by t c (g), where g then plays the role of labeling different types of crystals (to the extent that a one parameter description is usable). However, when y = g/2t 1 the potential term dominates and one is driven to the uniform

Fourier acceleration
Here we describe the Fourier accelerated hybrid Monte Carlo (HMC) algorithm that is used in our simulations [9,10]. The purpose of this modification of HMC is to reduce autocorrelations between configurations of φ(x) that are produced in the course of the simulation. As a result, a shorter simulation can be run, while producing comparable results and uncertainties. Previous work using Fourier accelerated HMC includes [11,12]. The HMC simulation proceeds as usual; the momentum field π(x) is drawn at random from a Gaussian distribution with unit variance. This corresponds to "integrating in" a Gaussian field in the partition function, leading to the "Hamilitonian" where as usual the Laplacian is discretized by In this notation ∂ µ is the forward difference operator in the µ direction, ∂ * µ is the backward difference operator, andμ is a unit vector in the µ direction. We work in lattice units, a = 1. The Hamiltonian H is evaluated to obtain H(0). Next, the fields π(x), φ(x) are evolved according to Hamilton's equationṡ for a trajectory of length τ. The numerical integration technique for this evolution in fictitious time should be reversible and area-preserving, where area refers to the functional integration measure [dπ(x) dφ(x)]. The standard method is the leapfrog algorithm. There is a step size dt in this integrator, and the Hamiltonian H(τ) at the end of the leapfrog trajectory will differ from H(0) due to error associated with not taking the dt → 0 limit. Thus at the end of the trajectory, we apply the Metropolis accept/reject step to obtain a "perfect" algorithm, which will sample the functional integral with the correct weight, the canonical distribution corresponding to H[π(x), φ(x)].
The Fourier acceleration enters into the leapfrog trajectory, where the Fourier modes of φ(x) and π(x) are integrated with a step size An equal number of steps N τ = τ/dt are taken for each mode k. The Fourier transform of the force −∂H/∂φ(x) must also be used in these equations. In (2.5), ∆(k) is the Fourier transform of −∆(x, y): for an L 1 × L 2 lattice. By integrating the longer wavelength modes with a larger step size, they are moved farther in configuration space. This tends to reduce autocorrelations, because it is precisely the long wavelength modes which are the source of this difficulty. We note that it is not necessary to take the k dependent step size dt(k) to be of the form (2.5). A follow-up study to the present one will be to optimize the choice of dt(k) using machine learning techniques. The figure of merit that will be maximized (i.e., the training goal) is the inverse of the integrated autocorrelation time. The present study is partly a preliminary step to this more extensive analysis of Fourier acceleration, using the sine-Gordon model as a working context. This is partly motivated by the fact that two-dimensional real scalar field theories are easily simulated on relatively small scale computing resources, and are thus well-suited to exploratory studies.

Thickness
Following [13] we study the thickness: Here V = L × L is the system size, "volume." The thickness provides a measure of the roughness, on average, for a given parameter pair (t, y). This can be seen from the fact that if the entire lattice sits in a single domain, with small fluctuations, then σ will be small. On the other hand, domain walls will contribute a nonzero result even in the absence of fluctuations, so ground state disorder will be picked up by the thickness observable. In order to not bias toward a particular ground state, unless otherwise stated we begin all of our simulations in studies of thickness with a random start. It is important to understand the autocorrelation for this quantity before taking any averages and estimating uncertainties. This is because we need several autocorrelation times in order to thermalize (i.e., to sufficiently advance the memory kernel beyond the initial conditions), and we need to know the separation between statistically independent samples-where the Monte Carlo simulation is effectively a Markov process. The autocorrelation for any observable O(t), where t here is the simulation time (measured in molecular dynamics time units), is given by and N is the total number of time steps in the simulation. In the present case, In Fig. 3 we show short, long and very long time scales. It can be seen that there is an initial rapid decay, but that a longer time component also contributes. In fact it takes O(10 3 ) or more updates to obtain a completely independent configuration. These results are for y = 0.1 with t values that bracket what will turn out to be the critical temperature, t c ≈ 18. In fact, it can be seen that t ≈ t c yields the longest autocorrelation times, which is to be expected. This is because critical fluctuations, which have very long wavelength, lead to significant slow-down in typical Monte Carlo algorithms. We see that the Fourier acceleration has not been entirely effective in alleviating this critical slowing down. On the other hand, it can be seen that monitoring the autocorrelation can be a surprisingly good way to locate the critical temperature. It is interesting to hold y = g/2t, the fugacity, fixed as we vary t. In Figs. 4 it can be seen that for small values of y and t, the thickness just behaves as a Gaussian variance directly proportional to t. It thus appears that the y coupling has essentially no effect in this regime, other than to determine the slope of the line. At larger values we begin to see nonlinear effects.
Hasenbusch et al. provide a perturbative prediction for the finite size scaling of the thickness above and below the transition temperature [13]. In our conventions it is given by t < t c (2.10) in the large L limit. They have also verified this behavior in Monte Carlo simulations using a cluster algorithm. We conduct the same study, but with the Fourier accelerated HMC. Fig. 5 shows precisely this behavior, with t c ≈ 18. We note that the y → 0 limit (fugacity expansion) predicts t c = 8π ≈ 25, so there is apparently some renormalization of t arising fromŷ = 0.1. If we fit the coefficient c of σ 2 c ln L + const. for the t = 22 data, we obtain c ≈ 6.7 ± 0.2, which is to be compared with the results of [13] which predict t/π = 7.0. We attribute the small discrepancy (1.5σ) to a possible underestimation of errors and a renormalization of t. Once differences in conventions are taken into account, our result of t c ≈ 18 is also in agreement with [13].

Outlook
In a forthcoming full length article [14] we will present further results of our simulations, including clustering of the field φ(x) into domains and the entropy of the system.
Along the line of fixed points that occurs at t > 8π, we will have a conformal field theory describing the infrared physics. Since this is a two-dimensional theory, the conformal group is infinite dimensional, represented in terms of the Virasoro algebra. One of the things that we would like to eventually  It can be seen that for t ≈ t c ≈ 18, the autocorrelation is the greatest, as is to be expected. It can be seen that even after 1,000 updates, autocorrelations are still significant for temperatures close to the critical temperature of t c ≈ 18, whereas for O(3000) time steps all memory of the initial configuration has vanished.
investigate is the emergence of this infinite dimensional symmetry algebra at long distances. Indeed, it is not a simple map to find the field operators of the ultraviolet theory that will correspond to the Virasoro generators. This is because the symmetry is not a property of the ultraviolet theory, due to the nonzero g. Various methods can be employed to improve the simulation. This includes embedding a cluster algorithm such as was done in [13], as well as parallel tempering to allow better sampling of the ground state degeneracy. These advances will appear in our future works.  . Thickness versus t for various y for an L = 64 lattice. It can be seen that a cross-over between different behaviors occurs for y ∼ 1, with a nonlinear dependence of σ 2 on t appearing for large y. On the right, we show thickness versus t. Some interesting features emerge at the smallest values of t seen in this resolution, though at the smallest value of t in the right-hand panel, numerical instabilities may be setting in.