Berezinskii-Kosterlitz-Thouless phase transition from lattice sine-Gordon model
BBerezinskii-Kosterlitz-Thouless phase transition from latticesine-Gordon model
Joel
Giedt ,(cid:63) and James
Flamino Department of Physics, Applied Physics and Astronomy, Rensselaer Polytechnic Institute, 110 8th Street,Troy NY 12180 USA
Abstract.
We obtain nonperturbative results on the sine-Gordon model using the lat-tice field technique. In particular, we employ the Fourier accelerated hybrid Monte Carloalgorithm for our studies. We find the critical temperature of the theory based autocor-relation time, as well as the finite size scaling of the “thickness” observable used in anearlier lattice study by Hasenbusch et al.
Many revolutions in theoretical physics occured in the early 1970s; one of these was the discoverythat phase transitions did not always associate themselves with spontaneous symmetry breaking andlong-range order. The physics of vortices, and the corresponding topological phase transitions, al-lowed 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 modeland the two-dimensional (2d) Coulomb gas, became proving grounds for ideas related to the Wilso-nian renormalization group. We illustrate the binding and unbinding of vortex – anti-vortex pairs oneither side of the transition in Fig. 1.The Berezinskii-Kosterlitz-Thouless (BKT) transition [4, 5] was originally connected to the su-perfluid transition in two dimensions, such as the He thin film. Later it was realized that it couldalso occur in superconducting thin films, in spite of being a charged superfluid. In this case the super-currents screen the fluctuations. However, the screening length is given by
Λ = λ / d , where λ is themagnetic penetration depth and d is the film thickness. Thus for large disorder (hence large λ ) and forsu ffi ciently thin films, Λ can still be large enough for the algebraic order of the XY universality classto emerge. The BKT transition is remarkably di ff erent from the usual Ginzburg-Landau transition. It is atransition without an order parameter. There is a line of conformal fixed points, rather than a singletemperature at which the theory is critical. This entire low temperature region displays algebraicordering, i.e., correlation functions that are power laws of the separation between operators. Indeed, (cid:63) Speaker, e-mail: [email protected]. JG was supported in part by the Department of Energy, O ffi ce of Science, O ffi ce ofHigh Energy Physics, Grant No. DE-SC0013496. Significant parts of this research were done using resources provided by theOpen Science Grid [1, 2], which is supported by the National Science Foundation award 1148698, and the U.S. Department ofEnergy’s O ffi ce of Science. We are also appreciative of XSEDE [3] resources (Stampede), where other significant computationsfor this study were performed. See for instance the review [6]. a r X i v : . [ h e p - l a t ] O c t igure 1. Illustration of the vortex binding andunbinding on the two side of the BKT transition. 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 classcan be regarded as a two-dimensional (2d) toy model for N = 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 a2d 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 thesuperfluid wave function.With the conventions that we use in this report, the Euclidean action for the theory that we studyis given by: S [ φ ] = t (cid:90) d x (cid:40)
12 [ ∂ µ φ ( x )] − g cos φ ( x ) (cid:41) (1.1)Here, t is the “sti ff ness,” which is identified with the inverse temperature in the corresponding XYmodel, while y ≡ g/ t is the fugacity of vortices. The small g behavior of this theory has beenknown for a long time. For t > π (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 > π corresponds to thelow temperature phase of the XY model, in terms of its long distance behavior. On the other hand,for t < π (the high temperature regime), the flow of g is g → ∞ , which leads to screening and anabsence of criticality.One should not confuse t here with the reduced temperature t red = ( T − T BKT ) / T BKT . Rather t mapsto an inverse temperature in the XY model. Nevertheless the BKT behavior ξ ∼ e a / √ t red translates inthe sine-Gordon theory into ξ ∼ e b / √ π − t for t < π . Thus we retain the essential singularity at thecritical temperature / sti ff ness.This sine-Gordon model falls into a class of solid-on-solid (SOS) models, which display the BKTtransition. In the context of the SOS models, this is a roughening transition. The relationship betweenroughening 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 ofa crystal. Hence above the critical t , where φ ( x ) becomes highly nonuniform, it is described as therough surface that grows at high temperatures. The critical line is described by t c ( g ), where g thenplays the role of labeling di ff erent types of crystals (to the extent that a one parameter description isusable). However, when y = g/ t (cid:29) y smooth rough Figure 2.
Sketch of hypothetical phase boundary forsine-Gordon model. vacuum states where φ ( x ) is frozen to a multiple of 2 π . Thus we expect the curve t c ( y ) to tend toinfinity as y is increased, as shown in the sketch, Fig. 2. Here we describe the Fourier accelerated hybrid Monte Carlo (HMC) algorithm that is used in oursimulations [9, 10]. The purpose of this modification of HMC is to reduce autocorrelations betweenconfigurations of φ ( x ) that are produced in the course of the simulation. As a result, a shorter simu-lation can be run, while producing comparable results and uncertainties. Previous work using Fourieraccelerated HMC includes [11, 12].The HMC simulation proceeds as usual; the momentum field π ( x ) is drawn at random from aGaussian distribution with unit variance. This corresponds to “integrating in” a Gaussian field in thepartition function, leading to the “Hamilitonian” H = (cid:88) x (cid:26) π ( x ) − t φ ( x ) ∆ φ ( x ) − g t cos φ ( x ) (cid:27) (2.1)where as usual the Laplacian is discretized by ∆ φ ( x ) = (cid:88) µ = ∂ ∗ µ ∂ µ φ ( x ) = (cid:88) µ = (cid:2) φ ( x + ˆ µ ) + φ ( x − ˆ µ ) − φ ( x ) (cid:3) (2.2)In this notation ∂ µ is the forward di ff erence operator in the µ direction, ∂ ∗ µ is the backward di ff erenceoperator, 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 equations˙ φ ( x ) = ∂ H ∂π ( x ) , ˙ π ( x ) = − ∂ H ∂φ ( x ) (2.3)or a trajectory of length τ . The numerical integration technique for this evolution in fictitious timeshould 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 inte-grator, and the Hamiltonian H ( τ ) at the end of the leapfrog trajectory will di ff er from H (0) due toerror associated with not taking the dt → / reject step P acc = min(1 , e − ∆ H ) , ∆ H = H ( τ ) − H (0) (2.4)to obtain a “perfect” algorithm, which will sample the functional integral with the correct weight, thecanonical 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 dt ( k ) = dt / ( ∆ ( k ) + m ff ) (2.5)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 ): ∆ ( k ) = (cid:88) µ = ( π k µ / L µ ) , k µ = , . . . , L µ − L × L lattice. By integrating the longer wavelength modes with a larger step size, they aremoved farther in configuration space. This tends to reduce autocorrelations, because it is precisely thelong wavelength modes which are the source of this di ffi culty.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 learningtechniques. The figure of merit that will be maximized (i.e., the training goal) is the inverse of theintegrated autocorrelation time. The present study is partly a preliminary step to this more extensiveanalysis of Fourier acceleration, using the sine-Gordon model as a working context. This is partlymotivated by the fact that two-dimensional real scalar field theories are easily simulated on relativelysmall scale computing resources, and are thus well-suited to exploratory studies. Following [13] we study the thickness: σ = V (cid:88) x ,y (cid:104) ( φ ( x ) − φ ( y )) (cid:105) (2.7)Here V = L × L is the system size, “volume.” The thickness provides a measure of the roughness, onaverage, for a given parameter pair ( t , y ). This can be seen from the fact that if the entire lattice sits ina single domain, with small fluctuations, then σ will be small. On the other hand, domain walls willcontribute a nonzero result even in the absence of fluctuations, so ground state disorder will be pickedup by the thickness observable. In order to not bias toward a particular ground state, unless otherwisestated 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 andestimating uncertainties. This is because we need several autocorrelation times in order to thermalize(i.e., to su ffi ciently advance the memory kernel beyond the initial conditions), and we need to knowhe separation between statistically independent samples—where the Monte Carlo simulation is e ff ec-tively a Markov process. The autocorrelation for any observable O ( t ), where t here is the simulationtime (measured in molecular dynamics time units), is given by C ( t ) = N N − t N − t − (cid:88) t (cid:48) = ( O ( t (cid:48) + t ) − (cid:104)O(cid:105) )( O ( t (cid:48) ) − (cid:104)O(cid:105) ) N = N N − (cid:88) t = ( O ( t ) − (cid:104)O(cid:105) ) , (cid:104)O(cid:105) = N N − (cid:88) t = O ( t ) (2.8)and N is the total number of time steps in the simulation. In the present case, O = V (cid:88) x ,y ( φ ( x ) − φ ( y )) (2.9)In Fig. 3 we show short, long and very long time scales. It can be seen that there is an initial rapiddecay, but that a longer time component also contributes. In fact it takes O (10 ) or more updates toobtain a completely independent configuration. These results are for y = . t values that bracketwhat will turn out to be the critical temperature, t c ≈
18. In fact, it can be seen that t ≈ t c yieldsthe longest autocorrelation times, which is to be expected. This is because critical fluctuations, whichhave very long wavelength, lead to significant slow-down in typical Monte Carlo algorithms. We seethat the Fourier acceleration has not been entirely e ff ective in alleviating this critical slowing down.On the other hand, it can be seen that monitoring the autocorrelation can be a surprisingly good wayto locate the critical temperature.It is interesting to hold y = g/ t , the fugacity, fixed as we vary t . In Figs. 4 it can be seen that forsmall values of y and t , the thickness just behaves as a Gaussian variance directly proportional to t . Itthus appears that the y coupling has essentially no e ff ect in this regime, other than to determine theslope of the line. At larger values we begin to see nonlinear e ff ects.Hasenbusch et al. provide a perturbative prediction for the finite size scaling of the thickness aboveand below the transition temperature [13]. In our conventions it is given by σ (cid:39) (cid:40) t π ln L + const. , t > t c const. t < t c (2.10)in the large L limit. They have also verified this behavior in Monte Carlo simulations using a clusteralgorithm. We conduct the same study, but with the Fourier accelerated HMC. Fig. 5 shows preciselythis behavior, with t c ≈
18. We note that the y → t c = π ≈ t arising from ˆ y = .
1. If we fit the coe ffi cient c of σ (cid:39) c ln L + const. for the t =
22 data, we obtain c ≈ . ± .
2, which is to be compared withthe results of [13] which predict t /π = .
0. We attribute the small discrepancy (1 . σ ) to a possibleunderestimation of errors and a renormalization of t . Once di ff erences in conventions are taken intoaccount, our result of t c ≈
18 is also in agreement with [13].
In a forthcoming full length article [14] we will present further results of our simulations, includingclustering of the field φ ( x ) into domains and the entropy of the system.Along the line of fixed points that occurs at t > π , we will have a conformal field theory describ-ing the infrared physics. Since this is a two-dimensional theory, the conformal group is infinite dimen-sional, represented in terms of the Virasoro algebra. One of the things that we would like to eventually C ( τ ) τ t = 14t = 16t = 18t = 20t = 22 -0.2 0 0.2 0.4 0.6 0.8 1 0 200 400 600 800 1000 C ( τ ) τ t = 14t = 16t = 18t = 20t = 22-0.2 0 0.2 0.4 0.6 0.8 1 0 500 1000 1500 2000 2500 3000 C ( τ ) τ t = 14t = 16t = 18t = 20t = 22 Figure 3.
The autocorrelation in the thickness observable on a 64 ×
64 lattice for y = .
1, and t that bracketthe phase transition. The three panels show progressively longer simulation time scales. 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 theVirasoro generators. This is because the symmetry is not a property of the ultraviolet theory, due tothe nonzero g .Various methods can be employed to improve the simulation. This includes embedding a clusteralgorithm such as was done in [13], as well as parallel tempering to allow better sampling of theground state degeneracy. These advances will appear in our future works. σ t y = 0.000244y = 0.000488y = 0.000977y = 0.00195y = 0.00390y = 0.00781y = 0.0156y = 0.0312y = 0.0625y = 0.125y = 0.25y = 0.5y = 1.0y = 2.0y = 4.0y = 8.0y = 16.0y = 32.0y = 64.0 0 200 400 600 800 1000 1200 1400 0 0.5 1 1.5 2 σ t y = 4.0y = 8.0y = 16.0y = 32.0y = 64.0y = 128.0 Figure 4.
Thickness versus t for various y for an L =
64 lattice. It can be seen that a cross-over between di ff erentbehaviors occurs for y ∼
1, with a nonlinear dependence of σ on t appearing for large y . On the right, we showthickness versus t . Some interesting features emerge at the smallest values of t seen in this resolution, though atthe smallest value of t in the right-hand panel, numerical instabilities may be setting in. σ L t = 14t = 16t = 18t = 20t = 22
Figure 5.
The thickness observable transitioningbetween the two behaviors (2.10) at t c ≈ y = . eferences [1] R. Pordes, et al., J. Phys. Conf. Ser. , 012057 (2007)[2] I. Sfiligoi, D.C. Bradley, B. Holzman, P. Mhashilkar, S. Padhi, F. Wurthwein, 2009 WRI WorldCongress on Computer Science and Information Engineering , 428 (2009)[3] J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither, A. Grimshaw, V. Hazlewood, S. Lathrop,D. Lifka, G.D. Peterson et al., Computing in Science & Engineering , 62 (2014)[4] V. Berezinskii, Soviet Physics JETP , 493 (1971)[5] J.M. Kosterlitz, D.J. Thouless, J. Phys. C6 , 1181 (1973)[6] L. Benfatto, C. Castellani, T. Giamarchi, in
40 years of Beresinskii-Kosterlitz-Thouless theory ,edited by J.V. José (World Scientific, Singapore, 2013), pp. 161–199, [7] S. Chui, J. Weeks, Phys. Rev.
B14 , 4978 (1976)[8] S. Chui, J. Weeks, Phys. Rev. Lett. , 733 (1978)[9] S. Duane, A. Kennedy, B. Pendelton, D. Roweth, Phys. Lett. , 216 (1985)[10] A. Ferreira, R. Toral, Phys. Rev. E47 , 3848 (1993)[11] D. Espriu, V. Koulovassilopoulos, A. Travesset, Phys. Rev.
D56 , 6885 (1997), hep-lat/9705027 [12] S. Catterall, S. Karamov, Phys. Lett.
B528 , 301 (2002), hep-lat/0112025 [13] M. Hasenbusch, M. Marcu, K. Pinn, Physica
A211 , 255 (1994), hep-lat/9408005hep-lat/9408005