Exploring the worldline formulation of the Potts model

We revisit the issue of worldline formulations for the q-state Potts model and discuss a worldline representation in arbitrary dimensions which also allows for magnetic terms. For vanishing magnetic field we implement a Hodge decomposition for resolving the constraints with dual variables, which in two dimensions implies self-duality as a simple corollary. We present exploratory 2-d Monte Carlo simulations in terms of the worldlines, based on worm algorithms. We study both, vanishing and non-zero magnetic field, and explore q between q = 2 and q = 30, i.e., Potts models with continuous, as well as strong first order transitions. 1 Introductory remarks The presentation of the Prokofev-Svistunov worm algorithm in [1] was a milestone for the developement of worldline Monte Carlo simulations in statistical mechanics and lattice field theory. Not only did worldline2 representations and suitable algorithms considerably increase the numerical efficiency of Monte Carlo simulations, but also helped to solve some complex action problems, in particular those emerging in lattice field theories at finite chemical potential (for some reviews see, e.g., [2–5]). It is important to note that there is no unique general approach to worldline techniques and suitable corresponding algorithms, but different types of systems have to be treated individually. Over the years new tools and techniques have been developed for worldline representations and now also examples of worldline representations for some non-Abelian symmetries are known, which pose a considerably more challenging task for a worldline approach. In this paper we explore worldline techniques for yet another system with a non-Abelian symmetry, the q-state Potts model [8] with the permutation group of q elements as its symmetry group (see [9] for a review). In the context of lattice field theories mainly the 3-state Potts model, which is equivalent to the Z3 spin model, has been analyzed [10] – [23]. Z3 is the center group of SU(3), the gauge group Member of NAWI Graz. We point out that in the context of gauge theories the objects that correspond to worldlines are actually discretized worldsheets described by occupation numbers placed on plaquettes. For notational convenience we continue to use ”worldlines” but remark that often very similar techniques (see, e.g., [6, 7]) lead from worldlines for spins systems or lattice field theories to the corresponding ”worldsheet” representations of gauge theories with the same symmetry group (which is of course a local symmetry in the gauge theory context). 1 ar X iv :1 91 1. 12 72 8v 1 [ he pla t] 2 8 N ov 2 01 9 of QCD, such that the 3-state Potts model serves as a simple model theory for studying aspects of the strong interaction. We here use a worldline representation for the Potts model at arbitrary q that also allows for the inclusion of magnetic terms, and work out the representation in arbitrary dimensions. The worldlines are represented by link-based fluxes jx,μ ∈ {0, 1, ... q − 1} that obey constraints which enforce flux conservation modulo q at each site of the lattice. For vanishing field the worldlines of flux thus consist of closed loops with flux conservation modulo q. The magnetic field gives rise to additional sinks and sources for the worldlines, such that in this case also open worldlines with magnetic terms at their endpoints are possible. For the case of vanishing magnetic field we implement as suitable form of the Hodge decomposition and show that the constraints can be obeyed by switching to dual variables which are plaquettes of flux and defect lines of flux winding around the compactified directions of the lattice. For the special case of two dimensions the self-duality of the model at all q follows directly form this dual representation of the worldlines and one can re-derive the critical couplings βc = ln(1 + √ q) as a simple corollary. We finally also present a first exploratory assessment of the worldline representation for Monte Carlo simulations by discussing a suitable generalization of the worm algorithm. We analyze the worm algorithm in 2-d simulations where we provide a determination of the dynamical critical exponent z for the q = 2 and q = 4 cases which have second order phase transitions and for q > 4 address aspects of the worm algorithm for the emerging strong first order transitions. 2 Worldline representation for the Potts model In this section we will first derive the worldline form of the partition sum and discuss some observables in this representation which later will be the basis for simulations with the worm algorithm. For the case of vanishing magnetic field we will then introduce dual variables, i.e., plaquette variables and defect fluxes to resolve the constraints. This representation may be used for a local dual update that we will use for cross-checking the worm results. Finally we consider the 2-dimensional case where the dual representation can be used to establish full duality in the Kramers-Wannier sense [24] that allows one to rederive the (known) results for the critical values βc for all q. 2.1 Derivation of the worldline representation We consider a d-dimensional hypercubic lattice with periodic boundary conditions. The sites of the lattice are denoted by x and on each site we place a spin sx ∈ {0, 1, ... q − 1} where q ≥ 2 is a free integer-valued parameter. The partition sum of the q-state Potts model is given by Z = ∑ {s} e β ∑ x ∑d μ=1 δ(sx,sx+μ̂) + η ∑ x δ(sx,0) = ∑ {s} [∏ x,μ e β δ(sx,sx+μ̂) ] [∏ x e η δ(sx,0) ] , (1) where the parameter β in front of the ferromagnetic nearest neighbor term is the inverse temperature and η denotes an external magnetic field where the corresponding term favors the spin orientation sx = 0. x is summed over all lattice sites, μ runs over all d directions and μ̂ is the unit vector in direction μ. ∑ {s} = ∏ x 1 q ∑q−1 sx=0 is the sum over all spin configurations and we use δ(s, s ′) ≡ δs,s′ to denote the Kronecker delta. In the second step of (1) we have rewritten the sums in the exponent into products of individual exponentials.


Introductory remarks
The presentation of the Prokofev-Svistunov worm algorithm in [1] was a milestone for the developement of worldline Monte Carlo simulations in statistical mechanics and lattice field theory. Not only did worldline 2 representations and suitable algorithms considerably increase the numerical efficiency of Monte Carlo simulations, but also helped to solve some complex action problems, in particular those emerging in lattice field theories at finite chemical potential (for some reviews see, e.g., [2][3][4][5]).
It is important to note that there is no unique general approach to worldline techniques and suitable corresponding algorithms, but different types of systems have to be treated individually. Over the years new tools and techniques have been developed for worldline representations and now also examples of worldline representations for some non-Abelian symmetries are known, which pose a considerably more challenging task for a worldline approach.
In this paper we explore worldline techniques for yet another system with a non-Abelian symmetry, the q-state Potts model [8] with the permutation group of q elements as its symmetry group (see [9] for a review). In the context of lattice field theories mainly the 3-state Potts model, which is equivalent to the Z 3 spin model, has been analyzed [10] - [23]. Z 3 is the center group of SU(3), the gauge group 1 Member of NAWI Graz. 2 We point out that in the context of gauge theories the objects that correspond to worldlines are actually discretized worldsheets described by occupation numbers placed on plaquettes. For notational convenience we continue to use "worldlines" but remark that often very similar techniques (see, e.g., [6,7]) lead from worldlines for spins systems or lattice field theories to the corresponding "worldsheet" representations of gauge theories with the same symmetry group (which is of course a local symmetry in the gauge theory context).
of QCD, such that the 3-state Potts model serves as a simple model theory for studying aspects of the strong interaction.
We here use a worldline representation for the Potts model at arbitrary q that also allows for the inclusion of magnetic terms, and work out the representation in arbitrary dimensions. The worldlines are represented by link-based fluxes j x,µ ∈ {0, 1, ... q − 1} that obey constraints which enforce flux conservation modulo q at each site of the lattice. For vanishing field the worldlines of flux thus consist of closed loops with flux conservation modulo q. The magnetic field gives rise to additional sinks and sources for the worldlines, such that in this case also open worldlines with magnetic terms at their endpoints are possible.
For the case of vanishing magnetic field we implement as suitable form of the Hodge decomposition and show that the constraints can be obeyed by switching to dual variables which are plaquettes of flux and defect lines of flux winding around the compactified directions of the lattice. For the special case of two dimensions the self-duality of the model at all q follows directly form this dual representation of the worldlines and one can re-derive the critical couplings β c = ln(1 + √ q) as a simple corollary.
We finally also present a first exploratory assessment of the worldline representation for Monte Carlo simulations by discussing a suitable generalization of the worm algorithm. We analyze the worm algorithm in 2-d simulations where we provide a determination of the dynamical critical exponent z for the q = 2 and q = 4 cases which have second order phase transitions and for q > 4 address aspects of the worm algorithm for the emerging strong first order transitions.

Worldline representation for the Potts model
In this section we will first derive the worldline form of the partition sum and discuss some observables in this representation which later will be the basis for simulations with the worm algorithm. For the case of vanishing magnetic field we will then introduce dual variables, i.e., plaquette variables and defect fluxes to resolve the constraints. This representation may be used for a local dual update that we will use for cross-checking the worm results. Finally we consider the 2-dimensional case where the dual representation can be used to establish full duality in the Kramers-Wannier sense [24] that allows one to rederive the (known) results for the critical values β c for all q.

Derivation of the worldline representation
We consider a d-dimensional hypercubic lattice with periodic boundary conditions. The sites of the lattice are denoted by x and on each site we place a spin s x ∈ {0, 1, ... q − 1} where q ≥ 2 is a free integer-valued parameter. The partition sum of the q-state Potts model is given by where the parameter β in front of the ferromagnetic nearest neighbor term is the inverse temperature and η denotes an external magnetic field where the corresponding term favors the spin orientation s x = 0.
x is summed over all lattice sites, µ runs over all d directions andμ is the unit vector in direction µ.
is the sum over all spin configurations and we use δ(s, s ) ≡ δ s,s to denote the Kronecker delta. In the second step of (1) we have rewritten the sums in the exponent into products of individual exponentials.
Each of these exponentials can assume only two different values, and we can write them as where α is a real parameter. In the second step we have written the Kronecker delta as a sum over the q-th roots of unity. The weight factors w α j are given by We will use this representation of the exponentials both for the nearest neighbor terms where we set α = β, as well as for the magnetic terms where α = η and s = 0. Using (2) in the expression for the partition sum we obtain The last step for obtaining the worldline representation of the Potts model is to sum over the original spin variables in (4), where we have defined Thus we find the worldline representation of the partition sum in the form where we introduced the discretized divergence ∇ · j x = µ [j x,µ − j x−μ,µ ].
These constraints require that the sum of the divergence ∇ · j x and the magnetic variable m x at each site x has to vanish modulo q. This condition has the nice geometrical interpretation, that the flux j x,µ is conserved modulo q and that the magnetic variables can act as sources and sinks modulo q. In Fig. 1 we show for q = 5 possible examples of flux and magnetic variables at a site x.
The admissible configurations thus are worldlines of flux j x,µ that is conserved modulo q and the worldlines either form closed loops or open worldlines with magnetic terms m x at their endpoints which serve as sources or sinks (again modulo q). Below we will discuss strategies for Monte Carlo updates of the flux configurations subject to the constraints (9) based on the worm algorithm [1].
We conclude this section with providing the expressions for some bulk observables in terms of the flux and magnetic variables. More specifically we consider the internal energy E and the heat capacity C which can be obtained as first and second derivatives of ln Z with respect to β, as well as the magnetization M and the susceptibility χ which are the first and second derivatives with respect to η. These derivatives of ln Z can be evaluated directly in the worldline representation (8) and one finds after a few steps of algebra The expectation values on the right hand sides of (10) are expectation values in the worldline formulation.  (11)). The rhs. plot gives an explicit example of the fluxes for p x,µν = 1 in the q = 5 model. w β jx,µ and w β jx,µ denote the first and second derivatives of w β jx,µ with respect to β, and w η mx and w η mx the first and second derivatives of w η mx with respect to η.

Resolving the constraints with dual variables
For the case of vanishing magnetic field, i.e., for η = 0, one may develop the worldline representation further by introducing dual variables such that the constraints are automatically fulfilled. For η = 0 the magnetic sinks and sources are absent in the worldline representation and the admissible configurations of flux correspond to closed loops, where at every site flux is conserved modulo q. The dual variables then consist of plaquette variables and defect fluxes that can be combined to build up all admissible configurations of the flux variables. The plaquette variables are integers p x,µν ∈ {0, 1, ... q − 1} placed on the plaquettes which we label as (x, µ, ν) with 1 ≤ µ < ν ≤ d. A value p x,µν of a plaquette variable generates flux on the four links of the plaquette as follows: The configuration of flux that is generated by some p x,µν = 0 obviously corresponds to the smallest possible closed loop with flux conserved modulo q. In Fig. 2 we illustrate the loop using the general form of (11) in the lhs. plot and with the example p x,µν = 1 for q = 5 in the rhs. plot.
At a given link (x, µ) the corresponding flux j x,µ of course receives contributions from the plaquette variables of all plaquettes that contain (x, µ). To indicate that in (11) we only show the contribution of a single plaquette we used the notation with "←" in (11). The full contribution of all plaquettes to a given link will then be summarized below.
It is important to note that the flux generated by the plaquette variables does not generate all admissible flux configurations, since loops of conserved flux may also close around the periodic boundaries we use. This can be taken into account by introducing defect fluxes defined as and The S ν (x) are the support functions for the coordinate axes in direction ν through the origin, and the parameters f ν , ν = 1, 2, ... d determine the flux that is introduced on the corresponding loops along the coordinate axes that close around the periodic boundary conditions. We now may represent all admissible configurations of the flux j x,µ in the form In a more abstract language, the representation (13) is the Hodge decomposition [25,26] of all flux configurations with ( ∇ · j x ) mod q = 0 ∀ x. In terms of the dual variables the partition sum now reads where the j x,µ in the link weight factors w β jx,µ are computed using (13). In the dual form the partition function is a sum over configurations of the plaquette variables p x,µν and the parameters f µ for the winding flux. These degrees of freedom no longer need to obey any constraints. We will see below, that the form (14) can be used for a local Monte Carlo update of the systems, while the worm update will be implemented directly in the worldline form (8).

Self-duality in two dimensions
As a small corollary of the worldline/dual form of the q-state Potts model we have presented here, we re-derive the self-duality of the model in two dimensions [9]. We work directly in the infinite volume limit, where the contributions of the defect fluxes can be neglected. Furthermore, in 2-d we have only one type of plaquettes, p x,12 , and we can label 3 these plaquettes by their lower left corner x, i.e., in 2-d we may denote the plaquette variables as p The flux on the links is evaluated according to (13), which in the absence of defect fluxes and in 2-d reduces to The values j x,µ determine which of the factors w β jx,µ given in Eq. (3) is taken into account. Note that (3) distinguishes only two values w β 0 and w β if p x and p x−μ are equal and w β 1 if p x and p x−μ are different. Thus we may write the corresponding interaction term in the form where we have defined Thus we may write the 2-d partition sum also as a sum over configurations of the plaquette variables (an overall irrelevant constant has been dropped) Comparing this form of the partition sum in terms of plaquette variables with the defining form (1) in terms of spins, we find that the partition function is self-dual. The original coupling β and the dual couplingβ are related via the relation (17).
If one now assumes that there is only a single critical point at β c , then we may setβ c = β c and use (17) to determine β c , Thus we have re-derived the well-known result for the critical coupling β c at all values of q as a simple corollary of the worldline representation presented here.

Monte Carlo simulation with worldlines
In this section we report some first exploratory studies of using the worldline formulation for Monte Carlo simulations. We discuss the implementation of the worm algorithm for the worldlines of the qstate Potts model and check its correctness by comparing its results to simulations in the standard spin representation and to results of a local update of the dual representation. Subsequently we consider the q = 2 and q = 4 models, that both have second order transitions, and provide a first estimate of the corresponding dynamical critical exponents for the worm algorithm. Finally we briefly discuss properties of the worm algorithm also for the first order case (q > 4).

Worm algorithm for the q-state Potts model worldline representation
The worm algorithm for our worldline representation of the Potts model is a generalization of the Prokofev-Svistunov worm algorithm [1]. It starts at some randomly selected site x 0 of the lattice and then propagates through the lattice attempting to change the flux on each link it visits by ±1. The worm finishes when it reaches the starting point x 0 and all constraints are intact again. Each step of the worm is accepted with a Metropolis decision [28].
In pseudocode notation the steps of the worm can be formulated as follows (rand() denotes a uniformly distributed random number in the range [0, 1]) : • Randomly select an increment ∆ ∈ {−1, +1} • Randomly select a starting site x 0 and set x ← x 0 • Repeat until worm terminates: -Select a direction µ ∈ {±1, ±2, ... ± d} -Proposal flux:j x,µ = (j x,µ + sign(µ)∆) mod q (with j x,µ ≡ j x−| µ|,|µ| for µ < 0) -Compute ρ = w β jx,µ / w β jx,µ . If (rand() < ρ ) then j x,µ ←j x,µ and x ← x +μ -If (x = x 0 ) terminate worm It is straightforward to see that the worm generates only admissible flux configurations and is ergodic. In order to study the case of non-vanishing magnetic terms one may start the worm with the insertion of a magnetic term (source) and then in each step of the worm also offer to close the worm with a second insertion of a magnetic term (sink). This is the version we test below, but remark that also a variant where the worm inserts two magnetic terms connected by a "large hop" is an interesting alternative option, in particular if the magnetic field is very small (see, e.g., the discussion [18] for the special case of q = 3). Finally we remark that also other variants of the worm, such as directed worms, geometric worms or worms with heat-bath steps [29][30][31] could further extend the toolbox for worldline simulations of q-state Potts models.
Before we come to a more detailed assessment of Monte Carlo simulations of the worldline representation of the q-state Potts model we briefly discuss the evaluation of the new representation and the worm algorithm. In the top row of plots of Fig. 3 we compare the results for the energy density as a function of β from three different simulations based on the conventional spin representation (circles), the local dual representation (14) (blue triangles pointing down) and the worm algorithm discussed in this section which updates the worldline representation (magenta triangles pointing up). We compare the models for q = 4, q = 5 and q = 30 (left to right) to assess second order, weak first-order and strong first order behavior. We choose a small lattice of size 8 × 8 where finite volume effects are still strong. This allows us to test also the non-trivial topological aspects that play a role in the dual representation via the defect fluxes, and in the worldline representation through worms that wind around the periodic boundary conditions. Obviously the results of all three simulations agree very well for the whole range of β and the different q we consider.
A similar picture results for the comparison of local spin and worm simulations for the case of finite magnetic field shown in the bottom row of plots in Fig. 3. Here the worm was started and terminated with magnetic flux insertions as discussed above. Again we compare the energy density as a function of β for q = 4, q = 5 and q = 30 (left to right) on a small 8 × 8 lattice to assess finite volume effects. We use η = 0.0, η = 0.2 and η = 0.5 and again find very good agreement of local spin simulation results (crosses) with those from the worm algorithm (other symbols).

Exploratory assessment of the algorithm for q ≤ 4 (2 nd order transitions)
We start the analysis of the worm algorithm with exploratory studies in the q = 2 and q = 4 cases where the Potts model has second order phase transitions at the critical couplings β c = ln(1 + √ q). The key challenge for Monte Carlo simulations is critical slowing down near β c that can be quantified by the dynamical critical exponent z. With finite volume scaling techniques for simulations done at β c we here estimate z for the worm algorithm in the q = 2 and q = 4 cases using the energy as observable. Typical statistics for these determinations are 10 6 to 10 11 configurations, which reflects the demanding nature of the determination of dynamical critical exponents. For completeness we summarize the conventions of autocorrelation functions, autocorrelation times and the dynamical critical exponent z. The normalized autocorrelation function for the energy is where E i is the value of the energy at the i-th step of the Monte Carlo time series. The autocorrelation function decays exponentially for sufficiently large time separations k and we may use this exponential decay to define the exponential autocorrelation time τ exp . Summing the autocorrelation function defines We compare the data from the local spin simulation (crosses) to those from the worm algorithm (other symbols). Again the lattice size is 8 × 8 for assessing finite volume effects and we show results for η = 0.0, η = 0.2 and η = 0.5.
the integrated autocorrelation time τ int , such that the defining equations for the two autocorrelation times we use here are given by [32], where τ int ≤ τ exp . Note that in principle k max should be considered in the limit k max → ∞, but in a practical determination one needs to cut off the sum that defines τ int . Such a cut-off is considered safe when the summation obeys k max ≥ 6 τ int (k max ) [32,33].
In the critical region of second order transitions the autocorrelation time diverges as a power of the correlation length ξ, where z is the dynamical critical exponent. On a finite lattice of volume V = L d (here d = 2), the correlation length ξ cannot diverge and will instead be cut off by the linear extent of the lattice L, which is expressed in the second equality of (22). For comparing autocorrelation times to the Metropolis algorithm or the Swendsen-Wang cluster algorithm, where one update sweep visits all degrees of freedom of the lattice, we need to define proper units for the autocorrelation times. We take this into account by normalizing the measured autocorrelation times with the factor worm length /(2V ), i.e., the average fraction of the degrees of freedom that was subject to a worm update.
In the lhs. plots of Fig. 4 we show results for the normalized autocorrelation function A(k) as a function of the Monte Carlo time separation k. We compare the autocorrelation functions for different volumes L 2 and consider q = 2 (top) and q = 4 (bottom). The inserted plots show the same data using a logarithmic scale for the vertical axis. We observe exponential behavior of the correlation function up to values of k where the autocorrelation function becomes smaller than the statistical error from the Monte Carlo simulation.
From the results for the autocorrelation function shown in Fig. 4 we now can determine τ int and τ exp according to (21). The rhs. plots of Fig. 4 illustrate our determination of τ int : We show the integrated autocorrelation time τ int (k max ) as a function of the cut-off k max and again compare q = 2 (top) and q = 4 (bottom). Obviously, when considered as a function τ int (k max ) of the cut-off k max , the integrated autocorrelation time approaches a constant for large k max , which is an estimate for τ int . Thus we estimate τ int by using the large-k max values of τ int (k max ) self-consistently for all volumes with the condition k max ≥ 6τ int (k max ) and denote the results as τ int,cut .
We explore a second strategy for determining τ int by using a fit of τ int (k max ) in the low-k max range up to where τ int (k max ) starts to saturate to the form [33] τ int (k max ) = τ int − c exp(−k max /τ exp ). The fit parameters are τ int at k max = ∞, the exponential autocorrelation time τ exp and an irrelevant constant c. This leads to estimates of both, the integrated and the exponential autocorrelation time which we denote by τ int,f it and τ exp,f it , respectively.
In the lhs. plots of Fig. 5 we show our results for τ int,cut , τ int,f it and τ exp,f it as a function of the lattice extent L, using a double logarithmic plot (top: q = 2, bottom: q = 4). We observe a slight deviation from perfect linear behavior expected according to (22). We attribute these deviations to finite size corrections, and in order to take them into account we fit the data for τ from consecutive pairs of volumes V 0 = L 2 0 and V = L 2 with extents L 0 < L with (22). These fits give rise to results for the dynamical critical exponents z which we study as a function of the larger extent L.
On the rhs. of Fig. 5   for z decrease with increasing spatial extent L, eventually reaching saturation which indicates that with our largest volumes we are reaching the infinite volume extrapolated result for z and all three determinations consistently give rise to values of z ∼ 0.63. Since the dynamical critical exponent usually is determined from the integrated autocorrelation time, we here quote the values z int,cut = 0.6352 (2) and z int,f it = 0.62325(2) which we take from the largest volume pair (the error is the statistical error). Combining these results gives our final result z int = 0.629 (6), where we quote the distance of the average to the two different determinations as our systematical error (the statistical error is much smaller). For completeness we also quote our final result for the exponential autocorrelation time, z exp = 0.62895 (8), which agrees well with the integrated autocorrelation time (error here is again the statistical error).
When inspecting the q = 4 results for the dynamical critical exponent (bottom right of Fig. 5) plotted against L on a double logarithmic scale, we do not observe a systematic variation with the spatial extent L and the results from determinations based on τ int,cut , τ int,f it and τ exp,f it agree within a a relative error of 5%. Quoting again the results from the largest volume pair we find z int,cut = 1.3448(7) and Rhs.: Results for the dynamical critical exponents obtained form fits of τ according to (22) using pairs L 0 < L of subsequent lattice sizes.
z int,f it = 1.34076 (8), which combine to a value of z int = 1.342 (2). The result for the dynamical critical exponent determined from the exponential autocorrelation time is z exp = 1.3478 (3), and also for q = 4 we find very good agreement of the results for z.
In summary, for the simple worm algorithm we find critical dynamical exponents z = 0.629(6) for q = 2 and z = 1.342(2) for q = 4. Obviously the simple worm studied here is still outperformed by the corresponding Swendsen Wang cluster algorithms with z < 0.3 for q = 2 and z ∼ 1.0 for q = 4 [36]. This may be attributed to effects such as links being visited several times by the worm [29]. For improved worm techniques [29,30] dynamical critical exponents of z < 0.3 were reported for q = 2 and we expect that such an increase in performance with improved worms should be possible also for the q = 3 and q = 4 models. Concerning the dependence on q the increase from z ∼ 0.63 for q = 2 to z ∼ 1.34 for q = 4 indicates that as q is increased and one approaches the first order regime, the numerical simulations become more demanding. We remark that a similar observation was also made for the Swendsen-Wang cluster algorithm [36]. Rhs.: Average worm length as a function of the spatial extent L. We compare the results for q = 4 (second order behavior) with q = 10, q = 20 and q = 30 where the first order behavior becomes increasingly stronger.

Challenges for q > 4 simulations (1 st order transitions)
For q > 4 the 2-d Potts model has first order transitions at β c = ln(1 + √ q) and for reasonably large volumes Monte Carlo algorithms face a severe sampling problem due to mixing of phases at β c . This issue typically can be addressed with various reweighting and histogram techniques (see, e.g., [34,35]), which of course can also be implemented for worm algorithms updating the worldline representation. We here do not attempt a detailed and systematic analysis of worm algorithms in the first order regime of the Potts models, but only briefly report some simple observations we made in small test simulations at q > 4, partly using large q where the first order transition is very strong.
The sampling problem is illustrated in the lhs. plot of Fig. 6, where we show histograms for the energy density E/V computed at β c . We compare the histograms for q = 10, q = 20 and q = 30 (note that due to the mentioned sampling problems a smaller volume was used for the latter two) that display the characteristic double peak structure of first order transitions. We observe, that as q is increased and the transition becomes harder, the region between the peaks becomes wider, or -in more physical terms -the latent heat grows. Thus it is increasingly harder for an algorithm that samples the un-modified Boltzmann weight to cross between the two peaks and simulations near β c get stuck in only a part of configuration space.
It is interesting to see how the sampling problem affects the worm algorithm for the worldline representation. This is illustrated in the rhs. plot of Fig. 6 where we show for simulations at β c the average worm length as a function of L for different q on a double logarithmic scale. For the second order case (q = 4) the worm length keeps growing with L (at least up to the volumes we considered), while for q = 20 and q = 30 it starts to drop for the larger values of L which we attribute to the fact that between the peaks the distribution in the histogram becomes suppressed exponentially and obviously the worm becomes confined in only one of the peaks. For q = 10 we observe a behavior that is similar to the q = 4 case, probably because at q = 10 the first order transition is not so strong yet. However, one may expect that also for q = 10 the worm length will start to drop for sufficiently large L and we conclude that the well-known sampling problems near first order transitions manifest themselves in a considerable decrease of the average worm length which reflects insufficient sampling.

Summary and discussion
In this paper we have explored a worldline representation for the q-state Potts model with magnetic term in arbitrary dimensions. The worldlines are described by link-based flux variables and the flux is conserved modulo q. For non-zero magnetic field the worldlines can start and end in magnetic source and sink terms. These are absent for vanishing magnetic field and admissible configurations are closed loops of flux that is conserved modulo q. For this case we show that one may resolve the constraints by introducing dual variables that consist of flux around plaquettes and defect fluxes that wind around the periodic boundaries, i.e., we implement the Hodge decomposition. Finally we show that in the 2-d case one may use the dual representation to establish self-duality of the q-state Potts models as a simple corollary and re-derive the known critical couplings β c = ln(1 + √ q).
Having established the worldline and the dual representation we address possible Monte Carlo simulations in the new representations. We present a worm algorithm for updating the worldlines and verify it against a conventional simulation in terms of spins and a local simulation of the dual representation. These tests were done for both, vanishing and finite external field and not only check the correct implementation of the worm algorithm, but also the correctness of the worldline and the dual representations.
In an exploratory study we determine the dynamical critical exponent of the worm algorithm in two dimensions for the q = 2 and q = 4 cases. The corresponding results are z = 0.629(6) and z = 1.342 (2). A brief look at the models with q > 4 indicates that the sampling problem near the corresponding first order transitions leads to a decrease of the worm length at larger volumes. This suggests that also the worm becomes trapped in only a part of the configuration space and suitable reweighting and histogram techniques need to be used.
It is obvious that the numerical assessment of the worldline simulations in this paper has only an exploratory character and several future directions will be interesting to follow. In particular using more advanced worm techniques for studies of the q > 2 cases in dimensions larger than d = 2, as well as extended studies of the system with non-vanishing magnetic field terms are planned for future work.