Graphical representations and worm algorithms for the O(N) spin model

Longxiang Liu ,Lei Zhang ,Xiaojun Tan and Youjin Deng2,

1 Simulation of Physical Systems Division,Beijing Computational Science Research Center,Beijing 100193,China

2 National Laboratory for Physical Sciences at Microscale,University of Science and Technology of China,Hefei,Anhui 230026,China

3 Department of Modern Physics,University of Science and Technology of China,Hefei,Anhui 230026,China

Abstract We present a family of graphical representations for the O(N) spin model,where N ≥1 represents the spin dimension,and N=1,2,3 corresponds to the Ising,XY and Heisenberg models,respectively.With an integer parameter 0 ≤? ≤N/2,each configuration is the coupling of ? copies of subgraphs consisting of directed flows and N -2? copies of subgraphs constructed by undirected loops,which we call the XY and Ising subgraphs,respectively.On each lattice site,the XY subgraphs satisfy the Kirchhoff flow-conservation law and the Ising subgraphs obey the Eulerian bond condition.Then,we formulate worm-type algorithms and simulate the O(N)model on the simple-cubic lattice for N from 2 to 6 at all possible ?.It is observed that the worm algorithm has much higher efficiency than the Metropolis method,and,for a given N,the efficiency is an increasing function of ?.Besides Monte Carlo simulations,we expect that these graphical representations would provide a convenient basis for the study of the O(N) spin model by other state-of-the-art methods like the tensor network renormalization.

Keywords: Markov-chain Monte Carlo algorithms,continuous spin models,graphical representations

1.Introduction

The O(N) spin model [1],also referred to as the O(N) nonlinear σ model in field-theoretic parlance,is central to the study of critical phenomena.It has been subjected to extensive studies over several decades via both field-theoretical approaches and statistical mechanical methods.Within the comprehensive family of O(N) spin models,the models for N=1,2,and 3 are of particularly significant interest as they correspond directly to the Ising,XY,and Heisenberg models,respectively.

To surpass the constraints of perturbation theory,the powerful numerical tool the Monte Carlo (MC) method [2] is widely used to simulate the O(N) spin model.However,the Metropolis scheme,which is frequently employed [3],undergoes significant ‘critical slowing down’,consequently impairing computational efficiency near the critical point [4].Addressing this issue,various non-local update algorithms have been put forward,encompassing multigrid techniques and cluster update algorithms [5].Among these,the cluster update algorithms prove to be the most effective for the O(N)spin model [6].On the flip side,for local update algorithms,the worm technique emerges as a promising contender.This algorithm,extensively deployed in both classical and quantum systems [7,8],exhibits superior efficiency at criticality,typically applied to graphical representations complying with Kirchhoff-like conservation laws.Fundamentally,a worm algorithm modifies configuration graphs by introducing two defects that infringe these conservation laws,moves them around according to the balance condition for MC simulations,and ultimately removes them once they reconvene.Hence,to develop worm-type algorithms for particular models,like the O(N) spin model,it is often required to first formulate their graphical representations.

Figure 1.Typical graphical configurations for the Ising model in representation (6) (a) and commonly used representation (7) (b).In both representations,a graph is constructed using undirected closed loops,with the number of lines on each bond representing the bond variable nij.Each of these lines contributes a weight factor WIs(nij) or Notably,while representation (b) allows a maximum of one line per bond,representation (a) can accommodate multiple lines on a single bond.

Graphical representations for the Ising (N=1) and XY(N=2) models are well established.For the Ising model,a specific graphical representation introduces a non-negative integer variable on each lattice bond,visually represented by the number of lines drawn on the corresponding bond(figure 1(a)).A configuration makes a non-zero contribution only when it satisfies the Eulerian condition,which requires each site to be connected by an even number of lines.On the other hand,a prevalent graphical representation for the XY model assigns two types of non-negative variables to each bond,symbolizing currents flowing in positive and negative directions,respectively.These currents,in their entirety,satisfy conservation at each site,reflecting the Kirchhoff flowconservation law (figure 2(a)).For the sake of convenience,configurations in these two representations are referred to as the Ising subgraphs and the XY subgraphs in this paper,respectively.

When considering an arbitrary value of N ＞2,studies on graphical representations have also made significant advances in the past decade.In 2010,Wolff put forward a graphical representation defined by a single set of bond variables,depicted by lines on the corresponding bonds [9].Each site incorporates a type of switch-board,creating pairwise connections among all neighboring lines.This approach implies that a range of configurations can be attributed to one specific set of bond variables,necessitating the inclusion of symmetry factors.A few years later,two distinct representations each depicted by N sets of bond variables were proposed.These representations view configurations as the combination of multiple subgraph copies [10,11].Specifically,a configuration in [10] is the union of one XY graph copy and N -2 Ising graph copies.Conversely,a configuration in [11] is the combination of N Ising graph copies.

In this work,we devise a systematic series of graphical representations for the O(N) spin model using a streamlined approach.This series is parameterized by an integer ?(0 ≤? ≤N/2).To establish these representations,we dissect the total spin vector into a weighted sum of ? copies of XY vectors and N -2? copies of Ising vectors.By incorporating currents (flows) for the XY vectors and undirected bonds for the Ising vectors,we obtain the graphical representation after integrating out the spin degree of freedom.The resulting configuration consists of ? copies of XY subgraphs and N -2? copies of Ising subgraphs.Notably,when ?=0 and 1,our representations reduce to those proposed in [10,11].Further,by introducing two defects in these graphical representations,we develop corresponding worm algorithms for any arbitrary N and potential ?.Using these algorithms,we conduct large-scale simulations on a simple-cubic lattice with a linear size of up to L=96 for N=2,3,4,5,and 6,across all possible ? values.We deduce and compare the dynamic exponents for each variant of the worm algorithms from the fitting of the integrated autocorrelation time for an ‘energylike’ observable.Our simulation results reveal that the algorithm efficiency increases with the number of XY subgraph copies in a single configuration,i.e.a larger ? promotes higher efficiency.In the best cases,dynamic exponents can reach z=0.20(1),0.32(2),0.22(1),0.26(3),and 0.16(2)for N=2,3,4,5,and 6,respectively.Consequently,the efficiency of worm-type algorithms is comparable with that of cluster update methods [12],significantly outperforming the Metropolis algorithm which has a dynamic exponent z ≈2.

In addition to MC simulations,researchers are exploring other advanced numerical methods for applications on the O(N) spin model,among which the tensor network renormalization (TNR) [13] is one of the most promising.It starts by representing the partition function or the ground state wavefunction as a tensor network state that is formed by taking the product of local tensors defined on the lattice sites.However,applying this method to the O(N) spin model is challenging because each spin has infinite degrees of freedom,and the conventional method of constructing the local tensors fails.For N=1 and 2,researchers proposed a novel scheme to construct the tensor representation using character expansions [14].These representations can actually be derived by slightly reformulating the conventional bondbased graphical representations for the Ising and XY models.Therefore,we expect that the development of graphical representations would provide a convenient foundation for investigating the O(N) spin model with TNR and other numerical methods.

This paper is structured as follows: in section 2,we provide a succinct introduction to the Hamiltonian and the partition function for the O(N) spin model,and then revisit the graphical representations for the Ising and XY models before deriving the family of graphical representations for the O(N)spin model.The worm algorithms are introduced in section 3,followed by the presentation of our simulation results in section 4.We conclude the paper with a discussion in section 5.

2.Hamiltonian and graphical representations

The O(N) spin model characterizes N-component spin vectors with inner-product interactions on a lattice,as represented by the Hamiltonian:

In the following parts,we begin by providing a brief overview of the well-established graphical representations for the Ising (N=1) and XY (N=2) models.Using these as a basis,we then construct a series of graphical representations adaptable to any given value of N.Note that,for the sake of simplicity,our illustrations in this section are confined to the square lattice.However,our derivation is applicable to any type of lattice in any spatial dimension.

2.1.N=1: the Ising model

The Ising model characterizes unit spin vectors on the lattice that can either point upwards or downwards.Consequently,its partition function is given by

where E(i) signifies the set of nearest neighboring sites of site i.As depicted in figure 1(a),if we draw a number of undirected lines on each bond equivalent to the corresponding bond variable nij,we get a graph for every summation term in the partition function equation (4),each with a different number of lines on each bond.The latter factor in each summation term implies that only when nijfulfills the constraint at any site i:

will the entire summation term make a non-zero contribution.This requirement implies that the degree of each site in the graph should be even,or alternatively,the graph can be seen as composed of closed loops.These non-zero contribution graphs are termed valid configurations,and if we use GIsto denote all of them,the partition function can be written in a graphical form:

where WIs(n) ≡Kn/n! stands for the weight factor contributed by each bond.

Consequently,we construct a graphical representation for the Ising model,in which configurations are composed of closed-loop graphs founded on the bond framework.In this configuration graph,each bond imparts a weight factor WIs(n),dictated by the number of lines,n,present on it.

The above graphical representation differs from the commonly employed one for the Ising model,as depicted in figure 1(b).In the conventional representation,the constraint of even degree at each site still applies,but the bond variables nijcan only adopt two values: 0 or 1.In graphical terms,loops in the graph should not share any bonds or traverse a bond more than once.This results in the following expansion for the partition function:

2.2.N=2: the XY model

The XY model characterizes unit planar vectors at each lattice site,with the partition function expressed as:

After integrating over all the angle variables φi,we can represent the partition function as:

This implies that the combination of positive-and negativedirected lines can be effectively treated as a conserved current(flow) defined on the bond structure.Consequently,we can represent valid configurations as graphs composed of directed closed loops.Furthermore,in these graphs,each bond can have multiple lines along both directions simultaneously.Let Gxydenote all such graphs,the partition function Zxycan then be expressed as:

where the bond variable mijis an integer taking values 0,±1,±2,…,the bond weight factor(m)is defnied as Im(K),corresponding to the modified Bessel function of the first kind with order m,anddenotes all graphs that also consist of directed closed loops.However,it imposes a constraint that a bond can only be traversed from a single direction at a time,implying that all lines on the same bond should follow the same direction: either positive or negative.In this representation,the quantity and direction of lines on a bond represent the absolute value and the sign of mij,respectively.A typical graphical configuration is illustrated in figure 2(b).Following a similar approach to the Ising model,we can apply a reformulation of the graphical representation in a tensor form as follows:

2.3.Arbitrary N

To obtain graphical representations for any integer N ≥1,we first analyze the structure of the N-dimensional spin vector S and discover that it can be decomposed into a sum of ? copies of two-dimensional (2D) and N -2? copies of one-dimensional (1D) orthogonal spin vectors,with 0 ≤? ≤N/2.This can be expressed as:

or more concisely written as ∣A∣=1,if we define an (N -?)-dimensional vector of all these magnitude coefficients A ≡(A1,…,A?,B1,…,BN-2?).

Given that all these decomposed low-dimensional vectors S(α)and s(β)are orthogonal to each other,the interaction between two spin vectors can be reduced to this decomposition:

Therefore,we can perceive a spin vector as a weighted sum of several 1D or 2D unit spin vectors.Moreover,the interaction between two spin vectors can be seen as a superposition of interactions between the low-dimensional spin vectors that constitute each original spin vector and their corresponding components in the other.The breaking down of spin vectors and interactions into these lower-dimensional components reminds us of the Ising and XY models,which also utilize 1D or 2D unit spin vectors and comparable interactions.It is logical to hypothesize that we can derive a representation in which each configuration is built from several instances of configurations specific to the Ising or XY models.

Here,G is the direct sum of all the graph setsandThis implies that a graph in G can be perceived as a superposition of ? copies of graphs for the XY model and N -2?ones for the Ising model.Mathematically,for any α or β,the bond variables must satisfy the following constraints for any site i:

The weight factor for each bond,denoted as W({m±(α)},{n(β)}),is the product of all weight factors from each decomposed graph.This can be expressed explicitly as:

Figure 3.Typical graphical configurations for the O(3) spin model at different ? values: (a) For ?=0,the configuration is composed of three copies of the Ising model graphs,as shown by orange,olive,and purple undirected lines,respectively;(b) for ?=1,the configuration includes one copy of the XY model graph,denoted by the cyan directed lines,and an additional copy of the Ising model graph,represented by orange undirected lines.In both instances,each bond or site contributes a weight factor as specified in equations (23) or (24),respectively.

It is important to emphasize that if we expandandinto their more conventional forms,as seen in equations (14) and (7),we will not be able to independently separate the magnitude coefficients from the bond weight factors.Consequently,it would not be possible to analytically integrate them out beforehand to calculate the site weight factors.

Thus,for each potential value of ?,we obtain a graphical representation.In this representation,every graph configuration comprises ? instances of XY graphs and N -2? instances of Ising graphs.Each bond and site within these configurations contribute a weight factor,as indicated in equations (23)and (24) respectively.For instance,when N=3,there are two viable values for ?: 0 and 1.As depicted in figure 3,when?=0,the graphical configuration can be interpreted as being comprised of three Ising graphs.Conversely,when ?=1,the configuration can be viewed as consisting of a single XY graph and one Ising graph.

Furthermore,the partition function in the graphical description can be represented in tensor form as follows:

where the tensor T(i)is defnied as

Consequently,each graphical representation yields a corresponding tensor representation,marking the initial stride towards the application of the TNR method to the O(N) spin model.

3.Worm algorithm

Based on the series of graphical representations derived in the previous section,we formulate worm algorithms to simulate the O(N) spin model for arbitrary N and ?.In these algorithms,starting from an arbitrary valid configuration,we introduce two defects named Ira and Masha into a randomly selected copy of an XY or Ising subgraph with a certain acceptance probability,causing the constraint in equation (22) to fail for the corresponding α or β of the selected subgraph at these two sites.In the defected subgraph,both Ira and Masha are connected by odd lines if it is an Ising subgraph,and one extra current flows in (out) for Ira (Masha) if it is an XY subgraph.To update the selected subgraph,we then move the defect Ira randomly by deleting (adding) a line or deleting(adding) a positive (negative) current,depending on whether it is an XY or Ising subgraph,until Ira meets Masha again.

We designate the space built by configurations with two defects as the ‘worm space’.Configurations within this space are weighted according to the two-point correlation function:

Due to the unitary modulus of the vector S,the correlation function Z(I,M)simplifeis to the partition function Z when I=M.Thus,we can equivalently view configurations with both Ira and Masha on the same site as physical configurations.In this context,we can measure the observables on these configurations,disregarding the presence of both defects.

The procedures we use to update the Ising or XY subgraphs are identical to those employed for the Ising or XY model.A comprehensive description of the entire algorithm for a specific ? graphical representation can be found in algorithm 1.

Algorithm 1.Worm algorithm for the O(N) spin model

4.Numerical simulations

4.1.Setup

To probe the dynamical critical behavior of our worm algorithms,we conduct extensive simulations at the critical points for each N at various ? values.First,we introduce the observable we measure:

This observable counts all the lines in the graph,irrespective of their types.Intriguingly,the thermodynamic average of this quantity is effectively the negative of the system’s energy,scaled by a temperature factor T,that is,E=-T〈N 〉.This relationship can be readily confirmed by evaluating the formulaE=-? lnZ?K.

Based on the time series of the quantity N in our Monte Carlo simulations,we compute the integrated autocorrelation time,τint,N,which is defined as:

Table 1.Critical couplings Kc we applied to simulate at for each N,and corresponding dynamical critical exponentszint,N at different ?which are roughly estimated via fitting the data by equation (32) at large system sizes (L ＞20).

Figure 4.Autocorrelation function ρN (t)at criticality versus tτint,N for N=2.At ?=1,ρN (t)is very close to a pure exponential decay,while at ?=0,it exhibits a crossover from a pure exponential to a complicated behavior as L increases.

In this study,we conduct simulations of systems at criticality for all possible values of ? across different N ranging from 2 to 6 on 3D cubic lattices with various system sizes L,extending up to 96.Periodic boundary conditions are implemented in all simulations.For every pair of N and L,we perform a total of 2 × 108measurements,with L/2 worm steps between two consecutive measurements.The critical couplings we use for each N are detailed in table 1.

4.2.Results

We frist examine the autocorrelation functionρN(t)at the critical points Kcacross various system sizes L for different N and ? values.As an illustrative example,figure 4 shows the variation ofρN(t)as a function of time t for various system sizes L when N=2.The figure reveals that for ?=1,the decay ofρN(t)nearly follows a pure exponential across all system sizes.However,for ?=0,while the decay is approximately exponential for smaller L,the trend deviates from this behavior as L increases,leading to an intersection of lines representingρN(t)at different L aroundtτint,N=3.For N=3,4,5,and 6,at all possible ? values,the behavior ofρN(t)is akin to that observed for N=2 and ?=0.The underlying mechanism driving this phenomenon,however,remains unclear at present.A plausible scenario is that the decay ofρN(t)is effectively the mixing of two exponential functions asρN(t)≈When exponents τ1＞τ2and amplitudes a1＜a2,the decaying behavior ofρN(t)is dominated by τ2for small t but is eventually governed by τ1for sufficiently large t.

Utilizing the formula in equation (30) and the truncation strategy,we derive the integrated autocorrelation timeτint,Nat criticality.Based on the patterns observed inρN(t),we select the windowing parameter c to be 8 for N=2,?=1,and 20 for the remaining cases.Figure 5 presents the relationship betweenτint,Nand system size for N=2,while figure 6 depicts the same for N=3,4,5 and 6.Across all N and ?values,τint,Nappears to increase approximately as a power law with growing L,thereby indicating the manifestation of the critical slowing down phenomenon.

Figure 5.Integrated autocorrelation timeτint,N at criticality versus system size L for N=2.Here,τint,N is in the unit of one system sweep which equals to the number of bonds in the system.Blue and red points are for ?=0 and 1,respectively.Dashed lines show the corresponding fitting results at large L with the fitting formula equation (32).The ranges of them also mark the points we applied in the fittings.

Figure 6.Integrated autocorrelation timeτint,N at criticality for N=3,4,5 and 6 at different ?,in the unit of system sweep.Points of different colors represent data for different ?,and dashed lines here have the same meaning as those in figure 5.

To quantify this,we fit the data forτint,Nat larger L values for all N and ? using the following formula:

The results of these fits are represented by the dashed lines in the corresponding figures,while the estimates forzint,Nat different N and ? are detailed in table 1.The findings suggest that for each N,as ? increases,zint,Ndecreases monotonically,implying a higher efficiency of the worm algorithm at larger ?.This observation can also be directly gleaned from the slopes of each line in figures 5 and 6.Even for the worst case (?=0),the dynamic exponent z is significantly smaller than that for the Metropolis algorithm,which has z ≈2.0,meaning that the worm algorithms are much more efficient than the Metropolis method.In comparison with the Swendsen–Wang (SW)cluster method,which has z=0.46(3) for the three-dimensional Ising model [12],the worm algorithm for the worst case (z ?0.5) seems to be comparable,and for the best case(z ?0.3) it seems to be somewhat more efficient,particularly for large N.Taking into account the existence of the ‘critical speeding-up’ phenomenon for susceptibility-like quantities in worm-type simulations [21],we conclude that the worm algorithm is at least as efficient as the SW cluster method.

5.Discussion

We present a successful reformulation of the O(N) spin model into a series of graphical representations,introduce corresponding worm algorithms to facilitate simulations and compare their dynamical critical behaviors.It is observed that the algorithm corresponding to the representation with a larger ?demonstrates a smaller dynamical critical exponent,implying a superior level of efficiency.Additionally,we employ these graphical representations to derive corresponding tensor representations,which we believe will be useful for TNR or other numerical methods.

The availability of multiple graphical representations allows us to investigate the underlying mechanisms contributing to the superior efficiency of these algorithms.This exploration provides valuable insights into their effective operation.To further enhance practical efficiency,our next step focuses on reducing the constant factor A of the autocorrelation time.We propose achieving this improvement by incorporating irreversible Monte Carlo techniques,specifically the lifting technique [22–25].Our preliminary test verifies the potential of this approach,motivating its adoption in our research.

In addition to providing a series of new graphical representations and corresponding algorithms,our work presents a general idea for reformulating systems describing vectors with multiple components into graphical representations.This involves breaking the vectors into the sum of analogs with fewer components,which are easier to represent graphically.The configurations can then be decomposed into a coupling of corresponding subgraphs,opening up prospects for future graphical representations of the CP(N -1) spin model.

Moreover,this approach can be applied to overcome the complex action problem when solving the O(N) spin model with a chemical potential μ coupled with conserved charges originating from the global O(N) symmetry [10].By binding the components related to the chemical potential as a copy of the XY spin vector during the decomposition of the total spin vector,the chemical potential term can be reduced to an extra weight factor that depends on the conserved charges in the graphical representation.In this new representation,the positive and negative currents for the chemical-potentialrelated XY spin vector are no longer conserved unless the extra integer variables defined on each site that stand for the conserved charges are taken into account.It is also noted that these graphical representations can provide a convenient and effective platform to study many physical problems like universal conductivity or halon physics [26–28].Overall,our approach offers a versatile tool to reformulate complex systems in a more accessible and intuitive graphical form.

Acknowledgments

We thank Nikolay Prokof’ev and Boris Svistunov for their valuable discussions.This work was supported by the National Natural Science Foundation of China (under Grant No.12275263),the Innovation Program for Quantum Science and Technology (under Grant No.2021ZD0301900),and the Natural Science Foundation of Fujian Province of China:2023J02032.

Appendix A.Computation of equation (24)

where θ1,θ2,…,θN-?-1∈[0,π/2].It is standard that the surface element can be expressed as

By using the identity that

where Γ(·) is the Gamma function,equation (24) can be obtained directly.

Appendix B.Acceptance probabilities