Sampling of General Correlators in Worm Algorithm-based Simulations

Using the complex $\phi^4$-model as a prototype for a system which is simulated by a worm algorithm, we show that not only the charged correlator $<\phi^{*}(x)\phi(y)>$, but also more general correlators such as $<|\phi(x)||\phi(y)|>$ or $<\text{arg}(\phi(x))\text{arg}(\phi(y))>$, as well as condensates like $<|\phi|>$, 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.


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 φ 4 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 system. 1 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. the standard way to do this would be to interpret e φ+λ φ 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), such that Z = ∞ n=0 w n (λ), (1.4) and then use the w n (λ) as probability weights for doing importance sampling with respect to the monomial order n, which is our new, integer-valued configuration variable, in terms of which the above observables read ∂ log(Z) ∂λ = n 1 + λ and ∂ 2 log(Z) ∂λ 2 = n 2 − n 2 − n (1 + λ) 2 . (1.5) This change of variables from φ ∈ R to n ∈ N 0 is of course only meaningful if the w n (λ) are real and non-negative for all n. 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 n. 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 φ 4 -model as a prototype for systems which can be simulated by worm algorithms, that not only the correlator φ * (x)φ(y) of the charged fields φ and φ * (which are the excitations to whose hopping the flux variables correspond), but also more general correlators such as |φ(x)||φ(y)| and arg(φ(x)) arg(φ(y)) 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 U(1) symmetry), the direct measurement of such one and two-point functions requires in general the introduction of nonvanishing source terms 2 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 U(1) 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 Secs. 2.2 and 2.5.
The paper is organized as follows: in the remainder of this section, we first define in Sec. 1.1 our lattice φ 4 action and the corresponding partition function, then we derive, following [9], the flux-variable representation of the latter. In Sec. 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.

Complex φ 4 with source terms
The Euclidean lattice action for a complex φ 4 field, coupled to a chemical potential μ and a (complex) source term s (which could be thought of as an external magnetic field), reads in d dimensions: (1.6) where κ and λ are dimensionless lattice parameters. 3 As (1.6) is in general complex for non-zero values of the chemical potential μ, the corresponding partition function, has a so-called sign problem, as a complex D[φ] e −S[φ] 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 e −S[φ] 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, and 3 We use the convention φ = 1 , which is convenient when generalizing the action to the O(N ) case. Note also that our κ is four times the "standard κ" used in other lattice formulations of φ 4 theory.
, the partition function (1.7) can then be written as The non-negative integers ξ x,ν , ξ x,ν and m x , m x are called flux and monomer variables, respectively: ξ x,ν counts the number of hoppings of particles from site x to site (x +ν) and of antiparticles from site (x +ν) to x, while ξ x,ν counts the corresponding inverse moves. It is therefore convenient [9] to introduce new variables k x,ν and l x,ν , such that ξ x,ν −ξ x,ν = k x,ν ∈ Z and ξ x,ν +ξ x,ν = |k x,ν | + 2 l x,ν , l x,ν ∈ N 0 , (1.14) where k x,ν counts the net charge flowing from site x to site (x +ν) and (|k x,ν | + 2 l x,ν ) counts the total number of particles and antiparticles hopping around between x and (x +ν).
Similarly m x and m x count the particle and anti-particle content of site x, and we define p x and q x as follows, so that p x counts the total charge content of site x and (|p x | + 2 q x ) its total particle content. Integrating now over the θ x in (1.13) and using (1.14) and (1.15) yields the final form of the φ 4 partition function (up to an irrelevant constant pre-factor): 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 with the radial, r s = √ 2 |s|, and the polar source, θ s , we see that locally this is indeed the case. But due to the delta-function constraint in (1.16), we have x p x = 0, 4 which implies x e i θ s p x = 1, and there is therefore no net complex phase. 5 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.

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 p and k variables would always vanish. For a non-zero source r s = √ 2|s|, one could sample the k and p variables simultaneously, such that the delta function constraint is always satisfied, but if |s| is small, this would be highly inefficient due to low acceptance rates. For |s| = 0, in order to satisfy the delta function constraint, one would have to update randomly chosen closed chains of k 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 k variables while sampling configurations that contribute to the charged correlator φ * (x) φ(y) instead of configurations that contribute to (1.16) itself, as will be explained in what follows.

Charged correlator
In order to define correlation functions, we consider from now on the source s in (1.16) as a local quantity s x , keeping in mind, that the partition function (or any observable) will actually only be evaluated when s x = s ∀x. The charged two-point function is then given by the following expression: Carrying out the formal derivative of Z with respect to s * x in (1.16), leads to (highlighting in red (color online) the changes from (1.16)): 4 Due to the delta function constraints in (1.16) we have p x = − ν (k x,ν − k x−ν,ν ). Furthermore we have x,ν (k x,ν − k x−ν,ν ) = 0, as each k-variable appears twice in this sum with opposite sign, and therefore x p x = 0. 5 Note that if s and s * were treated as independent variables, (1.18) would not be true and the overall complex phase would not vanish in general.

∂Z ∂s
, (1.20) which, after dividing by Z, can be interpreted as, Analogously the formal derivative of Z with respect to s y yields 1 Z ∂Z ∂s y = 1 s and for the mixed second derivative we find: where, we have set s x = s y = s as stated above. Although equations (1.21), (1.22) and (1.23) are correct, they are only well defined if |s| = 0. 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 k-variables.
To bring (1.20) into a more useful form, we can absorb the factor of ( 1 2 (|p x | − p x ) + q x ) into a shift of the p x and q x summation variables 6 : • Terms in (1.20) for which p x ≥ 0, such that ( 1 2 (|p x | − p x ) + q x ) = q x can, for q x > 0, be written as By introducing new variables p z and q z , such that we have |p x | = |p x | − 1 and obtain from (1.24): 6 To improve readability, we continue to highlight in red the important changes from one equation to the next.
. (1.26) Note that for terms with q x = 0, the right hand side of (1.25) would require q x = −1, which is not valid as, according to (1.15), all q-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 p x , q x . • Terms in (1.20) with p x < 0, for which 1 2 (|p x | − p x ) = |p x |, can be written as , (1.27) and by choosing p z and q z in order to satisfy (1.28) such that this time |p x | = |p x | + 1, we obtain from (1.27) again Dropping again the primes from the p and q variables, eq. (1.20) can therefore be written as W λ |p z | + 2 q z +δ x,z + ν |k z,ν | + |k z−ν,ν | + 2 (l z,ν + l z−ν,ν ) , (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 |s| = 0 and we can deduce that in the flux representation formulation of the φ 4 partition function, the right way to implement observables like is [5][6][7][8], to think of φ(x) and φ * (x) as external source and sink, respectively, where the effect of adding φ(x) to the system is a change in the local constraint at x, (1.36) and in the local weight at x, Similarly we have for the insertion of φ * (x): (1.38) and,

Worm algorithm
The idea behind the worm algorithm is now as follows (cf. [1,[7][8][9]): first propose to insert at some site x 0 an external source/sink pair φ(x 0 ), φ * (x 0 ). The −1 in the delta function constraint at x 0 , coming from the φ(x 0 ) insertion and the +1 from the insertion of φ * (x 0 ) cancel each other and the insertion of the external source/sink pair changes therefore only the argument of W λ at x 0 : 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 x 0 to a randomly chosen neighboring site, e.g. x = x 0 +ν. The latter update would change the product of local weights for the sites x 0 and x 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, k x,ν → k x,ν + 1 (assuming for the moment that the p variables are kept fixed) as depicted in Fig. 1. If this move Fig. 1. Starting from the configuration on the left-hand side, a standard worm-update consists of proposing to shift φ * (x) from x to some random neighboring site, x +ν, say, and to compensate for the charge-displacement by increasing the flux-variable k x,ν , as depicted in the right-hand figure. of φ * from x 0 to x and the simultaneous shift k x,ν → k x,ν + 1 are accepted, one can propose to move φ * further to a random nearest-neighboring site of x, and so on. In this way, one can update the k-variables while sampling the charged correlator φ(x 0 )φ * (x) . Whenever, during this procedure, the head of the worm hits the site x 0 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). 7 If |s| is non-zero, one also has to include the update of the p-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 k-variable, k x,ν → k x,ν + 1, required to compensate for the ±1 in the delta functions in (1.41), by the two shifts p x → p x + 1 and p x+ν → p x+ν − 1, which could be interpreted as compensating for the displacement from x to x +ν of the external charge φ * , by inserting a pair of oppositely charged (dynamical) monomers at x and x +ν. 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 x, and one can instead choose a random site y 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 k-variable, or to move the head to a random new site by changing two p-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 φ(x)φ * (y) , happens by measuring the fractional Monte Carlo time during which the external source/sink pair is present in the system and located at sites x and y 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.

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 1 x ∂s y . 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 x and y respectively.
Other observables on the other hand, can only be measured on "closed-worm" configurations, obtained in between successive worm updates. 8 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 length 9 varies as a function of the simulation parameters (κ, λ, s), 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", 10 as a function of μ, once in units of 2 d N worms 11 (left) and once in units of sweeps (right), where we now define a sweep, as suggested above, as the number L sweep 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 8 A 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. 9 The 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. 10 Note that we call E the "energy" just for simplicity: it is only one part of the full energy of the system. 11 The pre-factor 2 d N is just a normalization: d = 2 + 1 is the number of space-time dimensions, while N is the number of independent components of φ, i.e. in the present case, we have N = 2, as a complex scalar field can be understood as an O(2) field. The factor N reflects the fact that the worm can start with N different choices for the head and 2 d reflects the number of different directions in which the worm can start on each site. With increasing μ, the average worm-length increases, such that the average number of worms required to de-correlate two configurations decreases. In the right-hand figure, τ int is measured in units of sweeps, where a sweep is defined as the number of worm updates, such that the total number of local updates during the worm updates is of the order of the number of degrees of freedom in the system. In these units, the integrated auto-correlation time seems to depend only weakly on the volume, even at the pseudo critical points around μ ∼ 0.8. Fig. 4. The figure shows for a (2 + 1)-dimensional system of spatial size N s = 8, with couplings κ = 0.25, λ = 7.5 and source r s = 0.01, the average worm length divided by the system volume as a function of μ for four different values of N t ∈ {24, 32, 64, 128}. As can be seen, for μ > 0.8, when the system is in the ordered phase (broken symmetry), the different curves almost coincide, which means that the worm length is proportional to the volume. The latter is not true in the symmetric phase (unbroken symmetry), but as the worms are much shorter in this phase, the different curves cannot be distinguished on the given scale.
can be seen that the number of such worms required to de-correlate two configurations decreases with increasing μ, and that in the disordered phase, τ int 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 τ int 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 Z, as well as configurations that contribute to partition functions of the form Z 2 (x, y) = ∂ 2 Z ∂s * x ∂s y , which describe the behavior of the system in the presence of an external source φ at x and an external sink φ * at y for all possible (x, y). So the Markov process in fact samples configurations that contribute to the generalized partition function and not 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 Z) instead of a fixed number of micro-steps, we get the correct result without having to determine the correction factor Z gen /Z. For observables which depend on configurations in Z 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.
would in general be wrong. The introduction of Z gen (2.2) looks like an unnecessary complication. Indeed, in the standard worm algorithm (see e.g. [5]), one simply uses reweighting of Z 2 (x, x) to estimate Z. Here we need Z gen because, in the remainder of this section, we will further generalize to and This will allow us to measure also more general correlators like φ * (x)φ * (y) or φ(x)φ(y) , as well as the corresponding condensates φ and φ * during the worm updates, instead of on closed-worm configurations only.

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: We can then for example define the correlator for radial excitations by (2.11) The meaning of this correlator becomes clear when writing the original action (1.6) in terms of spherical coordinates (r x , θ x ): Taking the derivative of the partition function with respect to r s,z brings therefore down a factor r z cos(θ z − θ s ) under the integral, which for all θ z measures the magnitude of the projection of the field on the direction specified by θ s,z . Analogously, taking the derivative of the partition function with respect θ s,z would bring down a factor r s,z r z sin(θ z − θ s,z ) which, after dividing by r s,z , yields for all θ z the magnitude of the projection of the field on the direction perpendicular to the direction specified by θ s,z .
If we write Proceeding in the same way with the second derivative, we find: Of course, (2.15) could have been obtained in a simpler way by using that and therefore and so on, but it is instructive to repeat the argument with the shifts of the pand q-variables which give rise to the Kronecker deltas in the local constraints and local weights on the second and third lines of (2.15).
In a similar way, one finds the correlator for tangential or angular excitations: 5. New worm moves: for |s| > 0, starting from the configuration on the left-hand side, the (necessarily disconnected) head-changing worm-update consists of proposing to remove φ * (x) from x, and to insert instead φ(y) at some random other site y. The resulting charge-imbalance in the system and at x and y is compensated by increasing both, p x and p y by 1 (blue circle), as depicted in the right-hand figure. The inverse move works analogously: propose to replace φ(y) by φ * (x) and to simultaneously decrease p x and p y by 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 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 r s is zero or not. The other two pieces are non-zero only if r s > 0 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 x to some other site y and shifting p x → p x + 1 and p y → p y − 1, one can remove φ * (x) completely from the system, still compensate the absence of its charge at x by shifting p x → p x + 1, but now, instead of inserting at y again a φ * and compensating for its charge by shifting p y → p y − 1, we insert a φ and shift p y → p y + 1 (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 p-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 p-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 SU(2) principal chiral model instead of complex φ 4 , we could also have introduced moves that would allow us to sample correlators like π − (x)π 0 (y) , π + (x)σ (y) , 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.

Measuring condensates with a worm: why a source term is needed
The partition function for our complex φ 4 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 U(1) 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 U(1)-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, φ x , 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 U(1) 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 U(1) 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 12 : for example the value of the k-variable between two sites, which measures how well the fields on these two sites are "in phase", 13 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 k-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 |s| = 0 (which implies p x = 0 ∀x), 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 ±1 in the delta function at x) around by changing a sequence of k-variables, but without the insertion of a second external source or sink at which this sequence of changed k-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 |s| 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 p-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 x a single external source or sink, φ or φ * , and to compensate the corresponding 12 To 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 U(1) symmetry by setting the source s to a non-zero value. One manifestation of this is the fact that (1.16) does not depend on the phase θ s of the external source but only on its magnitude r s = √ 2|s|, 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 x p x = 0 for all valid configurations contributing to (1.16). However, by specifying a non-zero value for r s we can measure the U(1) invariant | φ | = | φ * | in polar coordinates (see also [15]). 13 Here "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. Fig. 6. The left-hand figure shows the situation where the insertion of a φ * (x) at some site x (indicated by the red circle) and a simultaneous shift, p x → p x − 1, (indicated by the filled blue dot) has been accepted. One can now propose either to remove φ * again from the system and to shift p x back to its original value, or one can propose to move φ * (x) to the site x +ν and to shift k x,ν → k x,ν + 1, as depicted in the middle figure. If the latter move is accepted one can again propose to either shift φ * (x +ν) further to a neighboring site of x +ν, or to remove φ * from the system and shifting p x+ν → p x+ν + 1, such that we would be left with the situation in the right-hand figure. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) ±1 appearing in the delta function at x by an appropriate shift p x → p x ∓ 1. If this insertion is accepted, we can either continue by proposing to remove the external source/sink again and to shift p x back to its original value, or we can propose to move the external source/sink to one of the neighboring sites of x, say x +ν and to update k x,ν 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 q-variables into the weight-function instead of sampling their values, i.e. one uses (2.19) in the worm algorithm instead of W λ (|p| + 2 q + A) itself, where A contains the sum over the k and l-variables, and the I ν (x) are modified Bessel functions of the first kind. The values of (2.19) can be pre-computed for a set of positive integers A and p and stored in a two-dimensional array. 14 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. . Both data sets were obtained from the same simulation. For the top row correlators, measurements were taken after every sweep and for the bottom row correlators the histograms were stored and reset after every 1000 sweeps. The error bars were then obtained with the jack-knife procedure (note that the correlators are displayed on a log-scale, whereas each error-bar is on a linear scale, which is set by the log-re-scaling factor of the corresponding data point). The dotted lines correspond to exponential fits which lead to the mass-values displayed in the upper right corner of each plot. Note the reduction in the error promoted by the new method.

Gain in efficiency
In Figs. 7 and 9 we compare results for the connected zero-momentum pieces of three different correlators, the charged one,  the improved method, every N sweeps = 1000 sweeps, the histograms for the correlators and the condensates, accumulated during the worm-updates, were normalized by 1/(L sweep × N sweeps ) (which is the Monte Carlo time over which the histograms were accumulated), stored and reset. The simulations had a length of 5 × 10 6 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 p and q-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 ∼ 5.
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 t = N t /2), before the disconnected piece is subtracted. In the disordered phase (μ < 0.8), the improved measurement reduces the error in this observable by about an order of magnitude. In the ordered phase on the other hand (μ > 0.8), the two methods seem to yield similar errors. The same is true for other components of the zero-momentum correlator, away from t = N t /2.
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 C 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 ∼e −m i t , which on a periodic lattice (and for neutral states) leads to terms like ∼ cosh(m i (t − N t /2)).
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 cosh-terms to take the excited states into account (which is usually difficult) or one has to restrict the fitting range to some interval [t min , N t − t min ] in which the excited states have already died away sufficiently. Fig. 10 now shows how the mass m, obtained by fitting (2.23) (without excited states) to the correlation function, changes as function of t min . For the improved correlator (black), it can be seen that the fitted mass approaches a clearly visible plateau with increasing t min . 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 t min , sometimes even before the excited states become negligible, such that the plateau cannot be observed.

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 φ * (x)φ(y) , then measurements for φ(x)φ(y) , φ * (x)φ * (y) and φ(x)φ * (y) can be obtained by appropriate reweighting factors, virtually replacing the external sources at head and tail of the worm (and appropriately shifting some p-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 twopoint functions which are sampled during the evolution of the worm, already leads to problems when dealing with just the standard two-point function, φ * (x)φ(y) , 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 ρ(x − y) in the sampling of the two-point function φ * (x)φ(y) , where ρ(x − y) can be thought of as an initial guess for φ * (x)φ(y) , such that the algorithm has to sample only the ratio , (2.24) which, for a good choice of ρ(x − y), should vary much less as function of (x − y) than φ * (x)φ(y) itself, leading to an exponential improvement in the signal-to-noise ratio. This ρ(x − y) could of course be generalized to a set of functions ρ f 1 ,f 2 (x − y), where f 1 and f 2 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 ρ f 1 ,f 2 (x − y), leading to equal statistics for all components and channels of the two-point function.

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 U(1) 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 φ x 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, 1 x k x,d , 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 r s is sufficiently large in order to give dominant weight to the φ-configurations where θ x ∼ θ s when integrating over the phases θ x . On the other hand, we would like r s to be as small as possible in order to reduce the associated bias. To get an approximate answer on how large r s 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 θ x variables and the dominant contribution to the integral over these variables will come from configurations where all θ x 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 with I ν (x) being the modified Bessel functions of the first kind. If we now evaluate (2.26) by using (2.28), we find where r is the mean field solution to (2.29) and V is the system volume. The ratio of Bessel functions in (2.29) vanishes if r s Vr is zero and asymptotically approaches 1 for (r s Vr → ∞) as shown in the left panel of Fig. 11 as a function of r s . In the thermodynamic limit (V → ∞), a non-zero value of r s leads to no matter how small r s 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 r in an infinite system at r s = 0, (2.29) reminds us that one first has to take the infinite volume limit before sending (r s → 0); otherwise the result will always be zero, since for every finite volume there is a minimal value for r s below which (r s Vr) will be in the region where I 1 (r s Vr) I 0 (r s Vr) is proportional to r s . The value of r s at which this change of behavior of I 1 (r s Vr) as could have been expected since r s couples to an extensive quantity in the action (2.12). In the full theory r lb s (V ) would be determined by the maximum in the susceptibility as a function of r s . The values of r s 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 r s r lb s (V min ), where V min 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: For the condensate we then find In the disordered phase the condensate is therefore expected to behave always like r ∝ r s unless r s 1 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 χ r defined in eq. (2.34) should (to lowest order in r s ) be identical to the condensate. In Fig. 12 we show for a (2 + 1)-dimensional system with κ = 0.25, λ = 7.5, N s = 8, N t = 24, 64, 128, 256 and for three different values of the chemical potential μ the condensate r (left) and the corresponding susceptibility χ r (right) as a function of the source r s . The top figures correspond to μ = 0.686, which is in the disordered phase, where condensate and susceptibility are both proportional to the source r s and independent of the system size. The other figures correspond to μ = 0.857 (middle) and μ = 1.2 (bottom) which are both in the ordered phase (since μ c ≈ 0.8) where the condensate develops a plateau as soon as r s 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 The same equation also holds for the correlator of radial excitations (2.15) (where for |s| = 0 only the first two terms in (2.15) are non-zero): The latter is particularly convenient: as we are usually simulating periodic finite systems, zeromomentum 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 N t /2, while for the zero-momentum piece of (2.37), its location can change as function of μ. However, since for vanishing source r s = √ 2|s| = 0 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 2 in order to get the correct value for the squared condensate r 2 .
Another possibility would be to extract |φ| or r from a fit of the form to the zero-momentum pieces of the correlators in (1.32) and (2.15), respectively, so that |φ| = C 1 and r = C 2 , where the factor 1/2 in front of the first term in (2.40) is again (left) and N t = 1024 (right) the condensate r as a function of μ, determined once by identifying r 2 with twice (see explanation in text below eq. (2.38)) the value of the zero-momentum piece of the measured (2.15) (black), and once by fitting to the same correlator the function (2.40) (red). As can be seen, for N t = 256 (left), the black and the red curve do not agree very well: the black curve shows too large values as the temperature is still too large and the zero-momentum correlator at N t /2 deviates too much from the constant background. The red curve from the fit shows however already an oscillatory signal. For N t = 1024 (right), the agreement between the two curves is much better and the oscillations as a function of μ are very pronounced. We have increased the number of sampling points between μ = 0.8 and μ = 0.9 in order to resolve also some of the fast oscillations (see inset). Between μ = 0.8 and μ = 1.2, the charge density increases from zero to slightly above one, which happens in N d−1 s = 8 2 = 64 small steps. Each valley in the condensate corresponds to one of these small steps. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) due to the fact that for r s = 0 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 r s . 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 N t ) (see also Fig. 14, right). In Fig. 14, we compare for two systems with different temporal extents, the result for the condensate r obtained with the fitting method with the one obtained by extracting the condensate from the minimum (at t = N t /2) of the zeromomentum piece of (2.15). The fitting method is more involved but works better for not too large N t , while for sufficiently large N t the two methods seem to work equally well.

Summary
Using the complex φ 4 model as an example, we first reviewed the relation between the wormalgorithm 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 nonvanishing 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 Algorithm 1 in terms of pseudo-code. They consist of 3 main routines: • a routine "WORM_UPDATE(x0,...)", which updates the k and l (and if the source is nonzero, also the p and q) 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 k, l, p and q 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 |s| (or equivalently r s ) is set to zero, the probability for executing next "WORM_UPDATE(x0,...)" is also 1/2. If |s| 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 p-variable, i.e. the insertion of an appropriate dynamical monomer of opposite charge).

Algorithm 1
The following pseudo-code describes our simulation program. The routine "SWEEP()" makes use of a worm-update routine "WORM_UPDATE()" which samples the two-point functions, allows for monomer insertions and head-changes of the worm (provided sr>0). If sr>0, it makes also use of a second worm-update routine "COND_WORM_UPDATE()" which samples the condensates. The input parameters are as follows: &sites is a reference to the two-dimensional array "sites", where: &wh is a reference to the array "wh" of weight ratios: wh[n]=W λ (n + 2)/W λ (n), where W λ (n) is defined in eq. (1.17), h=κ is the hopping parameter, mu=μ is the chemical potential, sr=s r 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 x k x,4 , &corrs and &corrt are references to the variables "corrs" and "corrt", containing the spatial and temporal (zero-momentum) propagators, e.g.  The l and q variables are updated along with the k and p variables in "WORM_UP-DATE(x0,...)" and "COND_WORM_UPDATE(x0,...)": at every step of the worm the type of variable (k, l, p or q variable) to be updated next is selected at random, in such a way that the probability of selecting an l variable (q variable) is the same as the probability of selecting a k variable (p variable). In for example [9] another updating strategy was used, where only the k variables were updated during the worm evolution while the l 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 l-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 l and q 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 x0) 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 x0.
Note also that "WORM_UPDATE(x0,...)" and "COND_WORM_UPDATE(x0,...)" could in principle be combined into a single worm, which, provided |s| 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 p-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 O(N ) sigma models (complex φ 4 corresponds to the linear O(2)-case and the also mentioned SU(2) chiral-effective model to the non-linear O(4)-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 C and C , 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 w(C)P C → C = w C P C → C , (A.1) where w(C) is the weight of the configuration C and P (C → C ) is the transition probability for going from configuration C to C . However, as during different stages of our simulation, the number of possible moves changes, it can happen, that the move C → C is chosen with a different probability than the inverse move C → C. We must then factor out from the transition probabilities the different move-choice probabilities, p, and write and similarly and then use only the reduced transition probabilities, P r , for the Metropolis acceptance test, i.e. 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 Fig. 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 A x = |p x | + 2 q x + ν (|k x,ν | + |k x−ν,ν | + 2 (l x,ν + l x−ν,ν )), and for the variable n c we have that n c = 1 if r s = √ 2|s| = 0, and n c = 2 otherwise. 15 15 The value of n c is different for r s = √ 2|s| = 0 and r s = √ 2|s| = 0 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 x 0 , i.e. inserting an external source/sink pair at this point, reads: from which it follows that: • If the head of the worm consists of a φ * , we define = 1, and if the head is a φ, we set = −1 instead. • For moving the head from site x to its nearest-neighbor in positive direction, say to x = x +ν, and changing k x,ν → k x,ν + =: k x,ν , the detailed balance equation reads: and if |k x,ν | < |k x,ν |: (A.10) • When moving the head from site x to a nearest-neighbor in a negative direction, e.g.
x = x −ν and shifting k x ,ν → k x ,ν − =: k x ,ν , the corresponding detailed balance equation becomes: • Moving the head from site x to another random site x and shifting p x → p x + =: p x , p x → p x − =: p x leads to the detailed balance equation, if |p x | < |p x |.
(A. 16) and and therefore: The reduced transition probabilities for the moves in "COND_WORM_UPDATE(...)" can be obtained in a completely analogous way.