Monte Carlo simulation with Wolff algorithm for scaling behavior of two dimensional XY model with KT phase transition

The property of XY model with its KT phase transition has been discussed in many aspects such as specific heat and magnetic susceptibility along with correlation function. Also, this is one of the original models for topological phase transition which gives a new type of phase transition without symmetry breaking. Since the calculation of these properties above can be achieved by Monte Carlo Method with Wolff algorithm, the visualization of vortex and antivortex is also an interesting task which is accomplished in this article. By successfully modeling the system and getting the step-by-step configuration, the relation between the generation of vortex and the variation of temperature can be concluded in a qualitative way. The calculation of spin correlation function shows the phase transition between disorder and quasi long-range order which represents the topological phase transition. As a reliable method, the Monte Carlo method with Wolff algorithm can be a choice of generating the training set for neural network which can calculated much more sophisticated case in a relatively short time and the transition between classical simulation and quantum simulation can be approached as mentioned in the end of this article.


Introduction
The KT phase transition which was first introduced by Berezinskii -Kosterlitz-Thoulessis, the one that cannot be described by order parameter [1][2][3].This kind of phase transition is also called topological phase transition which is explored in the melting of 2D crystal with dislocations as the very beginning [4].From then on, the research of topological matter in condensed matter physics has become a popular topic.Before the idea of KT phase transition, it is universally convinced that O(2) nonlinear sigma model performs no phase transition due to the protected continuous symmetry.Then the phase transition in 2D He-4 superfluid thin film demonstrates the existence of a kind of phase transition free of symmetry broken in O(2) system.In fact, He-4 thin film is perfectly elaborated by the theory of KT phase transition-its superfluid phase factor can be interpreted by XY model.For two-dimensional XY model, its Hamiltonian is given by:   = − ∑ <,>  ⃗  ⋅  ⃗  = − ∑ <,> (  −   ) Where < i, j > denotes nearest neighbor.In fact, it can be classified into nonlinear sigma model family whose O(N)-invariant partition function is where n is an O(N) isovector of unit length.This partition function gives the simplest nearestneighbor ferromagnetic interaction between the isovectors.If n has just one component, then the only choice of it is \pm 1 which is actually Ising model.In 2-dimensional case, the unit vector n can be described by θ(r ⃗) in which r ⃗ stands for lattice site.Then the partition function is in terms of θ(r ⃗) : This work aimed at illustrating the behavior of correlation function of spin in both low-temperature and high-temperature limit which separately corresponds to long-range power-law: (4) and short-range exponential limit: It can be seen that with the increasing of the temperature, the degree of disorder continues to become higher which results in the rapid loss of correlation between the spins on each lattice site.In addition, the renormalization group (RG) equation of Sine-Gordon model will be also discussed later, which gives a more exact prediction to the scaling law of XY model.This model is derived by Villian approximation, Hubbard-Stratonovich transformation from the initial XY model.By attaching Poisson summation formula, the Sine-Gordon model can be described as: The first term in the free energy of Sine-Gordon model is in the similar form from the free energy of XY model: Except the position of coupling constant K.This is a reflection of duality which can be retrospected from 2D Ising model's self-dual (the dual that gives the same partition function before and after).By using the technique of dual, people got the exact critical temperature of 2D Ising model even before Onsager.
As for the second term, it is added for the energy gap between excited states.In the continuum limit, this model can be written in the form of path integral: When y = 0, the excitation of ϕ is gapless which corresponds to low temperature case in XY model.On the contrary, when y ≫ 0, the minimum point of the second term in free energy is ϕ = n(n is integer).The expansion near one of the minimum points ϕ = n + δϕ can be illustrated as In other words, the low excitation δ got the energy gap which is proportional to √y.If the parameter y is relatively large, the system is always gaped which correspond to high temperature and its shortrange correlation function.
Guided by the theoretical approach above, people are interested in the scaling behavior of this system near the critical points.Also, the visualized simulation has been achieved by cluster algorithm which can also figure out the number of vortices and antivortices.In this vision, one can figure out the change in temperature which is motored by the generation of vortices and anti-ones.It is worth noticed that KT phase transition can't be described by Landau order parameter.Instead, the phase transition is described by topological invariant-the circle integration of the gradient of the angle of spin.
In fact, if one changes the form of angle as the function of position, this integration remains invariant.Because the value of I only depends on the number of vortices in the area of integral, I is called a topological invariant.So, the most direct way of observing the dynamics of KT phase transition in XY model is counting the number of vortices and antivortices which will be discussed in section 4.
The simulation of 2D XY model focused on its thermal properties and the behavior of correlation function which denotes its unique topological phase transition.By Monte Carlo method with Wolff algorithm, the critical phenomena of magnetization and energy has been distinguished diagrammatically.Also, the correlation functions simulated at different temperatures show its correspondence to the theoretical prediction in Sine-Gordon model solved by Wilson renormalization.

High temperature expansion and low temperature expansion
The first check of the spin correlation function can be achieved by taking low and high limit of temperature in order to derive the behavior of correlation function which reflects the degree of order [5].

High temperature expansion
For high temperature, β = 1/T is a relatively small parameter so that one can make the expansion of Boltzmann factor based on it, which is The spin correlation function takes the similar form: To calculate the partition function Z, brute force calculation is necessary.By inserting the term H XY , the form of partition function is discovered.
Then, one important property is about the integral by angle.
It figures out that if there are no enclosed path constructed by other lattices' spin than considered spin, the integral showed above will vanish.For closed loop for integral, the result is non-zero.Obviously, the m = 1 term in partition function cannot construct a closed loop.It shows that the non-zero term is at least two orders, i.e m = 2 as shown below: For the same reason as  1 ,  3 also vanishes for the non-existence of closed integral loop.Then by the same process, it can be derived that  4 ∝ K 4 .In conclusion, the partition function under high temperature condition can be written as  ≈ 1 + K 2 + K 4

+ ⋯
As for the calculation of numerator in correlation function, the extra two spin-angle needs to be considered compared with the calculation of partition function: The only condition for non-zero integral are the open paths which connects these two spins.Thus, the non-zero term of lowest order is the segment that connects two spins.
Combining the numerator and denominator, the final solution to the behavior of spin correlation function is Where the correlation length ξ is defined by ξ = 1/(ln(1/K) .This exponential decay in spin correlation function implies the disorder at high temperature.This phenomenon contributes a contract condition to the power-law correlation function which represents order state at low temperature.To find out the behavior of correlation function at low-temperature, the technique of low temperature expansion is introduced.

Low temperature expansion
In the limit of K → ∞, the only terms which should be taken into consideration is by the small variation of θ(r ⃗) which ensures the continuous limit where one can substitute the summation by integration.Then the partition function is given by: Where a represents lattice constant.Applying Fourier transformation () = ∫  2  (2) 2   ().One could get the Gaussian partition function which gives a solvable model by means of path integral.For a Gaussian action, the correlation function in momentum space is given by Where the small mass term r is attached for convergence.Inverse Fourier TR the correlation function in momentum space, the correlation function in real space is given by Then the correlation function of spin at low temperature is given by According to the proper behavior of correlation function at both high and low temperature, the next task is finding the critical temperature between the two phases.It is non-trivial for the arbitrary direction of the spins on the lattice sites which would motor the evolution of temperature.Recall that in 2D Ising model, people used renormalization group to find the critical temperature between ferromagnetic and para-magnetic phase by duality.For 2D XY model, the renormalization approach is also available for Sine-Gordon model which is the dual of 2D XY model.

RG equation of Sine-Gordon model
One of the main problems based on KT phase transition is the inefficiency of the theory of Landau mean field theory and the fix of Gaussian fluctuation.Also, the anomalous dimension which implies the deviation to free field theory by renormalization can be calculated.By Wilson perturbation, the RG equation of Sine-Gordon model is derived [6]: Where the high frequency part in renormalization process is defined by [Λ/b, Λ].
To analyze the situation near the critical point K = K * , the variables z = y ˆ2, x = 2 − πK << 1are defined, then the RG equation can be solved as Due to the initial value A, the solution to the system has three different forms.
At the critical point K = K * = 2/π, y = 0, the system is equivalent to free field theory.So one can substitute the K term in the low-temperature correlation function by K = K * = 2/π to get the spin correlation function at critical point Compared with the definition of a standard form correlation function with scaling exponents: In d = 2, the anomaly dimension η = 1/4 which is large enough to declare the difference between 2D XY model and free field theory with Gaussian fluctuation.For this reason, Landau's mean field theory is disable in this model and the author turned to the simulation by Monte Carlo method with Wolff algorithm.

Monte Carlo simulation of 2D XY model
Monte Carlo method is a proper method which can be applied in systems with complex Hamiltonian that is difficult to calculate.Wolff algorithm gives a thread of selecting reasonable hopping between different states-states of lower energy is always acceptable and higher energy should be decided due to Boltzmann factor.

Thermal properties
In equilibrium, the measurements can be derived by Boltzmann distribution through several average quantities through the simulation.The simulation includes magnetization, magnetic susceptibility, energy and specific heat into consideration.The initial state is set to be ferromagnetic with 12 × 12 size of periodical square lattice aims to study bulk property of this model.The condition of convergence is about both energy and magnetization which gives a state in equilibrium at certain temperature.From T = 5J to T = 0.13J , successful simulation of the critical behavior in this range of temperature is derived as shown in Figures 1-4.From the result of simulation, the critical phenomena near the theoretical T = π/2J is dramatic.The high magnetization given by the ferromagnetic initial state is broken near the phase transition point and gradually flows to zero with the increasing temperature.The similar situation occurs in the result of energy and thermal capacity.Rather than the systems with spins fixed at 1 or -1 , the arbitrary direction of spin can form more exotic diagram such as vortex.

Dynamical process of the evolution
To analyze the behavior of the system in spin space, the configuration of the spin system step by step in the process of Monte Carlo at certain temperature is selective drawn.To figure out the order and disorder, one low temperature T = 0.01J , one critical temperature T = π/2J and one high temperature T = 5J were chose to observe.For T = 0.01J,π/2J, 5J , the three configurations 11110 Monte Carlo step,112 Monte Carlo step and full Monte Carlo step) are shown in Figure 5, 6, 7, 8, 9, 10, 11, 12 and Figure 13.
Since the red and blue rectangle stands for vortex and antivortex respectively, one could figure out that the location and number of unit vortex is fixed after several steps of Monte Carlo.In comparison, the large vortex can be only observed in the condition of T = π/2J and T = 5J.At T = 0.01J which is a low temperature, the spins are aligned in the way of ferromagnetism.By visualization of the evolution of spins, the phase transition between quasi long-range order phase and disordered phase is clear.So called quasi long-range order cannot be regard as 'true' long-range order because Recall that the Mermin-Wagner theorem tells that in the symmetry breaking which generates massless Goldstone modes, the fluctuation leads to the destruction of long-range order at any finite temperature in case the dimension is no more than two [7].However, in this 2D model, this theorem seems to be wrong under the existence of quasi long-range order.Other than the normal phase transition, this kind of phase transition is called topological phase transition which is motivated by the condensation of topological defects.

Illustration of topological defect
As shown in 4.2, the symbol of topological phase transition is the generation of topological defects [8].Since the angle describing the orientation of spins is defined up to 2π, it is possible to gain a 2nπ rotation through a closed path.Here n is called topological charge which is robust under the continuous deformation of the form θ(r ⃗).In other words, it is impossible to find a continuous way of deformation to change the configuration to ground state which is ordered with zero topological charge.
The vortex which has unit topological charge discussed in XY model is one of the elementary topological defects.In a complete circle centered on the defect, the orientation of spin varies ±2π.For a large vortex in which the lattice structure can be ignored, the energy cost from a single vortex of charge n can be estimated at low temperature 1in case of small gradient for θ(r ⃗)) by Where the first term describes the energy contribution from the core area of vortex since the second and dominant term describes the contribution from the large outside which diverges logarithmically in the size of system L.These two areas is distinguished by a circle of radius a.Such a huge energy cost prevents the spontaneous formation of the vortex at low temperature.When taking the configurational entropy of possible location of vortex into consideration, the partition function of a configuration with single n charge vortex is approximated as: By partition function, one could easily check that when free energy is smaller than zero, i.e K < 2/(n 2 π) which is the situation that tends to generate the vortex.For XY model, n = ±1 so that the phase transition point of this system is T = 0.5πJ.This result corresponds to the result of Sine-Gordon model and is more concrete by its thermodynamic meaning of spontaneity.

Scaling behavior of spin correlation function
The definitive result obtained in the exploration of the quasi-long-range order and disorder is the behavior of spin correlation function shown in Figure 14.From the limitation of the lattice constant, the first point of cross-correlation function is located at 1 which is the nearest from original point on which the correlation function is auto-correlation function.For this reason, one could always observe a immediate decreasing at the beginning when the temperature is higher than phase transition point.Though the abnormal increasing at the distance larger than 8.49(6√2) is sometimes severe, consider the bar and confidence interval, most of them can be removed.Here T = 5J and T = 0.5πJ are chosen for exponential and power-law fitting in Figure 15 and Figure 16.

Figure 1 .
Figure 1.Magnetization-Temperature. Near the critical temperature, the magnetization losses rapidly which corresponds the transition to disorder.

Figure 2 .
Figure 2. Magnetic susceptibility-Temperature.Since the spins behave in a more stochastic way, the system is in paramagenetic state which responds to external field strongly.

Figure 3 .
Figure 3. Energy-Temperature.Figure 4. Specific heat-Temperature.From the result of simulation, the critical phenomena near the theoretical T = π/2J is dramatic.The high magnetization given by the ferromagnetic initial state is broken near the phase transition point and

Figure 4 .
Figure 3. Energy-Temperature.Figure 4. Specific heat-Temperature.From the result of simulation, the critical phenomena near the theoretical T = π/2J is dramatic.The high magnetization given by the ferromagnetic initial state is broken near the phase transition point and

Figure 5 .
Figure 5. T = 0.01 J with 1/10 Monte Carlo step.At this time, the whole system is actually in one single cluster.

Figure 8 .
Figure 8.  = /2 with 1/10 Monte Carlo step.The configuration consists merely one vortex.From the result of next two figures, by this time the system is approximately in equilibrium.

Figure 13 .
Figure 13.T = 5J with full Monte Carlo step.Multiple vortices can be observed in this configuration.In comparison, the large vortex can be only observed in the condition of T = π/2J and T = 5J.At T = 0.01J which is a low temperature, the spins are aligned in the way of ferromagnetism.By visualization of the evolution of spins, the phase transition between quasi long-range order phase and disordered phase is clear.So called quasi long-range order cannot be regard as 'true' long-range order

Figure 14 .
Figure 14.Spin correlation function-distance for  = 0.01, 0.5, 0.67, 1.2, 0.5, 5.Since the improper data with unacceptable anomaly is excluded by confidence interval algorithm along with the zero point in the fitting of power-law, the final result is satisfying-at T = 5J , the correlation length is 1.20 compared with the theoretical ln 5 ≈ 1.61; at T = 0.5πJ, the coefficient of power-law is -0.2388 compared with -0.25.Since the standard correlation function at critical point is ⟨ ⃗  ⋅  ⃗  ⟩ =  ∝ | ⃗  −  ⃗  | −(−2+)