LLattice sine-Gordon model
James Flamino ∗ and Joel Giedt † Department of Physics, Applied Physics and Astronomy,Rensselaer Polytechnic Institute, 110 8th Street, Troy NY 12065 USA (Dated: March 16, 2018)
Abstract
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. Wefind the critical temperature of the theory based on autocorrelation time, as well as the finite sizescaling of the “thickness” observable used in an earlier lattice study by Hasenbusch et al. We studythe 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 ofa toy models for Montonen-Olive duality in N = 4 super-Yang-Mills theory, since it relates solitonsto elementary field excitations. Our research lays a groundwork for such study on the lattice. ∗ fl[email protected] † [email protected] a r X i v : . [ h e p - l a t ] N ov . INTRODUCTION There were many important advances in theoretical physics in the early 1970s. Amongthese was the discovery of phase transitions that were not characterized by spontaneoussymmetry breaking and long-range order. The dynamics of vortices, and the correspondingtopological phase transitions, provided a way to have critical phenomena while still satisfyingthe theorems forbidding spontaneous symmetry breaking in two dimensions. The XY modeland the two-dimensional (2d) Coulomb gas provide a system where this physics of vorticesand topological order can be studied. They also played an important role in the developmentof the Wilsonian renormalization group. The binding and unbinding of vortex – anti-vortexpairs on either side of the transition is shown in Fig. 1.The Berezinskii-Kosterlitz-Thouless (BKT) transition [1, 2] was originally formulated todescribe the superfluid transition in two dimensions, such as the He thin film. It wassubsequently 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 attemptingto apply the BKT theory to this type of system. The screening length is given by Λ = λ /d ,where λ is the magnetic penetration depth and d is the film thickness. If the disorder islarge, then λ is also large. If in addition, the film is very thin, so that d is very small, Λ canbe large. Then we can approximately neglect the screening, and BKT can be applied. The BKT transition can be contrasted with the Ginzburg-Landau transition. A key dif-ference is the absence of a conventional order parameter, though in the case of the XY modelone can measure vorticity to distinguish the two phases. But one particularly interestingfeature of the XY model is that there is a line of conformal fixed points, rather than a singletemperature at which the theory is critical. Throughout this low temperature regime thereis algebraic ordering of the O (2) spin fields; i.e., correlation functions yield power laws of theseparation between operators. Thus one should view this system as a family of conformalfield theories (CFTs), since the anomalous dimensions (critical indices) are continuouslyvarying with temperatures. For this reason the XY universality class can be regarded as atwo-dimensional (2d) toy model for N = 4 super-Yang-Mills (SYM), which also has con-tinuously varying anomalous dimensions for (composite) operators depending on the gaugecoupling g (or more generally, the complexified coupling which incorporates the θ angle). See for instance the review [3]. IG. 1. Vortex binding and unbinding on thetwo side of the BKT transition. Figure takenfrom [4].
In this article we will study the sine-Gordon model. It is the effective theory of thevortices in the XY model, and has the Euclidean action: S [ φ ] = 1 t (cid:90) d x (cid:26)
12 [ ∂ µ φ ( x )] − g cos φ ( x ) (cid:27) (1.1)Here, t is the “stiffness.” Interestingly, it corresponds to the inverse temperature of thecorresponding XY model. The quantity y = g/ t is the fugacity of vortices. The small g ,or small y , behavior of the sine-Gordon theory is well understood. For t > π (the lowtemperature regime of the XY model), the renormalization group (RG) flow (toward theinfrared) of g is g →
0. We thus recover the theory with long range correlations (algebraicordering) in this part of the phase diagram. By contrast, for t < π (the high temperatureregime of the XY model), the flow of g is g → ∞ , which leads to screening and an absenceof criticality.Since t > π corresponds to the low temperature phase of the XY model, this is wherethe line of fixed points will occur. It is here that one has continuously varying criticalexponents, and conformal field theories will describe the infrared physics. Since this is atwo-dimensional theory, the conformal group is infinite dimensional, with generators L m satisfying the well-known Virasoro algebra. Note that this infinite dimensional symmetryalgebra emerges at long distances, associated with the flow y → y .As mentioned above, the temperature t of the sine-Gordon theory maps to an inverse tem-perature in the XY model. In the XY model, the correlation length ξ approaches criticalityaccording to ξ ∼ e a/ √ t red , where t red = ( T − T BKT ) /T BKT is the reduced temperature. In the Other parameterizations are discussed in Appendix C. mooth rough y t FIG. 2. Sketch of hypothetical phase boundary for sine-Gordon model. sine-Gordon theory this becomes ξ ∼ e b/ √ π − t for t < π , where ξ is the correlation length.Thus we continue to have an essential singularity in the critical temperature/stiffness. Aresult of this is that the transition is of infinite order, so that derivatives of the free energywill be smooth functions. We find results consistent with this in our study of the entropybelow, in Section VI.The sine-Gordon model belongs to the family of solid-on-solid (SOS) models which havecritical behavior consistent with the BKT transition. When one looks into the origin of theSOS models, one finds a roughening transition. The correspondence between rougheningtransitions in crystal facets and the XY universality class was first discussed in [5, 6]. Thevariable φ ( x ) of the sine-Gordon model is interpreted as a height variable above a two-dimensional (2d) surface—the facet of a crystal. Below, we display figures that visualize thispicture 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 temperaturegrowth of a rough surface. There is a critical line t c ( y ), and the fugacity y labels differenttypes of crystals. When y (cid:29)
1, the cosine potential term in the action (1.1) dominates andone is driven to the uniform ground states where φ ( x ) is frozen to a multiple of 2 π . Thuswe 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-4ang-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 elementaryexcitations. In the present case, the sine-Gordon theory is known to be dual to the massiveThirring model, with kink-like solitons in the sine-Gordon theory being dual to elementaryfermions in the Thirring theory [7]. A similar duality exists in N = 4 super-Yang-Mills: onthe Coulomb branch, W-bosons are dual to ’t Hooft-Polyakov monopoles, which are typeof soliton. Thus by studying duality on the lattice in the context of the sine-Gordon modeland Thirring model, we can lay groundwork for the eventual study of electric-magneticS-duality in N = 4 super-Yang-Mills. This has a practical advantage in that simulatingtwo-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 formulationof 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 topologicaltransition occurs, and find that it agrees with previous determinations. We also providesome evidence that the transition is of infinite order. Preliminary aspects of a part of thisstudy have previously appeared in [4].We now summarize the content of this article. In Section II we review Fourier accelerationof hybrid Monte Carlo. In Section III we review the relation between the 2d Coulomb gasand the sine-Gordon model, which also allows us to exhibit the fugacity expansion. Thisis relevant to various statements in later parts of the paper. In Section IV a few wordsabout the solitons of this model are given. These are again well-known facts. After thiswe move on to our results, beginning in Section V where we show our main simulationresults and the thickness observable. We find results consistent with those obtained in anearlier study by Hasenbusch et al. [8], which used a cluster algorithm. Next in Section VIwe show our explorations of entropy, and find a picture that is consistent with an infiniteorder transition. Section VII takes a closer look at the domains that form when the barrierheight is significant. We close the paper with Conclusions and a number of Appendices thatreview 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 generalmusings about the fugacity expansion. 5 I. 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 ofHMC is to reduce autocorrelations between configurations of φ ( x ) that are produced in thecourse of the simulation. Provided this is successful, a shorter simulation can be used, whileproducing equivalent results and uncertainties. Previous work using Fourier acceleratedHMC includes [11, 12].At the beginning of the HMC simulation, the momentum field π ( x ) is drawn at randomfrom a Gaussian distribution with unit variance. This corresponds to “integrating in” aGaussian field in the partition function, leading to the “Hamilitonian” H = (cid:88) x (cid:26) π ( x ) − t φ ( x )∆ φ ( x ) − gt cos φ ( x ) (cid:27) (2.1)where as usual the Laplacian is discretized by∆ φ ( x ) = (cid:88) µ =1 ∂ ∗ µ ∂ µ φ ( x ) = (cid:88) µ =1 [ φ ( x + ˆ µ ) + φ ( x − ˆ µ ) − φ ( x )] (2.2)Here, ∂ µ is the forward difference operator in the µ direction, ∂ ∗ µ is the backward differenceoperator, 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 evolvedaccording to Hamilton’s equations˙ φ ( x ) = ∂H∂π ( x ) , ˙ π ( x ) = − ∂H∂φ ( x ) (2.3)for a trajectory of length τ , which we typically set to 1. The numerical integration techniquefor this evolution in fictitious time should be reversible and area-preserving, where area refersto the functional integration measure [ dπ ( x ) dφ ( x )]. The standard method is the leapfrogalgorithm, 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 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 correctweight, the canonical distribution corresponding to H [ π ( x ) , φ ( x )].6he Fourier acceleration is introduced into the leapfrog trajectory, where the Fouriermodes of φ ( x ) and π ( x ) are integrated with a step size dt ( k ) = dt/ (∆( k ) + m ) (2.5)An equal number of steps N τ = τ /dt are taken for each mode k . The Fourier transform ofthe force − ∂H/∂φ ( x ) must also be used in these equations. In (2.5), ∆( k ) is the Fouriertransform of − ∆( x, y ):∆( k ) = (cid:88) µ =1 ( πk µ /L µ ) , k µ = 0 , . . . , L µ − L × L lattice. By evolving the longer wavelength modes with a larger step size, theyare moved farther in configuration space. This tends to reduce autocorrelations, because itis 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 willbe 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 theFourier 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 ishere that one can understand its dynamics in terms of a collection of charges. Concreteapplications include the physics of defects in liquid crystals, vortices in superconductingthin films, and vortices in 2d superfluids, such as a monolayer of He on a substrate. Suchfeatures can also be realized in cold atom arrays. It is worth noting that some of the mostimportant properties of BKT physics can be obtained from the multiple perspectives. Forinstance 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 arenot those of the elementary fields φ ( x ). Rather, we work with what high energy physi-cists would call “composite operators,” e ± iφ ( x ) . These have an interpretation as vortex andantivortex creation operators. 7 . Two-point function The basic two-point function to be considered is G ( x − y ) = (cid:104) e iφ ( x ) e − iφ ( y ) (cid:105) (3.1)which corresponds to the amplitude for the creation of a vortex at y and its subsequentannihilation at x . The action that determines this amplitude in the path integral has beengiven above in Eq. (1.1), leading us to evaluate G ( x ) = Z − (cid:90) [ dφ ( z )] e − S ( φ ) e i [ φ ( x ) − φ (0)] e ( g/t ) (cid:82) d y cos φ ( y ) Z = (cid:90) [ dφ ( z )] e − S ( φ ) e ( g/t ) (cid:82) d y cos φ ( y ) S ( φ ) = 12 t (cid:90) d x [ ∂ µ φ ( x )] (3.2)We proceed by expanding the exponential e ( g/t ) (cid:82) d y cos φ ( y ) in powers of g/t : ∞ (cid:88) n =0 n ! (cid:18) ( g/t ) (cid:90) d y cos φ ( y ) (cid:19) n = 1 + g t (cid:90) d y (cid:0) e iφ ( y ) + e − iφ ( y ) (cid:1) + g t (cid:90) d y d y (cid:0) e i [ φ ( y )+ φ ( y )] + e − i [ φ ( y )+ φ ( y )] + e i [ φ ( y ) − φ ( y )] + e − i [ φ ( y ) − φ ( y )] (cid:1) + . . . (3.3)In the correlation function, only the terms with an equal number of plus signs and minussigns in the exponent will survive. This further implies that only even powers of g willcontribute ( n even). The interpretation of the n th term in this sum n (cid:89) i =1 e i(cid:15) i φ ( y i ) , (cid:15) i = ± , n (cid:88) i =1 (cid:15) i = 0 , n ∈ Z (3.4)is the virtual pair production of n/ 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 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]). t of the sine-Gordon (roughening) model maps to an inverse temperature inthe corresponding 2d Coulomb gas model. Hence there is a sort of t → /t duality betweenthese equivalent descriptions, reminiscent of the Kramers-Wannier duality in the 2d Isingsystem or the strong-weak coupling duality g → /g (S-duality) in N = 4 super-Yang-Mills.Renormalizing the composite operators e i(cid:15) i θ ( x ) = ζ (cid:2) e i(cid:15) i θ ( x ) (cid:3) ren. , ζ = (Λ /µ ) − t/ π (3.5)we obtain renormalized averages of these operators in the theory with action S (i.e., thefree theory) (cid:42) n (cid:89) i =1 e i(cid:15) i φ ( x i ) (cid:43) ren. = (cid:89) i
9s one reason to perform the lattice discretization, because all quantities are under controland 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 thephysical system. The long distance singularity that occurs from infinite | x i − x j | can alsobe 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 inorder to avoid the bad behavior. B. Observations about renormalization
It may be that we would like to renormalize the theory, absorbing the UV divergencesin the low energy parameters. To accomplish this in renormalized perturbation theory, wewould cancel these infinities with counterterms; or, equivalently, we would rescale the barefields and re-express the bare parameters in terms of long distance quantities. The problemis whether redefinitions of t and g , together with φ → √ Zφ , suffices. After all, classically,the mass dimension of φ is d φ = ( d −
2) = 0, so we can write down an infinite number ofcounterterms in the renormalizable Lagrangian. The potential ( g/t ) cos( φ ( x )) is not uniqueaccording 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 expectthat only a finite number of operators will be relevant or marginal once quantum effects aretaken into account. Certainly the observed universality of XY type models suggests this. One constraint on the renormalized theory is the invariance of the bare theory under theshift 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 realnumber.The way that we would normally proceed is to enumerate all of the one-particle irre-ducible (1PI) functions and vertex functions that have a non-negative superficial degrees ofdivergence. That would tell us how many counterterms we need. But since we are working For instance, including next-to-nearest neighbor couplings will shift the critical value of the XY stiffness K , but not the critical exponents.
10n position space (due to the composite operators and cosine potential) and have integralsover x , it is not clear how to do the counting in a perturbation theory structured aroundthe fugacity expansion. Also, this is in terms of composite fields, and we are not workingwith vertices of elementary fields. In fact, the casual observer will not immediately see adiagramatic 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 bemodeled 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 describesthe vortex dynamics of the XY model. However it is interesting to think about how thisfeature might be incorportated into the effective field theory, and we will have some furthercomments on this direction of research in the conclusions. IV. SOLITONS
The minima of the potential V ( θ ) occur at φ ( x ) = 2 πn, n ∈ Z (4.1)The solitons occur for static configurations ∂φ ( t, x ) ∂t = 0 (4.2)Hence φ ( t, x ) = ϑ ( x ). The topologically nontrivial solitons asymptote to different vacua:lim x →−∞ ϑ ( x ) → πn − , lim x →∞ ϑ ( x ) → πn + , n − (cid:54) = n + (4.3)We are particularly interested in the configurations that are minima of the configurationenergy E = 1 t (cid:90) ∞−∞ dx (cid:26)
12 ( ∂ x ϑ ) + V ( ϑ ) (cid:27) (4.4)subject to the boundary conditions (4.3). These are solutions to the static equations ofmotion (saddle point equation in this Euclidean formulation) δEδϑ ( x ) = 0 ⇒ ∂ x ϑ = ∂V∂ϑ (4.5)11he solutions are well known; they are given for instance in the wonderful review of topo-logical 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 ϑ ( x ) = 4 arctan e mx , m = gt (4.6)Of course in our dynamical configurations this will not be exactly true (until we cool), and sowe would like to have a way of measuring the mass m eff of the soliton from such non-extremalconfigurations. We will explore this in future work.On the lattice we have an S × S geometry, so the asymptotic behavior described abovedoes not make sense. However, we can still have nontrivial topology by incorporating twistedboundary conditions in the spatial direction x = x : θ ( t, x + N x a ) = θ ( t, x ) + 2 πn, n ∈ Z (4.7)This is also left to future work. V. SIMULATION RESULTS
As we approach the transition temperature, an increasing number of degrees of freedomparticipate in the low energy theory, as can be seen by the onset of “algebraic ordering,” cor-responding to a divergent correlation length. This makes the fugacity expansion unreliable,and in fact is the reason that renormalization group methods were employed to understandthe transition. We therefore turn to a numerical study of the model where all features areregulated 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( φ ) indicatesthat there will be an infinity of degenerate vacua, φ ∼ πn , n ∈ Z . We expect to see thisfeature in the simulations, although the tunneling reflects the algorithm in addition to the12hysical barrier between these vacua. In Fig. 3 we show the history of the average φ avg = 1 V (cid:88) x φ ( x ) (5.1)where V = L L is the simulation volume. Time along the bottom axis is the simulationtime, measured in molecular dynamics time units (MDTU). It can be seen that for t = 1 . g = 2 . y = 1 . φ = 0. Clearly the barrier is too large to allow HMC tomove into other vacua. In Fig. 4 we reduce the barrier ten-fold, setting g = 0 . y = 0 . t , the effective tunneling rate will also dependon 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. 6shows 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 tendto cancel, giving a result that fluctuates about zero. It can be seen in Fig. 7 that domainshave 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 toa tendency toward coherence across the entire lattice. B. Thickness
In [8], Hasenbusch et al. introduce the thickness: σ = 1 V (cid:88) x,y (cid:104) ( φ ( x ) − φ ( y )) (cid:105) (5.2)In this formula, V = L × L is the system size, “volume.” They show that this is a usefulobservable for identifying the critical regime. We will use it here, and find results consistentwith theirs. The thickness is a measure of the roughness, on average, for a given parameterpair ( t, y ). For example, if the entire lattice sits in a single domain, with small fluctuations, 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 thecontinuum limit with chiral lattice fermions are facing serious problems with non-ergodicity, in that theyfail to sample all topological sectors. Perhaps it would be easiest to study this problem in the SG modelfirst. φ a v g MDTU t=1.0, g=2.0
FIG. 3. The average value of φ on a 4 × t = 1 . g = 2 .
0. It can be seen that the field fluctuates about the vacuum it was started in, φ ( x ) = 0 ∀ x .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 theabsence 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 beginall of our simulations in studies of thickness with a random start, as in Fig. 8, which has φ ∈ [ − ,
20] uniformly distributed.Before taking averages and estimating uncertainties, one should study and characterizethe autocorrelation of the thickness observable. This is because we need several autocorrela-tion times in order to thermalize, and we need to know the separation between statisticallyindependent 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 (measured14 φ a v g MDTU t=1.0, g=0.2
FIG. 4. The average value of φ on a 4 × t = 1 . g = 0 .
2. I.e., a barrier height 10 times lower than in Fig. 3. It can be seen that the field nowtunnels between vacua φ ( x ) = 2 πn ∀ x , n ∈ Z . With g = 0 .
2, the barrier to tunneling to othervacua is low enough for the Fourier accelerated HMC to allow it. Because the lattice is small, theentire field configuration moves coherently between the vacua. in MDTUs), is given by C ( t ) = 1 N N − t N − t − (cid:88) t (cid:48) =0 ( O ( t (cid:48) + t ) − (cid:104)O(cid:105) )( O ( t (cid:48) ) − (cid:104)O(cid:105) ) N = 1 N N − (cid:88) t =0 ( O ( t ) − (cid:104)O(cid:105) ) , (cid:104)O(cid:105) = 1 N N − (cid:88) t =0 O ( t ) (5.3)and N is the total number of time steps in the simulation, t = 0 , , . . . , N −
1. In the presentcase, O = 1 V (cid:88) x,y ( φ ( x ) − φ ( y )) (5.4)In Figs. 9, 10 and 11 we show short, long and very long time scales. It can be seen thatthere is an initial rapid decay, but that a longer time component also contributes. In fact15 φ a v g MDTU t=2.0, g=0.4t=4.0, g=0.4t=8.0, g=0.4
FIG. 5. The average value of φ on a 4 × t and g = 0 .
4. It can be seen that at smaller t the field rarely tunnels between vacua φ ( x ) = 2 πn ∀ x , n ∈ Z , whereas for larger t the tunneling is fairly rapid. This is because the heightof the barrier is proportional to 1 /t and the strength of fluctuations is proportional to √ t . it takes O (10 ) or more updates to obtain a completely independent configuration. Theseresults are for y = 0 . t values that bracket what will turn out to be the criticaltemperature, t c ≈
18. In fact, it can be seen that t ≈ t c yields the longest autocorrelationtimes, which is to be expected. This is because critical fluctuations, which have very longwavelength, lead to significant slow-down in typical Monte Carlo algorithms. We see thatthe Fourier acceleration has not been entirely effective in alleviating this critical slowingdown. Amusingly, monitoring the autocorrelation can be a surprisingly good way to locatethe critical temperature.We next turn to results where we hold the fugacity y = g/ t fixed as we vary t . In Fig. 12we show the thickness for y = 0 .
1, which is near the IR fixed point where y = 0. It can beseen that for these small t values, and at small y , the thickness just behaves as a Gaussianvariance directly proportional to t . It thus appears that the y coupling has essentially no16 φ a v g MDTU L=64, t=1.0, g=0.1
FIG. 6. The average value of φ on a 64 ×
64 lattice as a function of simulation time. Although thecouplings are the same as in Fig. 4, the lattice has 256 times more sites, opening up the possibilityof forming domains of the different vacua. If that was the case, they would tend to cancel in theaverage of φ , leading to values close to zero. Thus a finer-grained study of the lattice configurationsis necessary in order to resolve what is going on. In fact, Fig. 7 shows that this is not what isoccuring. Rather, the entire lattice collapses almost immediately into a single domain, and remainsthere, in contrast to what happens on the 4 × effect in this regime, other than to determine the slope of the line. A further ellaborationof the behavior of the thickness as a function of t and y is given in Figs. 13 and 14.In Fig. 15 we see that the approximately linear behavior σ ∼ t has a slope that riseswith L , if t is greater than some critical value. In fact it is this transition that will be usedto identify our critical temperature in the finite size scaling analysis.Hasenbusch et al. provide a perturbative prediction for the finite size scaling of the thick-ness above and below the transition temperature [8]. In our conventions it is given by σ (cid:39) tπ ln L + const. , t > t c const. t < t c (5.5)17 F r e q u e n c y L=64, y=0.10, tcoupl=1.00, meff=0.50, cfg=10000
FIG. 7. Histogram of the φ value for a simulation with lattice size L = 64, temperature t = 1 . g = 0 . y = 0 . in the large L limit. They have also verified this behavior in Monte Carlo simulations usinga 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 → t c = 8 π ≈
25, so there is apparently some renormalization of t arisingfrom y = 0 .
1. If we fit the coefficient c of σ (cid:39) c ln L + const. for the t = 22 data, we obtain c ≈ . ± .
2, which is to be compared with t/π = 7 .
0. We attribute the small discrepancy(1 . σ ) 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 operatorin the action is ∆ = t/ π (see for instance the discussion in Section 4.6 of [20]). Thus for t > π , the cosine operator is irrelevant and the theory flows to a free theory, which will beungapped and thus show the logarithmic divergence in the thickness variable. By contrastfor the t < π case it will flow to a heavily gapped theory that cuts off sensitivity to the18 t h i ck ne ss MDTU L=64, t=0.1, g=3.0L=64, t=1.0, g=1.0
FIG. 8. The thickness on a 64 ×
64 for two different choices of t and g , as function of simulationtime, beginning with the random start described in the text. For the case of a steep potential, y = g/ t = 15, the fields rapidly settle into minima and the thickness becomes small, as it onlygets contributions from domain walls between the minima. Note however, that ultimately thethickness is largers because that corresponds to small t = 0 .
1, whereas the t = 1 . 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 thefugacity. The entropy S is calculated using S = − (cid:18) ∂F∂t (cid:19) y , F = − t ln Z (6.1)We now show that S ( t ) = ln Z ( t ) + (cid:104) S ( t ) (cid:105) (6.2)19 C ( τ ) τ t = 14t = 16t = 18t = 20t = 22 FIG. 9. The autocorrelation in the thickness observable on a 64 ×
64 lattice for y = 0 .
1, and t that bracket the phase transition. This figure shows short time scales, with accompanying figuresshowing longer time scales. It can be seen that for t ≈ t c ≈
18, the autocorrelation is the greatest,as is to be expected. Figure previously appeared in [4]. where S ( t ) is the lattice action introduced above, which has an explicit factor of 1 /t . As afirst step we therefore write S ( t ) = 1 t ˆ S (6.3)where now ˆ S does not depend on t . It is then clear that S = ln Z + tZ ∂∂t (cid:90) [ dφ ] e − t ˆ S = ln Z + 1 Z (cid:90) [ dφ ] 1 t ˆ Se − S = ln Z + (cid:104) S (cid:105) (6.4)Next we must discuss the computation of ln Z ( t ). With Monte Carlo sampling, we canonly do this relative to some fixed value Z ( t ). Thus we use (cid:90) tt ddt ln Z ( t ) dt = ln Z ( t ) Z ( t ) (6.5)20 C ( τ ) τ t = 14t = 16t = 18t = 20t = 22 FIG. 10. The autocorrelation in the thickness observable, showing longer time scales than Fig. 9.It can be seen that even after 1,000 updates, autocorrelations are still significant for temperaturesclose to the critical temperature of t c ≈
18. Figure previously appeared in [4].
To find d ln Z ( t ) /dt we use the identity ddt ln Z ( t ) = 1 t (cid:104) S ( t ) (cid:105) (6.6)from above. Thus, we can obtain ln[ Z ( t ) /Z ( t )] from a straightforward expectation value—exactly what is easily calculated in the course of a Monte Carlo simulation.For a select set of t , the action is calculated for 50,000 configurations, where the datastarts to be recorded after 10,000 steps. This allows the system to thermalize. We performthis calculation for a wide range of t and fit (cid:104) S ( t ) (cid:105) /t to a high order polynomial. Thepolynomial is then integrated to obtain ln[ Z ( t ) /Z ( t )]. For the reference value t , we choose t = 1.An example of the entropy analysis can be seen in Fig. 17 where entropy is calculated fora temperature range of [0 ,
60] for fugacity y = 0 .
1. There appears to be a smooth logarithmicincrease in entropy with respect to t . This makes sense, as a BKT transition is a transitionof infinite order: a first derivative of the free energy is not expected to show a discontinuity,21 C ( τ ) τ t = 14t = 16t = 18t = 20t = 22 FIG. 11. The autocorrelation in the thickness observable, showing very long time scales, comparedto Figs. 9 and 10. It can be seen that by O (3000) time steps all memory of the initial configurationhas vanished. Figure previously appeared in [4]. and in fact should be smooth. To confirm this expectation, in Fig. 18 another exampleis shown for a higher fugacity, y = 1. Due to the higher potential barrier, tunneling wasa lot more erratic across temperatures, so the action required a higher order polynomialto reliably model the free energy. And as can be seen in the figure, the resulting curve issmooth with a monotonically increasing entropy.It is important to consider the effects of autocorrelation in this calculation of entropyas well. The observable is highly local, since it is a sum over the action density. Theresult is that the autocorrelation is significantly smaller than in the case of the thicknessmeasurements, which involve correlations over large scales. In order to quantify this smallerautocorrelation, in Fig. 19 we show the action for y = 1 generated for two separate randomseeds, which alters the initialization of φ ( x ) at the beginning of the simulation. The differencebetween the action on the two sets of configurations is well within the error bars, which takeinto account integrated autocorrelation time in the usual way.22 t h i ck ne ss t FIG. 12. Thickness versus t for y = 0 . L = 32 lattice. It can be seen that it is a linearfunction of t for this relatively weak value of y .As will be seen in subsequent figures, a more interesting behavior occurs for y > σ 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 FIG. 13. Thickness versus t for various y for an L = 64 lattice. It can be seen that a cross-overbetween different behaviors occurs for y ∼
1, with a nonlinear dependence of σ on t appearing forlarge y . σ t y = 4.0y = 8.0y = 16.0y = 32.0y = 64.0y = 128.0 FIG. 14. Thickness versus t , at rather small t , for various y for an L = 64 lattice. Some interestingfeatures emerge at the smallest values of t seen in this resolution. It can be seen that the steep risethat was observed in the previous plot simply becomes steeper, with its center moving to smaller t , as y increases. VII. CLUSTERING
We have analyzed the clustering of the field φ ( x ) into values around 2 nπ , n ∈ Z . Themethod for identifying a cluster is described in Appendix D. In Fig. 20 the various clustersfound by this algorithm are shown from a perspective where they are all visible. Fig. 20shows the same collection of clusters edge-on, where it is apparent that they are groupedinto layers. In order to produce results with clearly defined domains, we found that itwas helpful to use a random start with φ ( x ) drawn uniformly from [ − ,
20] and a verylarge fugacity of y = 15. We also set t = 0 . W and population P (the number of sites contained in the cluster) is show in Fig. 22, where50 lattice configurations were combined. It shows an approximate behavior of W ∼ (ln P ) (7.1)Note that the largest population and width is limited by the fact that this is an L = 64lattice. We are not aware of a theoretical argument explaining the behavior (7.1), but leave24 t h i ck ne ss t L = 8L = 12L=16 FIG. 15. The thickness observable rises essentially linearly with t for y = 0 .
1. However, the slopeincreases with L , as can be seen from these low L values. Note that this increase in slope onlyoccurs in a noticeable way for t greater than a critical value of around 10. A more precise estimateof t c ≈
18 will be obtained in our finite size scaling study below, summarized in Fig. 16. σ L t = 14t = 16t = 18t = 20t = 22
FIG. 16. The thickness observable transitioning between the two behaviors (5.5) at t c ≈
18. Here,we have used y = 0 . IG. 17. Entropy vs. temperature for fugacity y = 0 . y = 1. IG. 19. Average action vs. temperature for fugacity y = 1. it as a question for future investigations. VIII. CONCLUSIONS
One interesting avenue to explore is the possible existence of three phases in the XYtype systems, as shown in Fig. 24. Arguments for the existence of this more complicatedset of transitions were advanced in classic papers by Halperin, Nelson and Young [21–24].Although this feature has been argued for theoretically, and observed experimentally, it hasstill not been detected in numerical simulations to date. Thus, future work would be touse the more powerful and sophisticated simulation methods that are begun here to observethese types of behaviors. An important question is: “How would they be realized in the sine-Gordon model, if at all.” Does a theory of vortices contain enough information to capturethe full dynamics of the XY model in this regard? We leave this question open to futureresearch.We mentioned in the main body of the paper that in superconducting thin films one canunderstand the behavior as a superposition of BKT fluctuations (spin waves and vortices)27 [0]0 20 40 60 80 100 120 x[1]0 20 40 60 80 100 120 p h i [ x ] FIG. 20. The clustering of the field φ ( x ) into domains for a particular (thermalized) configurationwith t = 0 . y = 15 and L = 64. It can be seen that the bulk of the points cluster around φ = 0,but that other smaller clusters do form around the vacua at φ = ± π . x[0]020406080100120x[1]0 20 40 60 80 100 120 p h i [ x ] FIG. 21. The same configuration as in Fig. 20, but from a different perspective. and GL fluctuations. It would seem that the presence of GL fluctuations comes from thefact that the thin film is not truly two-dimensional, having a thickness d . This suggests anXY model with three dimensions, one of them being quite small. Alternatively, we couldformulate a SG model built on layers. One would like to establish the connection between28 Cluster Population (node count)505101520253035 C l u s t e r W i d t h ( s i z e a c r o ss c l u s t e r ) FIG. 22. Correlation between cluster population and cluster “width.” The width of a cluster isthe maximum size across in the two-dimensional domain. The population is the number of sites,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, andempirically is logarithmic, according to Eq. (7.1). these ideas. Currently we are developing numerical simulations that will analyze such layeredsystems, especially in the limit of weak interlayer interactions.
ACKNOWLEDGEMENTS
JG was supported in part by the Department of Energy, Office of Science, Office ofHigh Energy Physics, Grant No. de-sc0013496. Significant parts of this research were doneusing resources provided by the Open Science Grid [25, 26], which is supported by theNational Science Foundation award 1148698, and the U.S. Department of Energy’s Officeof Science. We are also appreciative of XSEDE [27] resources (Stampede), where othersignificant computations for this study were performed.29 F r e q u e n c y
20 15 10 5 0 5 10 15 20Phi(x)020406080100120140 F r e q u e n c y (a) (b)
20 15 10 5 0 5 10 15 20Phi(x)0100200300400500600 F r e q u e n c y
20 15 10 5 0 5 10 15 20Phi(x)0100200300400500 F r e q u e n c y (c) (d)FIG. 23. Distribution of field values on either side of the transition. Here, L = 64 and y = 0 . t = 14 and panel (b) is t = 22. These are the distributions for a single configura-tion, after 10000 thermalization updates. The two simulations were started from the same initialconditions with the same seed for the pseudorandom number generator, so that the differencesare entirely due to the value of t . It can be seen that t = 22 has longer tails, due to strongerfluctuations at the higher temperature. However, it is interesting that for this L , nothing dramaticoccurs as we cross the transition. In (c) and (d) we show the same analysis for L=128. ranslational orderorientational order translational disorderorientational order translational disorderorientational disorder FIG. 24. The three phases of algebraic order versus disorder, as advanced in the seminal papers ofHalperin, Nelson and Young (1978-79).
Appendix A: Free theory propagator1. Continuum
We generalize the free action to include a mass term, S ( θ ) = 1 t (cid:90) d x (cid:26)
12 [ ∂ µ φ ( x )] + 12 m φ ( x ) (cid:27) (A1)We must include a mass term to regulate the infrared divergence. Then∆( x, m ) = (cid:90) d p (2 π ) e ip · x p + m (A2)This can be evaluated with the help of results from Gradshteyn & Ryzhik (G&R) [28].Passing to polar coordinates and using G&R 3.339 (p. 357)14 π (cid:90) ∞ pdpp + m (cid:90) π dθ e ipx cos θ = 12 π (cid:90) ∞ pdpp + m J ( px ) (A3)Then using G&R 6.532 π K ( mx ) (A4)We are interested in the limit of small m , which is obtained from G&R 8.447 K ( z ) = − ln z ψ (1) + O ( z ) , ψ (1) = − γ (A5)31hus we finally obtain∆( x, m ) = − π ln mx − π γ + O ( m ) = − π (cid:20) ln m x γ (cid:21) + O ( m ) (A6)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 pathintegral (cid:42) n (cid:89) i =1 e i(cid:15) i φ ( x i ) (cid:43) = Z − (cid:90) [ dφ ( x )] e − S [ φ ( x )] n (cid:89) i =1 e i(cid:15) i φ ( x i ) , Z = (cid:90) [ dφ ( x )] e − S [ φ ( x )] (A7)by rewriting it in terms of a “source” J ( x ) = i n (cid:88) i =1 (cid:15) i δ ( x − x i ) (A8)so that (cid:42) n (cid:89) i =1 e i(cid:15) i φ ( x i ) (cid:43) = Z − (cid:90) [ dφ ( x )] exp (cid:18) − (cid:90) d x (cid:26) t ( ∂ µ φ ( x )) + m t φ ( x ) − J ( x ) φ ( x ) (cid:27)(cid:19) (A9)This can be computed, as usual, by completing the square, leading to (cid:42) n (cid:89) i =1 e i(cid:15) i φ ( x i ) (cid:43) = exp (cid:18) t (cid:90) d x d y J ( x )∆( x − y ; m ) J ( y ) (cid:19) (A10)Substituting (A8) we obtain (cid:42) n (cid:89) i =1 e i(cid:15) i φ ( x i ) (cid:43) = n (cid:89) i,j =1 (cid:18) m ( x i − x j ) (cid:19) t π (cid:15) i (cid:15) j (A11)To get closer to (3.6) we separate out factorslim δ → (cid:89) i (cid:18) mδ (cid:19) t π (cid:15) i · (cid:89) i The sine-Gordon model is supposed to capture the physics of vortices, where the fieldvariable φ ( x ) represents the angle of the U (1) (cid:39) O (2) phase of the order parameter: forinstance the condensate in a thin film superconductor, or the phase of the wavefunction ina superfluid. In a superconductor the vortices describe how magnetic field lines penetratethe bulk, and surrounding each vortex is an eddy current. If one considers two vortices suchas in Fig. 25, one can see that the eddy currents will tend to repel since the portions of thecurrents that are closest are going in opposite directions. Recall that unlike charges, for currents opposites repel. ppendix C: Map between different forms of the action An alternate form of the action that we have used early in our simulations is S (cid:48) [ φ ] = 1 t (cid:48) (cid:90) d x (cid:26) 12 [ ∂ µ φ (cid:48) ( x )] − g (cid:48) cos(2 φ (cid:48) ( x )) (cid:27) (C1)This can be related to the action in the body of the paper, Eq. (1.1), by defining φ ( x ) =2 φ (cid:48) ( x ) and relating ( t, g ) to ( t (cid:48) , g (cid:48) ) appropriately: t = 4 t (cid:48) , g = 4 g (cid:48) (C2)In [8] the following convention is used: S [ ϕ ] = (cid:90) d x (cid:26) β ( ∂ µ ϕ ) − z cos(2 πϕ ) (cid:27) (C3)The map is then ϕ = 12 π φ, β = t π , z = gt (C4)The thicknesses are related by σ H = 14 π σ J (C5)where on the l.h.s. there is the Hasenbusch et al. convention, and on the r.h.s. is ours. As aresult we end up with the prediction (5.5) from their result σ H (cid:39) βπ ln L + const. (C6) Appendix D: Domain analysis For every configuration, the values of φ tend to group together into distinct domainsduring the thermalization of the system. Plotting the system in terms of x, y, and φ yieldsa L × L scatter plot of particles as seen in Fig. 21. For this particular configuration thepotential barrier is low enough, and the temperature is high enough that sites are able totunnel to several domains above and below φ = 0, known as the threshold.The combination of fugacity and temperature for used to produce the configuration plot-ted in Fig. 21 makes it easy for an algorithm to recognize the clusters. In particular, it allowsK-means, a prototype-based cluster algorithm, to efficiently identify the clusters and theircentroids. One major issue that is present in using K-means, however, is that this algorithm34equires 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 methodis 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 usingthe k − dist values to define cluster radius and minimum point threshold. Once k number ofcluster radii (called Eps ) are identified, the algorithm uses K-means to classify those clus-ters. To aid in speeding up the process, centroids c ...c k are initialized within the local spaceof each Eps j , automatically assigning identified core points and border points to respectivecluster 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 toeven larger lattice sizes coupled with small deviations of φ . Many configurations of thisform can be too complex for this algorithm. A work-around for this is to analyze extremelylarge amounts of sites is to take the configuration data and project it onto a 2d planeusing principle component analysis (PCA). In PCA, a system with a certain amount ofunlabeled data is selected to form a s -dimensional sample space. The covariance matrix isthen computed as Σ = σ , σ , · · · σ ,j σ , σ , · · · σ ,j ... ... . . . ... σ i, σ i, · · · σ i,j (D1)Where σ i,i is the variance and σ i,j is the covariance. The eigenvalues and eigenvectorsare calculated from Σ, giving a set of potential axes (eigenvectors) that the data could beprojected upon. The j eigenvectors with the greatest eigenvalues are chosen and combinedwith the the s -dimensional sample space to form a s × j -dimensional matrix D , the reducedsubspace. The original samples X are then transformed onto the new subspace YY = D T × X (D2)In terms of our data, we considered the sample space of the configurations as a fibre bundleinstead 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 theanalysis in terms of a fibre bundle allows us to observe our PCA sample space with three35 FIG. 26. PCA of (thermalized) configuration 1000 for L = 256, y = 4, t = 1. features (x, y, and φ ) instead of just φ . This gives depth to our analysis, allowing for thecluster algorithm to identify the clusters quickly for a large number of increasingly largelattice sizes without any loss in accuracy. Fig. 26 shows the PCA of a configuration with alattice size of 256. Computation time is cut down by a factor of two and the cluster countand cluster width is not compromised. Appendix E: Hard core repulsion It is straightforward enough to implement the hard-core repulsion in the Coulomb gasmodel. For instance the fugacity expansion of the partition function as presented in themain text is Z CG = (cid:88) n ∈ Z y n n ! (cid:90) d x · · · d x n (cid:88) { (cid:15) } (cid:89) i 40 years of Beresinskii-Kosterlitz-Thoulesstheory , edited by J. V. Jos´e (World Scientific, Singapore, 2013) pp. 161–199, 1201.2307.[4] J. Giedt and J. Flamino, Proceedings, 35th International Symposium on Lattice Field The-ory (Lattice 2017): Granada, Spain, June 18-24, 2017 , EPJ Web Conf. , 14003 (2018),arXiv:1710.03188 [hep-lat].[5] S. Chui and J. Weeks, Phys. Rev. B14 , 4978 (1976).[6] S. Chui and J. Weeks, Phys. Rev. Lett. , 733 (1978). 7] S. R. Coleman, Phys. Rev. D11 , 2088 (1975), [,128(1974)].[8] M. Hasenbusch, M. Marcu, and K. Pinn, Physica A211 , 255 (1994), arXiv:hep-lat/9408005[hep-lat].[9] S. Duane, A. Kennedy, B. Pendelton, and D. Roweth, Phys. Lett. , 216 (1985).[10] A. Ferreira and R. Toral, Phys. Rev. E47 , 3848 (1993).[11] D. Espriu, V. Koulovassilopoulos, and A. Travesset, Phys. Rev. D56 , 6885 (1997), arXiv:hep-lat/9705027 [hep-lat].[12] S. Catterall and S. Karamov, Phys. Lett. B528 , 301 (2002), arXiv:hep-lat/0112025 [hep-lat].[13] J. M. Kosterlitz, J. Phys. C7 , 1046 (1974).[14] D. J. Amit, Y. Y. Goldschmidt, and G. Grinstein, J. Phys. A13 , 585 (1980).[15] J. Zinn-Justin, Quantum Field Theory and Critical Phenomena (4th Edition) (Oxford Uni-versity Press, New York, 2002).[16] S. Samuel, Physical Review D18 , 1916 (1978).[17] C. Mudry, Lecture notes on field theory in condensed matter physics (World Scientific, NewJersey, 2014).[18] B. I. Halperin and D. R. Nelson, Journal of Low Temperature Physics , 599 (1979).[19] G. ’t Hooft, Under the spell of the gauge principle (World Scientific, Singapore, 1994).[20] E. Fradkin, Field theory of condensed matter physics, 2nd edition (Cambridge UniversityPress, Cambridge, UK, 2013).[21] B. Halperin and D. Nelson, Phys. Rev. Lett. , 121 (1978).[22] B. Halperin and D. Nelson, Phys. Rev. B19 , 2457 (1979).[23] A. Young, Phys. Rev. B19 , 1855 (1979).[24] D. Nelson, Phys. Rev. B18 , 2318 (1978).[25] R. Pordes and et al., J. Phys. Conf. Ser. , 012057 (2007).[26] I. Sfiligoi, D. C. Bradley, B. Holzman, P. Mhashilkar, S. Padhi, and F. Wurthwein, 2009 WRIWorld Congress on Computer Science and Information Engineering , 428 (2009).[27] J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither, A. Grimshaw, V. Hazlewood, S. Lath-rop, D. Lifka, G. D. Peterson, R. Roskies, J. R. Scott, and N. Wilkins-Diehr, Computing inScience & Engineering , 62 (2014).[28] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series & Products, Fifth Edition , editedby A. Jeffrey (Academic Press, San Diego, California, 1994)., editedby A. Jeffrey (Academic Press, San Diego, California, 1994).