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 on autocorrelation time, as well as the finite size scaling of the"thickness"observable used in an earlier lattice study by Hasenbusch et al. We study the entropy, which is smooth across all temperatures, supportive of an infinite order transition. This system has a well-known duality with the massive Thirring model, which can play the role of a toy models for Montonen-Olive duality in N=4 super-Yang-Mills theory, since it relates solitons to elementary field excitations. Our research lays a groundwork for such study on the lattice.


I. INTRODUCTION
There were many important advances in theoretical physics in the early 1970s. Among these was the discovery of phase transitions that were not characterized by spontaneous symmetry breaking and long-range order. The dynamics of vortices, and the corresponding topological phase transitions, provided a way to have critical phenomena while still satisfying the theorems forbidding spontaneous symmetry breaking in two dimensions. The XY model and the two-dimensional (2d) Coulomb gas provide a system where this physics of vortices and topological order can be studied. They also played an important role in the development of the Wilsonian renormalization group. The binding and unbinding of vortex -anti-vortex pairs on either side of the transition is shown in Fig. 1.
The Berezinskii-Kosterlitz-Thouless (BKT) transition [1,2] was originally formulated to describe the superfluid transition in two dimensions, such as the 4 He thin film. It was subsequently applied to superconducting thin films, which are a sort of charged superfluid.
As a result, the supercurrents screen the fluctuations so care must be taken in attempting to apply the BKT theory to this type of system. The screening length is given by Λ = λ 2 /d, where λ is the magnetic penetration depth and d is the film thickness. If the disorder is large, then λ is also large. If in addition, the film is very thin, so that d is very small, Λ can be large. Then we can approximately neglect the screening, and BKT can be applied. 1 The BKT transition can be contrasted with the Ginzburg-Landau transition. A key difference is the absence of a conventional order parameter, though in the case of the XY model one can measure vorticity to distinguish the two phases. But one particularly interesting feature of the XY model is that there is a line of conformal fixed points, rather than a single temperature at which the theory is critical. Throughout this low temperature regime there is algebraic ordering of the O(2) spin fields; i.e., correlation functions yield power laws of the separation between operators. Thus one should view this system as a family of conformal field theories (CFTs), since the anomalous dimensions (critical indices) are continuously varying with temperatures. For this reason 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).  Figure taken from [4].
In this article we will study the sine-Gordon model. It is the effective theory of the vortices in the XY model, and has the Euclidean action: 2 Here, t is the "stiffness." Interestingly, it corresponds to the inverse temperature of the corresponding XY model. The quantity y = g/2t is the fugacity of vortices. The small g, or small y, behavior of the sine-Gordon theory is well understood. For t > 8π (the low temperature regime of the XY model), 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. By contrast, for t < 8π (the high temperature regime of the XY model), the flow of g is g → ∞, which leads to screening and an absence of criticality.
Since t > 8π corresponds to the low temperature phase of the XY model, this is where the line of fixed points will occur. It is here that one has continuously varying critical exponents, and conformal field theories will describe the infrared physics. Since this is a two-dimensional theory, the conformal group is infinite dimensional, with generators L m satisfying the well-known Virasoro algebra. Note that this infinite dimensional symmetry algebra emerges at long distances, associated with the flow y → 0 of the renormalization group. Thus it is challenging to obtain the composite field operators of the ultraviolet theory that will correspond to the Virasoro generators. This is because the conformal symmetry is not a feature of the ultraviolet theory, due to the nonzero y.
As mentioned above, the temperature t of the sine-Gordon theory maps to an inverse temperature in the XY model. In the XY model, the correlation length ξ approaches criticality according to ξ ∼ e a/ √ t red , where t red = (T − T BKT )/T BKT is the reduced temperature. In the sine-Gordon theory this becomes ξ ∼ e b/ √ 8π−t for t < 8π, where ξ is the correlation length.
Thus we continue to have an essential singularity in the critical temperature/stiffness. A result of this is that the transition is of infinite order, so that derivatives of the free energy will be smooth functions. We find results consistent with this in our study of the entropy below, in Section VI.
The sine-Gordon model belongs to the family of solid-on-solid (SOS) models which have critical behavior consistent with the BKT transition. When one looks into the origin of the SOS models, one finds a roughening transition. The correspondence between roughening transitions in crystal facets and the XY universality class was first discussed in [5,6]. The variable φ(x) of the sine-Gordon model is interpreted as a height variable above a twodimensional (2d) surface-the facet of a crystal. Below, we display figures that visualize this picture from actual simulations; see Figs. 20 and 21. As a result, above the critical t, where φ(x) becomes highly nonuniform, the sine-Gordon model describes the high temperature growth of a rough surface. There is a critical line t c (y), and the fugacity y labels different types of crystals. When y 1, the cosine potential term in the action (1.1) dominates and one is driven to the uniform ground states where φ(x) is frozen to a multiple of 2π. Thus we expect the curve t c (y) to tend to infinity as y is increased, as shown in the sketch, Fig. 2.
It was mentioned above that the sine-Gordon theory is a toy model for N = 4 super-Yang-Mills because it has a line of fixed points where nontrivial conformal field theories exist.
There is an additional correspondence, namely the duality between solitons and elementary excitations. In the present case, the sine-Gordon theory is known to be dual to the massive Thirring model, with kink-like solitons in the sine-Gordon theory being dual to elementary fermions in the Thirring theory [7]. A similar duality exists in N = 4 super-Yang-Mills: on the Coulomb branch, W-bosons are dual to 't Hooft-Polyakov monopoles, which are type of soliton. Thus by studying duality on the lattice in the context of the sine-Gordon model and Thirring model, we can lay groundwork for the eventual study of electric-magnetic S-duality in N = 4 super-Yang-Mills. This has a practical advantage in that simulating two-dimensional models is far less costly in terms of computer resources.
In this article we explore the behavior of two useful observables in a lattice formulation of the sine-Gordon model. We use the method of Fourier accelerated hybrid Monte Carlo, discussed in the next section. We are able to identify the value of t where the topological transition occurs, and find that it agrees with previous determinations. We also provide some evidence that the transition is of infinite order. Preliminary aspects of a part of this study have previously appeared in [4].
We now summarize the content of this article. In Section II we review Fourier acceleration of hybrid Monte Carlo. In Section III we review the relation between the 2d Coulomb gas and the sine-Gordon model, which also allows us to exhibit the fugacity expansion. This is relevant to various statements in later parts of the paper. In Section IV a few words about the solitons of this model are given. These are again well-known facts. After this we move on to our results, beginning in Section V where we show our main simulation results and the thickness observable. We find results consistent with those obtained in an earlier study by Hasenbusch et al. [8], which used a cluster algorithm. Next in Section VI we show our explorations of entropy, and find a picture that is consistent with an infinite order transition. Section VII takes a closer look at the domains that form when the barrier height is significant. We close the paper with Conclusions and a number of Appendices that review some details related to earlier discussion, including finer points of well-known results.
The appendices also discuss some of the technical points of our implementation, and general musings about the fugacity expansion.

II. FOURIER ACCELERATION
In this section we briefly discuss the Fourier accelerated hybrid Monte Carlo (HMC) algorithm that is the basis of our simulations [9,10]. The point this generalization of HMC is to reduce autocorrelations between configurations of φ(x) that are produced in the course of the simulation. Provided this is successful, a shorter simulation can be used, while producing equivalent results and uncertainties. Previous work using Fourier accelerated HMC includes [11,12].
At the beginning of the HMC simulation, 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 Here, ∂ µ is the forward difference operator in the µ direction, ∂ * µ is the backward difference operator, andμ is a unit vector in the µ direction. As usual, 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 τ , which we typically set to 1. 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, which is what we use. This integrator has step size dt, and the Hamiltonian H(τ ) at the end of the leapfrog trajectory will differ from H(0) due to dt not being zero.
To correct for this, 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 is introduced 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 evolving 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 that lead to long autocorrelation times.
It is certainly not necessary to take the k dependent step size dt(k) to be of the form (2.5). A further study that we will conduct at a future point is to generalize the choice of dt(k) and then optimize it 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.
Our current study that has the feature of being a first step in extensive investigations of the Fourier acceleration technique, using the sine-Gordon model as easily simulated toy model.

III. RELATION TO THE TWO-DIMENSIONAL COULOMB GAS
An important property of the SG model is its relation to the 2d Coulomb gas. It is here that one can understand its dynamics in terms of a collection of charges. Concrete applications include the physics of defects in liquid crystals, vortices in superconducting thin films, and vortices in 2d superfluids, such as a monolayer of 4 He on a substrate. Such features can also be realized in cold atom arrays. It is worth noting that some of the most important properties of BKT physics can be obtained from the multiple perspectives. For instance the renormalization group equations that were originally obtained by Kosterlitz [13] can also be obtained from the sine-Gordon model [14].
The correlation functions that are the focus of this aspect of the sine-Gordon model are not those of the elementary fields φ(x). Rather, we work with what high energy physicists would call "composite operators," e ±iφ(x) . These have an interpretation as vortex and antivortex creation operators.
The basic two-point function to be considered is which corresponds to the amplitude for the creation of a vortex at y and its subsequent annihilation at x. The action that determines this amplitude in the path integral has been given above in Eq. (1.1), leading us to evaluate We proceed by expanding the exponential In the correlation function, only the terms with an equal number of plus signs and minus signs in the exponent will survive. 3 This further implies that only even powers of g will contribute (n even). The interpretation of the nth term in this sum is the virtual pair production of n/2 vortex-antivortex pairs that interact with the vortex at y and antivortex at x as well as with each other. The connection between the 2d Coulomb gas, the interpretation in terms of vortices of the XY model, and the sine-Gordon model has been 3 It turns out that n must be even because of a charge neutrality condition requiring cancelations of signs, due to the logarithmic divergence of the long distance 2d Coulomb potential (see e.g. [15]). explained in many places; see for instance [3,[15][16][17]. One thing that is interesting is that the temperature t of the sine-Gordon (roughening) model maps to an inverse temperature in the corresponding 2d Coulomb gas model. Hence there is a sort of t → 1/t duality between these equivalent descriptions, reminiscent of the Kramers-Wannier duality in the 2d Ising system or the strong-weak coupling duality g → 1/g (S-duality) in N = 4 super-Yang-Mills.
Renormalizing the composite operators 4 we obtain renormalized averages of these operators in the theory with action S 0 (i.e., the From this we find that Here i, j = 0, . . . , n + 1 with x 0 = y, x n+1 = z and 0 = 1, n+1 = −1 so that they are not summed over, hence the primed summation symbol. Z SG is the same expression but without the i, j = 0 and n + 1 terms (vacuum diagrams only). I.e., with i, j = 1, . . . , n only.
Note that there are singularities when i j < 0 and x i → x j under the integration. It is not an integrable singularity for t ≥ 4π. So it is not the case that we have rendered the correlation function finite to all orders simply by renormalizing the composite fields by the factor ζ. All that has been accomplished is to make (3.6) finite, but the subsequent integration of this expression in the sine-Gordon model leads to further difficulties. 6 This 4 See for instance Section 32.1 of Zinn-Justin [15] for further details. Whereas the free theory is obviously finite and free of UV divergences, the same is not true of composite operators in that theory. They can have (trivial) short-distance singularities, as is the case of the operator discussed here. 5 Details of this calculation are given in Appendix A. 6 Although we have emphasized the short distance singularities x i → x j at the beginning of this paragraph, there will also be long distance divergences coming from i j > 0 if we work in infinite 2d space.
is one reason to perform the lattice discretization, because all quantities are under control and finite.
Hard-core repulsion can be inserted into the above formula, thereby avoiding x i = x j .
Of course, this does not follow from the path integral of the sine-Gordon model directly, but is an add-on to control the singularity and model the regulation of this infinity in the physical system. The long distance singularity that occurs from infinite |x i − x j | can also be moderated by periodic boundary conditions or some other finite system size regulation, which is again a matter of inserting physical constraints on the mathematical calculation in order to avoid the bad behavior.

B. Observations about renormalization
It may be that we would like to renormalize the theory, absorbing the UV divergences in the low energy parameters. To accomplish this in renormalized perturbation theory, we would cancel these infinities with counterterms; or, equivalently, we would rescale the bare fields and re-express the bare parameters in terms of long distance quantities. The problem is whether redefinitions of t and g, together with φ → √ Zφ, suffices. After all, classically, the mass dimension of φ is d φ = 1 2 (d − 2) = 0, so we can write down an infinite number of counterterms in the renormalizable Lagrangian. The potential (g/t) cos(φ(x)) is not unique according to this power-counting based around the Gaussian fixed point. However, the field φ(x) will acquire an anomalous mass dimension, as will operators built upon it; so we expect that only a finite number of operators will be relevant or marginal once quantum effects are taken into account. Certainly the observed universality of XY type models suggests this. 7 One constraint on the renormalized theory is the invariance of the bare theory under the shift symmetry φ(x) → φ(x) + 2π. Certainly this forbids a host of terms in the potential.
However, it allows for terms of the form cos[pφ(x) + α] where p is any integer and α is a real number.
The way that we would normally proceed is to enumerate all of the one-particle irreducible (1PI) functions and vertex functions that have a non-negative superficial degrees of divergence. That would tell us how many counterterms we need. But since we are working 7 For instance, including next-to-nearest neighbor couplings will shift the critical value of the XY stiffness K, but not the critical exponents.
in position space (due to the composite operators and cosine potential) and have integrals over x, it is not clear how to do the counting in a perturbation theory structured around the fugacity expansion. Also, this is in terms of composite fields, and we are not working with vertices of elementary fields. In fact, the casual observer will not immediately see a diagramatic expansion in the above formula. So the usual tools seem to be off the table.
We will address this conundrum in analysis that follows.
From the work of [18] we know that in superconducting thin films, the behavior can be modeled as a superposition of BKT fluctuations and Ginzburg-Landau (GL) fluctuations.
In our sine-Gordon model we will be spared of this complication, because it only describes the vortex dynamics of the XY model. However it is interesting to think about how this feature might be incorportated into the effective field theory, and we will have some further comments on this direction of research in the conclusions.

IV. SOLITONS
The minima of the potential V (θ) occur at The solitons occur for static configurations Hence φ(t, x) = ϑ(x). The topologically nontrivial solitons asymptote to different vacua: We are particularly interested in the configurations that are minima of the configuration The solutions are well known; they are given for instance in the wonderful review of topological solitons by 't Hooft, Chapter 4.1 of [19]. For instance, in the case of n − = 0 and n + = 1, the soliton is described by the first branch of Of course in our dynamical configurations this will not be exactly true (until we cool), and so we would like to have a way of measuring the mass m eff of the soliton from such non-extremal configurations. We will explore this in future work.
On the lattice we have an S 1 × S 1 geometry, so the asymptotic behavior described above does not make sense. However, we can still have nontrivial topology by incorporating twisted boundary conditions in the spatial direction x = x 1 : This is also left to future work.

V. SIMULATION RESULTS
As we approach the transition temperature, an increasing number of degrees of freedom participate in the low energy theory, as can be seen by the onset of "algebraic ordering," corresponding to a divergent correlation length. This makes the fugacity expansion unreliable, and in fact is the reason that renormalization group methods were employed to understand the transition. We therefore turn to a numerical study of the model where all features are regulated by the lattice's UV and IR cutoffs.
A. φ distribution and vacuum degeneracy As we have seen in the previous section, the form of the potential V (φ) ∼ cos(φ) indicates that there will be an infinity of degenerate vacua, φ ∼ 2πn, n ∈ Z. We expect to see this feature in the simulations, although the tunneling reflects the algorithm in addition to the physical barrier between these vacua. 8 In Fig. 3 we show the history of the average where V = L 0 L 1 is the simulation volume. Time along the bottom axis is the simulation time, measured in molecular dynamics time units (MDTU). It can be seen that for t = 1.0 and g = 2.0 (fugacity of y = 1.0), no tunneling between the various vacua occurs; here, the simulation was initiated at φ = 0. Clearly the barrier is too large to allow HMC to move into other vacua. In Fig. 4 we reduce the barrier ten-fold, setting g = 0.2 (fugacity of y = 0.1), and now there is significant migration between the vacua. Of course since the barrier height is inversely proportional to t, the effective tunneling rate will also depend on this quantity. In Fig. 5 we show what occurs for three different values of t, all with g = 0.4. As per expectations, tunneling is significantly suppressed by decreasing t. Fig. 6 shows a much larger lattice, 64 × 64 (L = 64), where it is possible that domains are forming, consisting of the different vacua. In this case the different vacuum values of φ would tend to cancel, giving a result that fluctuates about zero. It can be seen in Fig. 7 that domains have not formed, but rather on the larger lattice it is simply stuck near the φ = 0 vacuum, rather than tunneling. Apparently, on the larger volumes, tunneling is more difficult due to a tendency toward coherence across the entire lattice.

B. Thickness
In [8], Hasenbusch et al. introduce the thickness: In this formula, V = L × L is the system size, "volume." They show that this is a useful observable for identifying the critical regime. We will use it here, and find results consistent with theirs. The thickness is a measure of the roughness, on average, for a given parameter pair (t, y). For example, if the entire lattice sits in a single domain, with small fluctuations, 8 In fact, it is an interesting problem to design an HMC-type algorithm that easily moves over these barriers, because of the analogy to topological charge transitions in QCD. In that context, simulations close to the continuum limit with chiral lattice fermions are facing serious problems with non-ergodicity, in that they fail to sample all topological sectors. Perhaps it would be easiest to study this problem in the SG model first. The barrier to tunneling to other vacua is too great for the HMC to allow it. then σ will be small. However, domain walls will contribute a nonzero result even in the absence of fluctuations, so ground state disorder will increase 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, as in Fig. 8, which has φ ∈ [−20, 20] uniformly distributed.
Before taking averages and estimating uncertainties, one should study and characterize the autocorrelation of the thickness observable. This is because we need several autocorrelation times in order to thermalize, 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 MDTUs), is given by and N is the total number of time steps in the simulation, t = 0, 1, . . . , N − 1. In the present case, In Figs. 9, 10 and 11 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. Amusingly, monitoring the autocorrelation can be a surprisingly good way to locate the critical temperature.
We next turn to results where we hold the fugacity y = g/2t fixed as we vary t. In Fig. 12 we show the thickness for y = 0.1, which is near the IR fixed point where y = 0. It can be seen that for these small t values, and at small y, the thickness just behaves as a Gaussian variance directly proportional to t. It thus appears that the y coupling has essentially no In Fig. 15 we see that the approximately linear behavior σ 2 ∼ t has a slope that rises with L, if t is greater than some critical value. In fact it is this transition that will be used to identify our critical temperature in the finite size scaling analysis.
Hasenbusch et al. provide a perturbative prediction for the finite size scaling of the thickness above and below the transition temperature [8]. In our conventions it is given by 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. 16 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 y = 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 t/π = 7.0. We attribute the small discrepancy (1.5σ) to a possible underestimation of errors and a renormalization of t.
There is also an RG interpretation of these results. The dimension of the cosine operator in the action is ∆ = t/4π (see for instance the discussion in Section 4.6 of [20]). Thus for t > 8π, the cosine operator is irrelevant and the theory flows to a free theory, which will be ungapped and thus show the logarithmic divergence in the thickness variable. By contrast for the t < 8π case it will flow to a heavily gapped theory that cuts off sensitivity to the finite size, and thus gives a fixed result for the thickness as a function of L.

VI. ENTROPY
Another important observable is the entropy of the system for different values of the fugacity. The entropy S is calculated using We now show that where S(t) is the lattice action introduced above, which has an explicit factor of 1/t. As a first step we therefore write where nowŜ does not depend on t. It is then clear that Next we must discuss the computation of ln Z(t). With Monte Carlo sampling, we can only do this relative to some fixed value Z(t 0 ). Thus we use 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. Figure previously appeared in [4].
To find d ln Z(t)/dt we use the identity For a select set of t, the action is calculated for 50,000 configurations, where the data starts to be recorded after 10,000 steps. This allows the system to thermalize. We perform this calculation for a wide range of t and fit S(t) /t to a high order polynomial. The polynomial is then integrated to obtain ln[Z(t)/Z(t 0 )]. For the reference value t 0 , we choose t 0 = 1.
An example of the entropy analysis can be seen in Fig. 17 where entropy is calculated for a temperature range of [0, 60] for fugacity y = 0.1. There appears to be a smooth logarithmic increase in entropy with respect to t. This makes sense, as a BKT transition is a transition of infinite order: a first derivative of the free energy is not expected to show a discontinuity, and in fact should be smooth. To confirm this expectation, in Fig. 18 another example is shown for a higher fugacity, y = 1. Due to the higher potential barrier, tunneling was a lot more erratic across temperatures, so the action required a higher order polynomial to reliably model the free energy. And as can be seen in the figure, the resulting curve is smooth with a monotonically increasing entropy.
It is important to consider the effects of autocorrelation in this calculation of entropy as well. The observable is highly local, since it is a sum over the action density. The result is that the autocorrelation is significantly smaller than in the case of the thickness measurements, which involve correlations over large scales. In order to quantify this smaller autocorrelation, in Fig. 19 we show the action for y = 1 generated for two separate random seeds, which alters the initialization of φ(x) at the beginning of the simulation. The difference between the action on the two sets of configurations is well within the error bars, which take into account integrated autocorrelation time in the usual way. As will be seen in subsequent figures, a more interesting behavior occurs for y > 1.

VII. CLUSTERING
We have analyzed the clustering of the field φ(x) into values around 2nπ, n ∈ Z. The method for identifying a cluster is described in Appendix D. In Fig. 20 the various clusters found by this algorithm are shown from a perspective where they are all visible. Fig. 20 shows the same collection of clusters edge-on, where it is apparent that they are grouped into layers. In order to produce results with clearly defined domains, we found that it was helpful to use a random start with φ(x) drawn uniformly from [−20, 20] and a very large fugacity of y = 15. We also set t = 0.1 for this simulation, which tends to freeze the configuration into domains (very low temperature). The correlation between cluster "width" W and population P (the number of sites contained in the cluster) is show in Fig. 22, where 50 lattice configurations were combined. It shows an approximate behavior of W ∼ (ln P ) 2 (7.1) Note that the largest population and width is limited by the fact that this is an L = 64 lattice. We are not aware of a theoretical argument explaining the behavior (7.1), but leave  it as a question for future investigations.

VIII. CONCLUSIONS
One interesting avenue to explore is the possible existence of three phases in the XY type systems, as shown in Fig. 24. Arguments for the existence of this more complicated set of transitions were advanced in classic papers by Halperin, Nelson and Young [21][22][23][24].
Although this feature has been argued for theoretically, and observed experimentally, it has still not been detected in numerical simulations to date. Thus, future work would be to use the more powerful and sophisticated simulation methods that are begun here to observe these types of behaviors. An important question is: "How would they be realized in the sine- or "nodes," that have been identified as belonging to the cluster. Naturally the wider the cluster, the more nodes that it contains. However, the precise relationship is not so easily predicted, and empirically is logarithmic, according to Eq. (7.1).
these ideas. Currently we are developing numerical simulations that will analyze such layered systems, especially in the limit of weak interlayer interactions. Thus we finally obtain which agrees with ZJ (32.6) [15].
Next we consider the correlation functions of composite operators as appears in (3.6), making use of (A6) to prove this identity. The trick (e.g. [15]) is to compute the path by rewriting it in terms of a "source" This can be computed, as usual, by completing the square, leading to To get closer to (3.6) we separate out factors where various factors of 2 have been accounted for in the exponents; in particular the doublecounting of pairs in ij versus i<j . The δ → 0 limit occurs because of the coincident "diagonal" terms x i = x j in ij , and will in practice be regulated by lattice expressions for the propagator. I.e., they correspond to wherek µ = 2 sin(k µ /2). Because of the compact domain of integration (finite lattice spacing), there is no singularity associated x i → x j . The second factor of (A12) agrees with (3.6), in terms of the dependence on |x i − x j |; everything else is swept into the renormalization factor ζ, while introducing the renormalization scale µ.

Appendix B: Physics of vortices
The sine-Gordon model is supposed to capture the physics of vortices, where the field variable φ(x) represents the angle of the U (1) O(2) phase of the order parameter: for instance the condensate in a thin film superconductor, or the phase of the wavefunction in a superfluid. In a superconductor the vortices describe how magnetic field lines penetrate the bulk, and surrounding each vortex is an eddy current. If one considers two vortices such as in Fig. 25, one can see that the eddy currents will tend to repel since the portions of the currents that are closest are going in opposite directions. 9 9 Recall that unlike charges, for currents opposites repel.
requires a predefined cluster count, which is not easily recognizable for large lattice sizes; and due to the nature of the field-versus-site scattering, the commonly-used elbow method is not inherently helpful. We therefore used an alternative: density-based cluster definitions.
Here the cluster count is defined using the DBSCAN method, following the format of using the k − dist values to define cluster radius and minimum point threshold. Once k number of cluster radii (called Eps) are identified, the algorithm uses K-means to classify those clusters. To aid in speeding up the process, centroids c 1 ...c k are initialized within the local space of each Eps j , automatically assigning identified core points and border points to respective cluster c j . For any points missed, we simply follow argmin D(x i , c j ), where D is distance, assigning the noise point x i to closest cluster c j .
The algorithm is quick and efficient, but is still not powerful enough when applied to even larger lattice sizes coupled with small deviations of φ. Many configurations of this form can be too complex for this algorithm. A work-around for this is to analyze extremely large amounts of sites is to take the configuration data and project it onto a 2d plane using principle component analysis (PCA). In PCA, a system with a certain amount of unlabeled data is selected to form a s-dimensional sample space. The covariance matrix is then computed as Σ =         σ 1,1 σ 1,2 · · · σ 1,j σ 2,1 σ 2,2 · · · σ 2,j . . . . . . . . . . . .
Where σ i,i is the variance and σ i,j is the covariance. The eigenvalues and eigenvectors are calculated from Σ, giving a set of potential axes (eigenvectors) that the data could be In terms of our data, we considered the sample space of the configurations as a fibre bundle instead of just using the singular φ data generated by the simulation. A fibre bundle is, for our simple theory, the product of the base space and the target space. Considering the analysis in terms of a fibre bundle allows us to observe our PCA sample space with three features (x, y, and φ) instead of just φ. This gives depth to our analysis, allowing for the cluster algorithm to identify the clusters quickly for a large number of increasingly large lattice sizes without any loss in accuracy. Fig. 26 shows the PCA of a configuration with a lattice size of 256. Computation time is cut down by a factor of two and the cluster count and cluster width is not compromised.
this theory is dual to in the absence of the hard-core repulsion? One path to formulating an answer is to represent the Heaviside unit step function as where as usual > 0 is an infinitesmal number that provides a pole prescription. Thus we obtain the expression If we define ω ij = ( i j t/2π)ω ij (E5) then the expression becomes more amenable to intepretation: Thus we posit a differential operator D(ω) that has Green's function G(ω; x, y) = iω|x i − x j | + ln |x i − x j | (E7) Obviously it is some generalization of the 2d Laplacian.