{textblock*}100pt(399pt,-96pt) CERN-TH-2016-056 Sampling of General Correlators in Worm-Algorithm Based Simulations
Abstract
Using the complex -model as a prototype for a system which is simulated by a worm algorithm, we show that not only the charged correlator , but also more general correlators such as or , as well as condensates like , can be measured at every step of the Monte Carlo evolution of the worm instead of on closed-worm configurations only. The method generalizes straightforwardly to other systems simulated by worms, such as spin or sigma models.
Keywords: Monte Carlo; Worm algorithm; Dual variables; Flux-representation; General correlators.
1 Introduction
In a paper from 2001 [1], Prokof’ev and Svistunov proposed the worm algorithm as an alternative to cluster algorithms [2] to overcome critical slowing-down in the simulation of classical spin models like the Ising or the 3 state Potts model, but also of the complex model.
Critical slowing-down usually occurs at a second order phase transition where the typical length-scale over which the degrees of freedom (e.g. the spins of an Ising system) are correlated, diverges and the system develops long-range order. Local update algorithms then usually become very inefficient due to the high energy barrier that has to be overcome to flip individual spins against the field of its nearest neighbors and to the corresponding low acceptance rate of such updates. Cluster algorithms overcome this problem by being non-local: typical structures of correlated sites, the so called clusters, are grown and then flipped as a whole. Because the typical size of these clusters grows like the correlation length, critical slowing-down is averted. Unfortunately, the typical cluster size usually grows even further when leaving the critical region and going deeper into the ordered phase, resulting in a loss of efficiency of the cluster updates which becomes particularly dramatic when the typical cluster size approaches the size of the whole system111When the typical clusters consist of almost the whole system, the cluster updates essentially just flip the whole system back and forth.. The more relevant problem with cluster algorithms is however that they are in general not useful if the system couples to an external field, as the large energy barrier that has to be overcome to flip a whole cluster against the external field quickly leads to very low acceptance rates even for moderate cluster sizes.
Worm algorithms on the other hand are based on local updates and remain relatively efficient even in the ordered phase and in the presence of external fields. For bosonic systems, worm algorithms emerge quite naturally when expressing the partition function in terms of a particular kind of integer-valued "dual" variables, the so-called flux-variables (see [1, 3, 4, 5, 6, 10] and Sec. 1.1). The working principle behind this dualization process is as follows: consider the simple system described by the following partition function:
| (1.1) |
If we wanted to use Monte Carlo to compute for example expectation value and variance of the field , i.e.
| (1.2) |
the standard way to do this would be to interpret as probability weight for doing importance sampling with respect to . To obtain the dual formulation, we instead write the integrand of (1.1) in a power series in , carry out the integral for each monomial (which can usually be done analytically),
| (1.3) |
such that
| (1.4) |
and then use the as probability weights for doing importance sampling with respect to the monomial order , which is our new, integer-valued configuration variable, in terms of which the above observables read
| (1.5) |
This change of variables from to is of course only meaningful if the are real and non-negative for all . In the above mentioned dual representation of bosonic systems, the flux variables, which can be interpreted as representing the hopping of particles/excitations between neighboring sites, play a similar role as . However if the excitations carry a conserved charge, these flux variables are in general subject to a so-called closed loop constraint, which makes it impossible to update them with ordinary local Metropolis. The worm algorithm deals with this problem by temporarily violating the constraint on two sites, referred to as head and tail of the worm, and updates the flux variables while the head is moved around from site to neighboring site using importance sampling. The constraint is restored as soon as head and tail meet each other again. This reformulation of the partition function in terms of flux variables also avoids the sign problem in many significant cases, notably in those where the sign problem is related to the introduction of a non-zero chemical potential (see [3, 4, 9, 10, 15] and Sec. 1.1 below).
As pointed out for example in [5, 7, 8], the sampling of the location of the head of the worm (relative to its tail) during the worm update corresponds to the sampling of the correlator for the particle to whose hopping the updated flux variables correspond, and provides therefore an easy and very efficient way to estimate this correlator by recording a histogram for the head-tail distances realized during the worm updates.
In the present paper we show, using the complex -model as a prototype for systems which can be simulated by worm algorithms, that not only the correlator of the charged fields and (which are the excitations to whose hopping the flux variables correspond), but also more general correlators such as and as well as condensates like can be measured at every step of the Monte Carlo evolution of a worm. However, as the field is integrated out exactly (and with it the symmetry), the direct measurement of such one and two-point functions requires in general the introduction of non-vanishing source terms222As an alternative to source terms, one could also introduce a background field: for the complex model discussed in this paper, one would change variables where is the background field which would be set to the vacuum expectation value. to the action, to which, as suggested in [8], correspond additional dual variables representing so-called monomers (see Sec. 1.1). Note that the necessity of a source term in order to obtain a non-vanishing result for the condensate is not a flaw of the dual formulation but a consequence of the definition of the partition function in terms of an Euclidean path-integral, which explicitly sums over all possible field configurations. This means that even in the ordered phase where the minimum of the free energy develops a degeneracy and one usually assumes (motivated by the fact that in an infinite system, vacua that correspond to different minima of the free energy do not interact) that the system undergoes spontaneous symmetry breaking and picks just one of the corresponding degenerate vacua, the partition function explicitly always sums homogeneously over all of these vacua, unless one specifies a source term in the action that gives dominant weight to just one of them. This will be discussed in more detail in Sec. 2.2 and 2.5.
The paper is organized as follows: in the remainder of this section, we first define in part 1.1 our lattice action and the corresponding partition function, then we derive, following [9], the flux-variable representation of the latter. In part 1.2 we then discuss the relation between the charged correlator and the worm algorithm. Sec. 2 is intended to demonstrate how some simple modifications of the worm algorithm allow one to measure also more general correlators as well as condensates during the worm evolution, provided a non-vanishing source term is added to the action. The section also contains a short discussion on why such a source term is necessary. A short summary follows in Sec. 3. A detailed description of our simulation algorithm is given in appendix A.
1.1 Complex with Source Terms
The Euclidean lattice action for a complex field, coupled to a chemical potential and a (complex) source term (which could be thought of as an external magnetic field), reads in dimensions:
| (1.6) |
where and are dimensionless lattice parameters333We use the convention , , which is convenient when generalizing the action to the case. Note also that our is four times the ”standard ” used in other lattice formulations of theory.. As (1.6) is in general complex for non-zero values of the chemical potential , the corresponding partition function,
| (1.7) |
has a so-called sign problem, as a complex has no probabilistic meaning and can therefore not be used for importance sampling of configurations within a Monte Carlo simulation.
To overcome this problem we follow [9] and first write in (1.7) as a product of individual exponentials,
| (1.8) |
and then write these exponentials, except the ones that correspond to the potential part, in terms of their power-series,
| (1.9) |
| (1.10) |
and
| (1.11) |
| (1.12) |
After writing in polar form , the partition function (1.7) can then be written as
| (1.13) |
The non-negative intergers , and , are called flux and monomer variables, respectively: counts the number of hoppings of particles from site to site and of antiparticles from site to , while counts the corresponding inverse moves. It is therefore convenient [9] to introduce new variables and , such that
| (1.14) |
where counts the net charge flowing from site to site and counts the total number of particles and antiparticles hopping around between and . Similarly and count the particle and anti-particle content of site , and we define and as follows,
| (1.15) |
so that counts the total charge content of site and its total particle content. Integrating now over the in (1.13) and using (1.14) and (1.15) yields the final form of the partition function (up to an irrelevant constant pre-factor):
| (1.16) |
where
| (1.17) |
At first sight it seems that the source terms in (1.16) give rise to a complex phase. By writing them in the form
| (1.18) |
with the radial, , and the polar source, , we see that locally this is indeed the case. But due to the delta-function constraint in (1.16), we have 444Due to the delta function constraints in (1.16) we have . Furthermore we have , as each -variable appears twice in this sum with opposite sign, and therefore ., which implies , and there is therefore no net complex phase555Note that if and were treated as independent variables, (1.18) would not be true and the overall complex phase would not vanish in general.. All other terms in (1.16) are manifestly real and non-negative, and thus the flux representation (1.16) of the partition function (1.7) is sign-problem free and can be sampled by Monte Carlo.
1.2 Charged Correlator and Worm Algorithm
Due to the delta function constraints, configurations contributing to (1.16) cannot be sampled directly by a local Metropolis algorithm: the transition probabilities for local updates of the and variables would always vanish. For a non-zero source , one could sample the and variables simultaneously, such that the delta function constraint is always satisfied, but if is small, this would be highly inefficient due to low acceptance rates. For , in order to satisfy the delta function constraint, one would have to update randomly chosen closed chains of variables, which is in general also very inefficient. The alternative is to use a worm algorithm, which, roughly speaking, is based on the idea [1, 7, 8, 9] to update the variables while sampling configurations that contribute to the charged correlator instead of configurations that contribute to (1.16) itself, as will be explained in what follows.
1.2.1 Charged Correlator
In order to define correlation functions, we consider from now on the source in (1.16) as a local quantity , keeping in mind, that the partition function (or any observable) will actually only be evaluated when . The charged two-point function is then given by the following expression:
| (1.19) |
Carrying out the formal derivative of with respect to in (1.16), leads to (highlighting in red the changes from (1.16)):
| (1.20) |
which, after dividing by , can be interpreted as,
| (1.21) |
Analogously the formal derivative of with respect to yields
| (1.22) |
and for the mixed second derivative we find:
| (1.23) |
where, we have set as stated above. Although equations (1.21), (1.22) and (1.23) are correct, they are only well defined if . Also these expressions so far do not affect in any way the delta function constraint in (1.16) and will therefore not help in finding an efficient updating algorithm for the -variables.
To bring (1.20) into a more useful form, we can absorb the factor of into a shift of the and summation variables666To improve readability, we continue to highlight in red the important changes from one equation to the next.:
-
•
Terms in (1.20) for which , such that can, for , be written as
(1.24) By introducing new variables and , such that
(1.25) we have and obtain from (1.24):
(1.26) Note that for terms with , the right hand side of (1.25) would require , which is not valid as, according to (1.15), all -variables have to be non-negative. But these terms are anyway zero in (1.20) and it is therefore completely fine that they can not be expressed in terms of valid primed variables , .
- •
Dropping again the primes from the and variables, eq. (1.20) can therefore be written as
| (1.30) |
and in a similar way we find
| (1.31) |
and finally:
| (1.32) |
Equations (1.30), (1.31) and (1.32) are now well defined even if and we can deduce that in the flux representation formulation of the partition function, the right way to implement observables like
| (1.33) |
| (1.34) |
and
| (1.35) |
is[5, 6, 7, 8], to think of and as external source and sink, respectively, where the effect of adding to the system is a change in the local constraint at ,
| (1.36) |
and in the local weight at ,
| (1.37) |
Similarly we have for the insertion of :
| (1.38) |
and,
| (1.39) |
1.2.2 Worm Algorithm
The idea behind the worm algorithm is now as follows (cf. [1, 7, 8, 9]): first propose to insert at some site an external source/sink pair , . The in the delta function constraint at , coming from the insertion and the from the insertion of cancel each other and the insertion of the external source/sink pair changes therefore only the argument of at :
| (1.40) |
If this insertion is rejected, the proposal counts as an attempted worm update. If the insertion is accepted, the next step could either consist of proposing to remove the source/sink pair again from the system, or to move from to a randomly chosen neighboring site, e.g. . The latter update would change the product of local weights for the sites and as follows:
| (1.41) |
and in order to get a non-zero result on the right hand side of the arrow in (1.41), this update has to be combined with a simultaneous shift of the corresponding flux variable, (assuming for the moment that the variables are kept fixed) as depicted in Fig. 1. If this move of from to and the simultaneous shift are accepted, one can propose to move further to a random nearest-neighboring site of , and so on. In this way, one can update the -variables while sampling the charged correlator . Whenever, during this procedure, the head of the worm hits the site where also the worm’s tail is located, one can in addition to the head-shift also propose to remove the external source/sink pair again. If this latter proposal is accepted, the worm terminates: we have completed a whole worm-update and are left with a new closed-worm configuration which contributes to the partition function (1.16)777As an alternative to the insertion and removal of the external source/sink pair, one can do importance sampling with respect to the location of the pair and just reweight it away when evaluating observables on closed-worm configurations..
If is non-zero, one also has to include the update of the -variables into the worm algorithm in order to correctly sample the charged correlator. This can in principle be done by adding a new move to the worm-update algorithm, in which one replaces the shift in the -variable, , required to compensate for the in the delta functions in (1.41), by the two shifts and , which could be interpreted as compensating for the displacement from to of the external charge , by inserting a pair of oppositely charged (dynamical) monomers at and . However, in this case, there is no need to restrict the set of possible locations to which the head, , of the worm can be shifted to the nearest-neighboring sites of , and one can instead choose a random site as the target location (see Fig. 2). If such a move is accepted, the worm just continues from the new location as before proposing either to move the head to a nearest-neighboring site by changing a -variable, or to move the head to a random new site by changing two -variables. But still, the worm can only be terminated if the head hits again the site where its tail is located and the removal of the external source/sink pair is accepted.
The measurement of the two-point functions, e.g. of the component , happens by measuring the fractional Monte Carlo time during which the external source/sink pair is present in the system and located at sites and respectively.
So far, the algorithm is analogous to the closed-worm algorithm described in [11] for the case of a Potts model in the presence of a chemical potential and an external magnetic field. In the next section we supplement this algorithm with additional moves which allow for the sampling of more two-point functions as well as of one-point functions. For the latter the moves will be analogous to the open-worm moves of [11], but unlike [11], where simulations are performed using closed-worm or open-worm a priori, here the choice between closed-worm and open-worm will be performed stochastically.
2 Sampling General Correlators and Condensates During the Worm
In the previous section we have seen that the worm algorithm can be interpreted as sampling the two-point function . The latter can therefore be measured in a very efficient way during the worm-update by keeping track of the fractional Monte Carlo time the external source/sink pair is present in the system and source and sink are located at and respectively.
Other observables on the other hand, can only be measured on "closed-worm" configurations, obtained in between successive worm updates888A worm update consists of either a failed attempt to insert an external source/sink pair into the system, or, if the insertion succeeds, of all local updates between insertion and removal of the source/sink pair.. The auto-correlation times of such observables do not depend on the number of worm updates but on the number of "micro-steps"[7] done during these worm updates (i.e. the number of local updates of the flux variables). It therefore makes sense to choose the number of worm updates between successive measurements (on closed-worm configurations) so that the average number of local updates is of the order of the number of degrees of freedom in the system, which corresponds to the usual definition of a sweep used in ordinary non-worm based Monte Carlo simulations. As the average worm length999The average worm length is the average number of attempted local updates between start and end of a worm, i.e. between insertion and removal of the external source/sink pair. varies as a function of the simulation parameters (, , ), so does the optimal number of worm updates which are necessary to obtain a certain average number of local updates between successive measurements.
In Fig. 3 we show the integrated auto-correlation time for the average "energy"101010Note that we call the ”energy” just for simplicity: it is only one part of the full energy of the system.,
| (2.1) |
as a function of , once in units of worms111111The pre-factor is just a normalization: is the number of space-time dimensions, while is the number of independent components of , i.e. in the present case, we have , as a complex scalar field can be understood as an field. The factor reflects the fact that the worm can start with different choices for the head and reflects the number of different directions in which the worm can start on each site. (left) and once in units of sweeps (right), where we now define a sweep, as suggested above, as the number of worms which are necessary on average to proceed as many local updates as there are degrees of freedom in the system. In the first case, where the auto-correlation time is given in units of worm-updates, it can be seen that the number of such worms required to de-correlate two configurations decreases with increasing , and that in the disordered phase, is proportional to the system size, while in the ordered phase no volume-dependency is visible. This has of course to do with the fact that the average number of local updates required to complete a worm, increases with and becomes proportional to the system volume in the ordered phase (see Fig. 4). In the second case, where is measured in terms of sweeps, no volume dependency is visible, except at the pseudo-critical point. And even there, the volume dependency seems to be rather mild.
It is highly convenient to use the above definition of a sweep as measure for the time between successive measurements, because the computational cost for such sweeps is almost independent of (at least for a reasonably broad interval of values) and it can therefore also efficiently be used in combination with parallel tempering [13], where it avoids long waiting times among the replicas before swap moves can take place.
Note that although we tune the number of worms per sweep so that the average number of local updates or micro-steps per sweep is of the order of the system size, the sweep itself is defined in terms of a fixed number of worm updates and not in terms of a fixed number of micro-steps. By allowing the insertion and removal of external sources/sinks, we have extended the configuration space explored by the Markov process, so that it generates both, configurations that contribute to the usual partition function , as well as configurations that contribute to partition functions of the form , which describe the behavior of the system in the presence of an external source at and an external sink at for all possible . So the Markov process in fact samples configurations that contribute to the generalized partition function
| (2.2) |
But we want for example,
| (2.3) |
and not
| (2.4) |
and therefore, only if we define a sweep in terms of a fixed number of worm updates (which, since worm updates always start and end on closed-worm configurations, equals the number of configurations that contribute to ) instead of a fixed number of micro-steps, we get the correct result without having to determine the correction factor . For observables which depend on configurations in only, this correction factor is irrelevant. Also the masses obtained from fits to (2.3) and (2.4) would of course be identical. But for example the magnetic susceptibility obtained from (2.4), i.e.
| (2.5) |
would in general be wrong.
The introduction of (2.2) looks like an unnecessary complication. Indeed, in the standard worm algorithm (see e.g. [5]), one simply uses reweighting of to estimate . Here we need because, in the remainder of this section, we will further generalize to
| (2.6) |
where
| (2.7) |
and
| (2.8) |
This will allow us to measure also more general correlators like or , as well as the corresponding condensates and during the worm updates, instead of on closed-worm configurations only.
2.1 General Correlator in Terms of External Sources: New Worm Moves
To demonstrate how one can measure general correlators during the worm updates (without the use of reweighting), we use (1.18) to write the source terms in the partition function (1.16) in polar form:
| (2.9) |
We can then for example define the correlator for radial excitations by
| (2.10) |
where
| (2.11) |
The meaning of this correlator becomes clear when writing the original action (1.6) in terms of spherical coordinates :
| (2.12) |
Taking the derivative of the partition function with respect to brings therefore down a factor under the integral, which for all measures the magnitude of the projection of the field on the direction specified by . Analogously, taking the derivative of the partition function with respect would bring down a factor which, after dividing by , yields for all the magnitude of the projection of the field on the direction perpendicular to the direction specified by .
If we write
| (2.13) |
we recognize, that by using again shifts of the form (1.25), (1.28), we can write (2.11) as
| (2.14) |
Proceeding in the same way with the second derivative, we find:
| (2.15) |
Of course, (2.15) could have been obtained in a simpler way by using that
| (2.16) |
and therefore
| (2.17) |
and so on, but it is instructive to repeat the argument with the shifts of the - and -variables which give rise to the Kronecker deltas in the local constraints and local weights on the second and third line of (2.15).
In a similar way, one finds the correlator for tangential or angular excitations:
| (2.18) |
The first two pieces inside the bracket on the last line of (2.15) and on the right hand side of (2.18) can directly be measured during the worm update described above, no matter whether is zero or not. The other two pieces are non-zero only if and in order to measure them during the worm update, we have to extend the worm algorithm by an additional type of move which allows the worm to change the external source/sink at its head. This can be done by modifying the disconnected worm update described above in Fig. 2 in such a way, that instead of moving the external source from to some other site and shifting and , one can remove completely from the system, still compensate the absence of its charge at by shifting , but now, instead of inserting at again a and compensating for its charge by shifting , we insert a and shift (see Fig. 5). If this move is accepted, the worm samples from now on the last instead of the first piece from inside the brackets on the last lines of (2.15) and (2.18). The move which changes the external charge at the head of the worm back from to works completely analogously, except that the -variables have to be shifted in the opposite direction. The complex phases appearing in these expressions get always canceled by the complex phases coming from the shifts of the -variables during the head-changing moves. Note also that with these additional head-changing moves, the worm can still be terminated only if its head and tail are located on the same site, and in addition, if the head’s external charge is again the same as at the beginning of the worm.
The head-changing moves are completely general: if we had chosen for example the principal chiral model instead of complex , we could also have introduced moves that would allow us to sample correlators like , , etc. during the worm [15]. In [16] we used this method to verify the existence of oscillatory two-point functions in the three-state Potts model if the real and imaginary parts of the field are coupled to different external fields.
2.2 Measuring Condensates with a Worm: Why a Source Term is Needed
The partition function for our complex model is defined in terms of an Euclidean path integral which explicitly sums over all possible field configurations. This means that without an explicit breaking of the global symmetry (e.g. by a source term in the action), the condensates and have to vanish identically even in the ordered phase, as the sum of all the condensates corresponding to different -degenerate vacua, is zero. In the context of "ordinary" Monte Carlo simulations, based on local Metropolis or heat-bath updates of the original field variables, , this might not always be so obvious, as in such simulations, with increasing system size, the Markov chain usually has more and more problems to tunnel between different degenerate vacua. The reason for this is of course the local update strategy and the fact, that the different vacua should indeed become non-interacting in the thermodynamic limit. This means that any intermediate configuration required for the tunneling from one vacuum into another, becomes energetically more and more expensive with increasing system size and therefore more and more unlikely to be realized, so that the Markov chain tends to remain for a long time in the same vacuum, leading to apparently non-zero condensates. This is of course an illusion: it just reflects the fact that the algorithm is non-ergodic and therefore fails to generate a representative subset of configurations that contribute to the path integral.
In the flux-representation (1.16) the situation is qualitatively different: a non-zero source term is mandatory for a non-zero condensate, even for very large systems! The reason for this is the following: as we have integrated out the complex phases of the -field analytically in order to obtain the weights for the dual configurations in terms of flux variables, we also integrated out the global symmetry. Therefore every admissible flux-variable configuration, which could be understood to describe a configuration of particle world-lines, stands in fact for the superposition of the simultaneous realizations of these particle world-line configurations in each of the degenerate vacua. So, without adding source terms to the system, a spontaneous breaking of the global symmetry can only be observed in quantities which look the same in all degenerate vacua, or which are non-local in the sense that they depend on at least two sites
121212To be more precise: as we have integrated out the complex phase of the field to arrive at the partition function (1.16), we cannot even break the global symmetry by setting the source to a non-zero value. One manifestation of this is the fact that (1.16) does not depend on the phase of the external source but only on its magnitude , as shown at the end of Sec. 1.1. Another manifestation is that which can be seen from (1.21) and (1.22), together with the fact that for all valid configurations contributing to (1.16). However, by specifying a non-zero value for we can measure the invariant in polar coordinates (see also [15]).:
for example the value of the -variable between two sites, which measures how well the fields on these two sites are "in phase"131313Here ”in phase” effectively means ”aligned”, but for the case where one sums over all possible orientations in which the fields can align in the internal space, i.e. how dominant the "in phase" contribution to the integral over the individual phases of the fields on the two sites is. If the average value of the -variables increases, this indicates that the system gets more and more ordered.
By looking at equations (1.30) and (1.31), we can see that if (which implies ), there is no way to satisfy the local delta function constraints in (1.30) or (1.31) simultaneously for all sites: one can move the defect (introduced by the in the delta function at ) around by changing a sequence of -variables, but without the insertion of a second external source or sink at which this sequence of changed -variables could terminate (which would then correspond to measuring a two-point function instead of a condensate), there will always be a vanishing delta function factor somewhere in the weight of the configuration. So there exist no configurations which could support a non-zero measurement of a condensate if the source is set to zero.
However, if is non-zero, a second external source or sink is not required, as then the defect can be compensated by a shift of one of the -variables in the system. It is precisely this property which we can use to efficiently measure condensates (and therefore also the disconnected piece of correlators) even in the case of small sources (see Fig. 6): we propose for example to insert at some site a single external source or sink, or , and to compensate the corresponding appearing in the delta function at by an appropriate shift . If this insertion is accepted, we can either continue by proposing to remove the external source/sink again and to shift back to its original value, or we can propose to move the external source/sink to one of the neighboring sites of , say and to update in order to compensate for the resulting charge-displacement (completely analogous to the "connected move" in the "charged worm" described in the previous section). This procedure is repeated until the update which removes the external source/sink is accepted. Here, our algorithm is analogous to the open-worm algorithm from [11]. The expectation values and are then obtained by measuring the fraction of the Monte Carlo time, for which the external source/sink is present in the system.
Yet another advantage of this method for measuring condensates is, that it provides a simple means of variance reduction. One can absorb the summation over the -variables into the weight-function instead of sampling their values, i.e. one uses
| (2.19) |
in the worm algorithm instead of itself, where contains the sum over the and -variables, and the are modified Bessel functions of the first kind. The values of (2.19) can be pre-computed for a set of positive integers and and stored in a two-dimensional array141414More precisely: one can even pre-compute directly the relevant weight ratios and .. This speeds up the simulation and reduces the noise in the correlator and condensate measurements further.
In Appendix A, we provide pseudo-code describing our implementation of the concepts introduced in this and the previous sections, together with some comments and a derivation of the transition probabilities required to ensure detailed balance during the Monte Carlo simulation.
2.3 Gain in Efficiency
In figures 7 and 9 we compare results for the connected zero-momentum pieces of three different correlators, the charged one,
| (2.20) |
the one for radial,
| (2.21) |
and the one for angular excitations,
| (2.22) |
obtained with two different measurement techniques: once with the naive approach, based on measuring (1.21), (1.22) and (1.23) on closed-worm configurations, and once with the improved method, by measuring the correlators and condensates during the worm-updates according to (1.30), (1.31) and (1.32), as well as (2.15) and (2.18). The results from both methods were obtained during the same simulation and thus are directly comparable. For the naive approach, measurements were taken after every sweep (with a sweep defined as below eq. (2.1)). For the improved method, every sweeps, the histograms for the correlators and the condensates, accumulated during the worm-updates, were normalized by (which is the Monte Carlo time over which the histograms were accumulated), stored and reset. The simulations had a length of sweeps (plus thermalization). The errors were then obtained with the jack-knife method. As can be seen, the correlators obtained with the improved method provide a much cleaner and less noisy signal than the ones obtained with the naive procedure (based on the values of the and -variables in closed worm configurations). The large relative errors in the correlators shown in Fig. 9, are due to the subtraction of the disconnected piece in (2.20) and (2.21), which is large in the ordered phase. Overall, the new method yields smaller errors on the fitted mass, by a factor .
In Fig. 8 we show, for the two measurement methods, as a function of , the errors in the smallest component of the zero momentum piece of (2.21) (i.e. for ), before the disconnected piece is subtracted. In the disordered phase (), the improved measurement reduces the error in this observable by about an order of magnitude. In the ordered phase on the other hand (), the two methods seem to yield similar errors. The same is true for other components of the zero-momentum correlator, away from .
If one is just interested in extracting a mass from (2.21), the disconnected piece is irrelevant as it can be taken into account by simply adding a constant to the fitting function:
| (2.23) |
Still, in Fig. 10 we show how the reduced errors in the correlation functions, shown in Fig. 8, affect such a mass-fit: more precisely, Fig. 10 shows how the fitted mass depends on the fitting range. The zero-momentum correlation functions measured during a simulation contain contributions from several states which have different masses. They all decay exponentially like , which on a periodic lattice (and for neutral states) leads to terms like . In order to extract the correct ground-state mass by fitting (2.23) to the measured correlation function, one therefore either has to include additional -terms to take the excited states into account (which is usually difficult) or one has to restrict the fitting range to some interval in which the excited states have already died away sufficiently. Fig. 10 now shows how the mass , obtained by fitting (2.23) (without excited states) to the correlation function, changes as function of . For the improved correlator (black), it can be seen that the fitted mass approaches a clearly visible plateau with increasing . For the non-improved correlator (red), the plateau would of course also be there, but due to the bigger noise in this correlator measurement, the ground state signal gets quickly lost with increasing , sometimes even before the excited states become negligible, such that the plateau cannot be observed.
2.4 Alternatives and Possible Improvements
Despite the obvious improvement in the signal-to-noise ratio provided by the method described above in Sec. 2.1 for measuring different channels of general two-point functions during different stages of a generalized worm-algorithm, this measuring at different stages of the algorithm has the potential draw-back that the different channels of the two-point function are not necessarily determined with the same accuracy: with the current algorithm, the statistics in the measurement of any component of the two-point function is proportional to its expectation value. This is fine as long as there are not too many channels in the two-point functions or if one is only interested in the dominant propagating modes. However, if one is for example interested in the full mass spectrum of the theory, including sub-dominant modes/channels, this can be problematic. A naive attempt to deal with this problem would be to use reweighting to obtain always measurements for all channels during all stages of the worm-algorithm: if the worm is for example currently sampling , then measurements for , and can be obtained by appropriate reweighting factors, virtually replacing the external sources at head and tail of the worm (and appropriately shifting some -variables) and vice versa. It is however clear that there will in general be a big overlap problem between the different channels of the correlator.
The proportionality between statistics and true expectation value for the components of two-point functions which are sampled during the evolution of the worm, already leads to problems when dealing with just the standard two-point function, , if one is interested in very low temperatures or if the masses are all large, such that the correlator decays very quickly. To overcome this problem it was proposed in [6] to use an additional weight in the sampling of the two-point function , where can be thought of as an initial guess for , such that the algorithm has to sample only the ratio
| (2.24) |
which, for a good choice of , should vary much less as function of than itself, leading to an exponential improvement in the signal-to-noise ratio. This could of course be generalized to a set of functions , where and label the different possible choices for the external charges at head and tail of the worm. Furthermore, one could then use a method similar to the Wang-Landau sampling [14] in order to find the optimal , leading to equal statistics for all components and channels of the two-point function.
2.5 Effect of Source Term
As already mentioned in Sec. 2.2, the Euclidean path integral that defines the partition function for our theory sums over all possible field configurations and thereby integrates out the global symmetry. In the dual formulation in terms of flux variables this feature of the partition function is manifest as during the dualization process, the integration over the field has been carried out exactly. Spontaneous symmetry breaking can therefore only be observed in our simulations either when looking at non-local observables as for example two-point functions or the charge density, , or after the introduction of a non-vanishing source and the choice of an appropriate coordinate system for the internal space, in which one can then observe spontaneous symmetry breaking also in condensates.
In this section we now address the question of how big the source has to be chosen and how one can extract results for the theory at zero source from simulations carried out in the presence of a non-zero source. For this discussion, we will use the partition function (1.7) in terms of the original field variables in polar form:
| (2.25) |
In these variables the condensate for radial excitations (2.14) then reads
| (2.26) |
and the cosine in the observable under the integral makes it clear that this will only yield the expected result in the ordered phase if is sufficiently large in order to give dominant weight to the -configurations where when integrating over the phases . On the other hand, we would like to be as small as possible in order to reduce the associated bias. To get an approximate answer on how large should be chosen, we can proceed as follows: we know that in the ordered phase, the exponential of the interaction term ("interaction" in the sense of coupling fields on neighboring sites),
| (2.27) |
dominates over the entropy in the configuration space of the variables and the dominant contribution to the integral over these variables will come from configurations where all are parallel or, as we called it earlier, "in phase", i.e. we can assume that in the ordered phase, the partition function (2.25) behaves like
| (2.28) |
with being the modified Bessel functions of the first kind. If we now evaluate (2.26) by using (2.28), we find
| (2.29) |
where is the mean field solution to (2.29) and is the system volume. The ratio of Bessel functions in (2.29) vanishes if is zero and asymptotically approaches for as shown in the left panel of Fig. 11 as a function of .
In the thermodynamic limit , a non-zero value of leads to
| (2.30) |
no matter how small is. However, in practice we can simulate only finite systems and the thermodynamic limit has to be obtained by extrapolation (finite volume scaling). If we are interested in results for in an infinite system at , (2.29) reminds us that one first has to take the infinite volume limit before sending ; otherwise the result will always be zero, since for every finite volume there is a minimal value for below which will be in the region where is proportional to . The value of at which this change of behavior of happens is given by the peak in the quantity
| (2.31) |
when plotted as a function of and we can define the volume dependent lower bound for :
| (2.32) |
which scales like
| (2.33) |
as could have been expected since couples to an extensive quantity in the action (2.12). In the full theory would be determined by the maximum in the susceptibility
| (2.34) |
as a function of . The values of used for simulations that should yield results to be used in a scaling analysis to extract information on an infinite system at zero source, should therefore be chosen such that , where is the volume of the smallest system involved in the analysis.
Deep in the disordered phase on the other hand, entropy dominates over the exponential of the interaction term (2.27) and we can approximate the partition function (2.25) as follows:
| (2.35) |
For the condensate we then find
| (2.36) |
In the disordered phase the condensate is therefore expected to behave always like unless as there is now no volume factor in the argument of the Bessel functions. Due to this direct proportionality between condensate and source, it also follows that the susceptibility defined in eq. (2.34) should (to lowest order in ) be identical to the condensate.
In Fig. 12 we show for a (2+1)-dimensional system with , , , and for three different values of the chemical potential the condensate (left) and the corresponding susceptibility (right) as a function of the source . The top figures correspond to , which is in the disordered phase, where condensate and susceptibility are both proportional to the source and independent of the system size. The other figures correspond to (middle) and (bottom) which are both in the ordered phase (since ) where the condensate develops a plateau as soon as reaches a volume-dependent minimal value.
It should be mentioned at this point that if one is only interested in the magnitude of the condensate, an alternative to the introduction of a source for measuring it would be to use that
| (2.37) |
The same equation also holds for the correlator of radial excitations (2.15) (where for only the first two terms in (2.15) are non-zero):
| (2.38) |
The latter is particularly convenient: as we are usually simulating periodic finite systems, zero-momentum correlation functions always have a minimum, which is the point where the deviation of the value of the correlation function from the square of the corresponding condensate is minimal. For the zero-momentum piece of (2.38) the minimum is always located at , while for the zero-momentum piece of (2.37), its location can change as function of . However, since for vanishing source only the first two terms of (2.15) can be measured, while the other two terms (the ones that cannot be measured) would contribute similarly to the constant piece of the zero-momentum correlator, one has to keep in mind that the measured background should be multiplied by a factor in order to get the correct value for the squared condensate .
Another possibility would be to extract or from a fit of the form
| (2.39) |
or
| (2.40) |
to the zero-momentum pieces of the correlators in (1.32) and (2.15), respectively, so that and , where the factor in front of the first term in (2.40) is again due to the fact that for only the first two terms of (2.15) can be measured. In Fig. 13, we show an example of such measurements and compare them to results obtained by measuring the condensate directly, as described in Sec. 2.2, for four different, non-zero values of the source . Note that the slight oscillations in the data coming from the systems with vanishing source are not just due to systematic errors in the fitting process, but are a manifestation of the phenomenon described in [12], which is due to the quantization of the charge, visible in a finite system at low temperature (i.e. large ) (see also Fig. 14, right). In Fig. 14, we compare for two systems with different temporal extents, the result for the condensate obtained with the fitting method with the one obtained by extracting the condensate from the minimum (at ) of the zero-momentum piece of (2.15). The fitting method is more involved but works better for not too large , while for sufficiently large the two methods seem to work equally well.
3 Summary
Using the complex model as an example, we first reviewed the relation between the worm-algorithm and the charged correlator, how the former can be derived from the latter and how the standard worm can be thought of as sampling the charged correlator.
We then showed how the worm-algorithm can be generalized in order to sample also more general correlators and condensates and explained why this requires the inclusion of a non-vanishing source term in the action. In general, all correlators probing different internal space components of the field can be sampled during the worm Monte Carlo evolution by generalizing the standard worm-algorithm along the lines we presented.
In the Appendix we explain in some detail how our algorithm works and how the transition probabilities are computed in order to satisfy detailed balance.
Appendix A Worm Algorithm
In this appendix we describe our implementation of the generalized worm-algorithm introduced in the text. We explain why the code is organized in this particular way and show how the transition probabilities are obtained as to satisfy detailed balance.
A.1 Organization of Program
The main components of our program are described in Alg. 1 in terms of pseudo-code. They consist of 3 main routines:
-
•
a routine "WORM_UPDATE(x0,…)", which updates the and (and if the source is non-zero, also the and ) variables while sampling the full two-point function, starting and ending at "x0",
-
•
a routine "COND_WORM_UPDATE(x0,…)", which, if the source is non-zero, also updates the , , and variables while sampling the and condensates, starting at "x0" and ending at some other site (which will then be used as new value for "x0"),
-
•
and the routine "SWEEP()", which, for a fixed number of times, randomly executes either "WORM _UPDATE(x0,…)" or "COND_WORM_UPDATE(x0,…)" (provided the source is non-zero), or picks a new random location for "x0".
In "SWEEP()", the probability for choosing next a new random location for "x0" is always 1/2. If (or equivalently ) is set to zero, the probability for executing next "WORM_UPDATE(x0,…)" is also 1/2. If is non-zero, "WORM_UPDATE(x0,…)" and "COND_WORM_UPDATE(x0,…)" share the 1/2 probability for being executed next and are therefore executed with probability 1/4 each. The two worms can always start in two different ways: for "WORM_UPDATE(x0,…)" we can either have as the head and as the tail of the worm or vice versa, and similarly for "COND_WORM _UPDATE(x0,…)", we can either attempt to start by inserting a or a (together with an appropriate shift of a -variable, i.e. the insertion of an appropriate dynamical monomer of opposite charge).
The and variables are updated along with the and variables in "WORM_UP-DATE(x0,…)" and "COND_WORM_UPDATE(x0,…)": at every step of the worm the type of variable (, , or variable) to be updated next is selected at random, in such a way that the probability of selecting an variable ( variable) is the same as the probability of selecting a variable ( variable). In for example [9] another updating strategy was used, where only the variables were updated during the worm evolution while the variables were updated in separate sweeps, executed periodically after a fixed number of worms. This can become problematic as soon as the system develops long-range order and the average worm-length becomes of the order of the system size: the worm then evolves for a long time in a fixed -background, which slows down de-correlation and gives rise to larger errors in the two-point functions measured during the worm evolution. Furthermore, alternating between different types of updates in a fixed order strictly speaking breaks detailed balance (although this is usually harmless). Our updating scheme avoids both of these problems: the and variables do not form a fixed background while the worm evolves and detailed balance is satisfied at every step of the simulation.
Note that if one is interested in two-point functions only, one could in principle avoid the removal and re-insertion of the external source/sink pair (and the uniform sampling of ) whenever the worm closes, and instead just do importance sampling with respect to the location of the pair, which would be slightly more efficient as it would require less calls to the random-number generator. Measurements of observables which have to be obtained on closed-worm configurations would then have to be corrected by a reweighting factor, compensating for the presence of the external source/sink pair at .
Note also that "WORM_UPDATE(x0,…)" and "COND_WORM_UPDATE(x0,…)" could in principle be combined into a single worm, which, provided is non-zero, would not only sample the location and type of the head of the worm, but also sample the state of the worms tail: whether it should consist of an external source or sink as in "WORM_UPDATE(x0,…)" or if it should be replaced by an appropriate shift of the -variable (i.e. by a dynamical monomer) at the location of the tail as is the case in "COND_WORM_UPDATE(x0,…)". The reason why we decided to use separate worms is that our full simulation program is implemented to work for arbitrary linear and non-linear sigma models (complex corresponds to the linear -case and the also mentioned chiral-effective model to the non-linear -case) with arbitrary source terms, in which case the implementation with two different worms for sampling two- and one-point functions turned out to be simpler.
A.2 Detailed Balance and Transition Probabilities
In order to ensure that our Markov chain really samples the partition, one- and two-point function, we have to ensure detailed balance between each pair of successive configurations and , i.e. configurations contributing to the partition, one- or two-point function, which can be turned into each other by a single update. The detailed balance equation takes the general form
| (A.1) |
where is the weight of the configuration and is the transition probability for going from configuration to . However, as during different stages of our simulation, the number of possible moves changes, it can happen, that the move is chosen with a different probability than the inverse move . We must then factor out from the transition probabilities the different move-choice probabilities, , and write
| (A.2) |
and similarly
| (A.3) |
and then use only the reduced transition probabilities, , for the Metropolis acceptance test, i.e.
| (A.4) |
and
| (A.5) |
The reason for this is, that in order for the Metropolis test to take place, the corresponding move has to be chosen already and the probability for this to happen should therefore not be taken into account a second time for the final decision whether the transition to the new configuration should take place or not. This is of course always the case, but if the move choice probabilities are the same for all the moves that are possible during the simulation, in particular for all pairs of move and inverse move, the move choice probabilities just cancel out in (A.4) and (A.5).
The diagram in figure 15 illustrates the different move-choice probabilities which occur at different stages of the simulation and in order to demonstrate how the corresponding reduced transition probabilities are determined, we show the derivation of these probabilities for the moves occurring in "WORM_UPDATE(…)". In what follows we use the abbreviation , and for the variable we have that if , and otherwise151515The value of is different for and because ”WORM_UPDATE(x0,…)” is executed with different probabilities within ”SWEEP()” in these two cases, as explained at the beginning of Sec. A.1..
-
•
The detailed balance equation for starting the worm at site , i.e. inserting an external source/sink pair at this point, reads:
(A.6) from which it follows that:
(A.7) -
•
If the head of the worm consists of a , we define , and if the head is a , we set instead.
-
–
For moving the head from site to its nearest-neighbor in positive direction, say to , and changing , the detailed balance equation reads:
(A.8) from which it follows that if :
(A.9) and if :
(A.10) -
–
When moving the head from site to a nearest-neighbor in a negative direction, e.g. and shifting , the corresponding detailed balance equation becomes:
(A.11) such that for :
(A.12) and if :
(A.13) -
–
Moving the head from site to another random site and shifting , leads to the detailed balance equation,
(A.14) such that
(A.15) where
(A.16) and
(A.17) - –
-
–
-
•
Finally, if head and tail are again located both on site and consist of opposite external charges, one can propose to terminate the worm, for which the detailed balance equation is given by
(A.18) and therefore:
(A.19)
The reduced transition probabilities for the moves in "COND_WORM_UPDATE(…)" can be obtained in a a completely analogous way.
&sites is a reference to the two-dimensional array "sites", where:
sites[][0]=,
sites[][1]=sites[x][0] (for a model with sites[][1] would be different from sites[][0]),
sites[][2]=, sites[][3]=,
&bonds is a reference to the 3-dimensional array "bonds", where for nu<d, with =nu+1:
bonds[][nu][0]=,
bonds[][nu][1]=,
and for nud, with =(nu-d)+1:
bonds[][nu][0]=&bonds[][nu-d][0],
bonds[][nu][1]=&bonds[][nu-d][1],
&xnbr is a reference to the nearest-neighbor lookup table "xnbr", where for nu<d, with =nu+1, we have xnbr[][nu]=, and for nu>=d, with =(nu-d)+1: xnbr[][nu]=.
&wh is a reference to the array "wh" of weight ratios:
wh[n]=, where is defined in eq. (1.17),
h= is the hopping parameter,
mu= is the chemical potential,
sr= is the radial source,
Ns: number of sites per spatial dimension,
Nt: number of sites in the temporal direction,
d: number of dimensions,
sweeplen: number of attempted worms or shifts of x0 per sweep, can be tuned to obtain sweeps with fixed average number of local updates.
&ni: reference to the variable "ni", containing the total charge ,
&corrs and &corrt are references to the variables "corrs" and "corrt", containing the spatial and temporal (zero-momentum) propagators, e.g. corrt[t].
The notation x%N for x, N, refers to and the function randi(N) returns a random integer in the range while rand() is understood to return a random floating point number in the semi-open interval .
References
- [1] N. Prokof’ev, B. Svistunov, Worm Algorithm for Classical Statistical Models, Phys. Rev. Lett. 87 (2001) 160601, arXiv:cond-mat/0103146
- [2] R. H. Swendsen, J.-S. Wang, Nonuniversal Critical Dynamics in Monte Carlo Simulations, Phys. Rev. Lett. 58 (1987) 86
- [3] M. G. Endres, Avoiding the boson sign problem at finite chemical potential, PoS (LATTICE 2006) 133, arXiv:hep-lat/0609037
- [4] M. G. Endres, Method for simulating O(N) lattice models at finite density, Phys. Rev. D 75 (2007) 065012, arXiv:hep-lat/0610029
- [5] U. Wolff, Simulating the All-Order Strong Coupling Expansion I: Ising Model Demo, Nucl. Phys. B 810 (2009) 491-502, arXiv:0808.3934 [hep-lat]
- [6] U. Wolff, Simulating the All-Order Strong Coupling Expansion III: O(N) sigma/loop models, Nucl. Phys. B 824 (2010) 254-272, arXiv:0908.0284 [hep-lat]
- [7] T. Korzec, I. Vierhaus, U. Wolff, Performance of a worm algorithm in theory at finite quartic coupling, Comp. Phys. Comm. 182 (2011) 1477, arXiv:1101.3452 [hep-lat]
- [8] C. Gattringer, T. Kloiber, Spectroscopy in finite density lattice field theory: An exploratory study in the relativistic Bose gas, Phys. Lett. B 720 (2013) 210–214, arXiv:1212.3770 [hep-lat]
- [9] C. Gattringer, T. Kloiber, Lattice study of the Silver Blaze phenomenon for a charged scalar field, Nucl. Phys. B869 (2013) 56 arXiv:1206.2954 [hep-lat]
- [10] F. Bruckmann, C. Gattringer, T. Kloiber, T. Sulejmanpasic, Dual lattice representation for and models with a chemical potential, Phys. Lett. B 749 (2015) 495-501, arXiv:1507.04253 [hep-lat].
- [11] Y. Delgado, H. G. Evertz, C. Gattringer, Worm algorithms for the 3-state Potts model with magnetic field and chemical potential, Comp. Phys. Commun. 183 (2012) 1920-1927, arXiv:1202.4293 [hep-lat]
- [12] D. Banerjee, S. Chandrasekharan, Finite size effects in the presence of a chemical potential: A study in the classical non-linear O(2) sigma-model, Phys. Rev. D 81 (2010) 125007, arXiv:1001.3648 [hep-lat]
- [13] C. J. Geyer, Markov chain Monte Carlo maximum likelihood, Computing Science and Statistics: Proc. 23rd Symp. Interface (1991) 156-163
- [14] A. G. Cunha-Netto, A. A. Caparica, S. Tsai, R. Dickman, D. P. Landau, Improving Wang-Landau sampling with adaptive windows, Phys. Rev. E 78 (2008) 055701, arXiv:0808.1210 [cond-mat.statmech]
- [15] T. Rindlisbacher, P. de Forcrand, Lattice simulation of the chiral model at zero and non-zero pion density, PoS (LATTICE 2015) 171, arXiv:1512.05684 [hep-lat]
- [16] O. Akerlund, P. de Forcrand, T. Rindlisbacher, Oscillating propagators in heavy-dense QCD, arXiv:1602.02925 [hep-lat]