Simulation strategies for the massless lattice Schwinger model in the dual formulation

The dual form of the massless Schwinger model on the lattice overcomes the complex action problems from two sources: a topological term, as well as non-zero chemical potential, making these physically interesting cases accessible to Monte Carlo simulations. The partition function is represented as a sum over fermion loops, dimers and plaquette-surfaces such that all contributions are real and positive. However, these new variables constitute a highly constrained system and suitable update strategies have to be developed. In this exploratory study we present an approach based on locally growing plaquette-surfaces surrounded by fermion loop segments combined with a worm based strategy for updating chains of dimers, as well as winding fermion loops. The update strategy is checked with conventional simulations as well as reference data from exact summation on small volumes and we discuss some physical implications of the results.


Introduction
Calculations at finite density are considered to be one of the great challenges for lattice field theory Monte Carlo simulations. The quest for an ab-initio understanding of non perturbative features of quantum field theories, such as phase diagrams at non-vanishing chemical potentials, have given rise to a considerable number of innovative approaches to the problem and reviews at the yearly lattice conferences [1] provide an overview about ideas that were pursued in recent years. At the core of the challenge is the so called "complex action problem" or "sign problem": Many theories have a complex action S when a chemical potential, or also a topological term (theta term, vacuum term) are coupled. In these cases the Boltzmann factor exp(−S) in the path integral has a complex phase and cannot be interpreted as a probability in a Monte Carlo simulation.
An elegant approach that entirely overcomes the complex action problem is to rewrite a lattice field theory in terms of new degrees of freedom, often referred to as "dual variables". In a successful dualization the theory is exactly rewritten in terms of the dual variables such that in the new representation the partition sum has only real and positive contributions. The dual variables are world-lines for matter fields and world-sheets for gauge fields. The world-sheets can either form closed surfaces or surfaces with a boundary formed by a matter loop. For bosonic abelian theories an impressive body of work based on the dual approach was presented in recent years (for examples see [2,3]).
For relativistic fermions the situation is complicated further by additional signs from the Grassmann nature of the fermion fields and the commutators of the Dirac matrices. Thus only relatively few real and positive dual representations for systems with fermions can be found in the literature, often in the limit of strong coupling or the sign quenched approximation (for some examples see, e.g., [4,5,6,7]). A case where a real and positive representation is known for arbitrary couplings is 2-dimensional QED, i.e., the Schwinger model, with chemical potential and a topological term [8]. Both, the chemical potential, as well as the topological term lead to complex action problems and the dual representation in terms of world-lines, dimers and world-sheets given in [8] solves both these problems in principle by providing an exact representation where all contributions to the partition sum are real and positive.
In this exploratory paper we present strategies for simulating the dual representation [8]. More specifically we consider the following two cases: 1) the one flavor case with a topological term (at vanishing chemical potential) and, 2) the case of two oppositely charged flavors with a non-zero chemical potential (no topological term coupled). In both cases one has to update the fermion loops together with the world-sheets that are bounded by the world-lines. In addition one has to sample the dimer contributions from the fermion integral which do not couple to the gauge degrees of freedom. The fermion loops and dimers constitute a highly constrained system and here we propose and test strategies to update them in accordance with the constraints coming from the gauge symmetries of the conventional representation. We test our algorithms against exact calculations on small lattices and conventional simulations without topological term or chemical potential. The results presented here constitute the first simulation of the Schwinger model at arbitrary, i.e., weak couplings with a finite chemical potential and a topological term.
We expect that the explorative study of simulation strategies presented here will be useful for simulating the highly constrained world-line/world-sheet representations of other fermionic systems. Furthermore, the reference data for the Schwinger model which can be obtained from the dual approach tested here, will be useful also for assessing other techniques that are explored for overcoming complex action problems: In particular methods based on complexification (see, e.g., the reviews [9]) or new strategies for simulating systems at finite vacuum term [10] can be cross-checked against results from the world-line/world-sheet methods presented here.
As announced in the introduction, in this paper we study two variants of the massless Schwinger model on the lattice: the one flavor case with a topological term, and the case of two oppositely charged fermions with chemical potentials. This section serves to present the two models in their conventional representation and to summarize the dual formulation presented in [8], which is the basis for the dual simulation discussed here.

The one flavor model with topological term
For the one flavor case the partition sum is given by where U ν (n) = exp(iA ν (n)), A ν (n) ∈ [−π, π] are the U(1)-valued link variables for the gauge fields. By n = (n 1 , n 2 ) we denote the sites of a 2-dimensional N S × N T lattice and for a technical simplification in the dual representation (see [8]) we restrict ourselves to lattices where N T and N S are multiples of 4. The index ν = 1, 2 runs over the Euclidean space (ν = 1) and time (ν = 2) directions, and V = N S N T is used to denote the total number of lattice points. For the fermions we use one-component Grassmann valued variables ψ(n) and ψ(n). In the conventional representation (1) For the gauge action S G [U ] we use the Wilson form, where the plaquettes U p (n), labelled by the coordinate n of their lower left corner, are the products U p (n) = U 1 (n) U 2 (n +1) U 1 (n +2) U 2 (n) and β is the inverse gauge coupling. The topological charge Q[U ] is introduced using the field theoretical definition which in the continuum limit (β → ∞) goes over into the integer valued topological charge of the continuum. The massless staggered fermions are described by the action with the staggered sign function γ ν (n) given by γ 1 (n) = 1 , γ 2 (n) = (−1) n 1 . We can write the partition function as where we have introduced the abbreviations η = β 2 − θ 4π , η = β 2 + θ 4π . The conventional representation (6) has a complex action problem for non-zero values of the vacuum angle θ, since then η = η , such that the Boltzmann factor has a complex phase and cannot be used as a probability weight in a Monte Carlo update. We remark, that the physically interesting continuum limit is reached for β → ∞ where both parameters η and η are positive. As we will see below it is sufficient to restrict β to values β > |θ|/2π, which ensures a particularly simple form of the dual representation.
To overcome the complex action problem an exact rewriting of the partition function in terms of new degrees of freedom, so-called dual variables, was presented in [8]. In this dual representation each term in the partition sum is real and positive, such that one can access finite values of θ with Monte Carlo simulations. The transformation is based on an expansion of the Boltzmann factors and a subsequent exact integration of the original field variables. The dual form of the partition sum is given by (compared to the form given in [8] we have dropped an overall irrelevant factor of 2 −V ) The sum {l,d,p} in (7) runs over all admissible configurations of the dual variables, non-intersecting oriented fermion loops l, dimers d and plaquette occupation numbers p(n). A configuration is admissible when each site of the lattice is either run through by a fermion loop or is the endpoint of a dimer. The fermion loops introduce fluxes along their contour. These fluxes have to be compensated by filling the fermion loop with plaquette occupation numbers p(n) ∈ Z. A value p(n) = +k (or −k) with k > 0 gives rise to k units of flux around the plaquette at site n with mathematically positive (negative) orientation. If two neighboring plaquettes are not separated by a segment of fermion loop, then they must have the same plaquette occupation numbers, such that the flux along the link where they touch is compensated.
In the dual representation all boundary conditions are periodic. An example of an admissible configuration on an 8×8 lattice is shown in Fig. 1. Note that we could also increase or decrease all plaquette occupation numbers by the same integer and still have an admissible configuration. Each configuration (l, d, p) comes with a weight factor n I |p(n)| (2 √ ηη) ( η/η ) p(n) that depends only on the plaquette occupation numbers p(n) which are restricted by the admissibility condition of the configuration. By I p we denote the modified Bessel functions, and obviously the weights are real and positive (assuming that we choose β > |θ|/2π such that η and η are positive). Thus the dual representation solves the complex action problem completely.
Note that a non-zero vacuum angle θ gives different weights to positive and negative plaquette occupation numbers via the factor ( η/η ) p(n) . For θ > 0 we have η < η and thus η/η < 1, which implies that for θ > 0, negative plaquette occupation numbers p(n) have a larger weight and thus are favored (and vice versa).
In this study we focus on bulk observables which can be obtained as derivatives of ln Z. In particular we compute the plaquette expectation value U p and its susceptibility χ p , as well as the topological charge density q and the topological susceptibility χ t , These derivatives can easily be evaluated also in the dual representation and give rise to weighted sums of the modified Bessel functions and their moments.  [8]). The lattice is filled with oriented self-avoiding loops (single lines with arrows) and dimers (double lines), such that each site is either run through by a loop or is the endpoint of a dimer. Flux introduced by loops has to be compensated by placing plaquette occupation numbers, which are specified by the integers we write inside each plaquette.

Two flavor model with chemical potential
The conventional representation of the lattice partition sum for the two flavor model with chemical potential reads (we omit the topological term here) ψ, ψ and χ, χ are two fermion flavors with opposite charge. Like the measure D[ ψ, ψ] given in (2), also the measure D[χ, χ] for the second flavor χ is a product over Grassmann measures at all sites. The fermionic action now contains chemical potentials µ ψ and µ χ for both flavors: Since the second flavor χ has negative charge, in S χ [U, χ, χ] the link variables U ν (n) are interchanged with their complex conjugate U ν (n) * . This ensures overall electric neutrality of the two flavor system as required by Gauss' law. The chemical potentials µ ψ and µ χ for the two flavors give different weights to temporal (ν = 2) forward and backward hops. We remark that for the two flavor case considered here physics will only depend on the sum of the two chemical potentials µ ψ + µ χ due to Gauss' law. However, for more than two flavors (keeping overall electric neutrality) physics will depend also on non-trival combinations of the individual chemical potentials (see, e.g., [11]). The dualization for more than two flavors is straightforward [8], and with this possible generalization in mind we find it instructive to explicitly show the dependence on the individual chemical potentials.
For the dual representation of the two flavor model we now have a second set of dual variables, l, d for the fermion loops and the dimers of the second fermion flavor. Together with l and d for the first flavor and the plaquette occupation numbers p(n) this constitutes the set of dual variables of the two flavor case. The dual form of the partition function of the two flavor case is then given as a sum over the configurations of all dual variables (compared to [8] we again drop an irrelevant overall factor) In an admissible configuration the fermion constraints have to be obeyed for both flavors, i.e., for both, the l, d and the l, d variables each site has to be either the endpoint of a dimer or run through by a fermion loop. At each link the combined flux of the fermion loops from both flavors has to be compensated by activated plaquettes. Since the second flavor has negative charge the flux from the l-loops is counted with a negative sign, i.e., here the flux from the loop l and from the plaquettes p(n) have to be equal and run in the same direction for cancellation. Also equal fluxes from l and l along a link that run in the same direction saturate the gauge constraints on that link. The rhs. of Fig. 2 shows such a double fermion loop (we use dashed lines for the l loops and the d dimers).
The two fermion flavors interact with each other only on those plaquettes which are connected to fermion loops of both flavors, such that the weight I |p(n)| (β) has a different p(n) from what it would have for only a fermion loop of one flavor. We remark that also in the two flavor model one can add the topological term, simply be replacing I |p(n)| (β) with I |p(n)| (2 √ ηη ) ( η/η ) p(n) (see [8]). However, in (13) we use the simpler dual form for θ = 0, since this is the case we actually simulate.
In the dual representation the chemical potentials µ ψ and µ χ couple to the total temporal winding numbers W (l) and W (l) of the loops l and l that correspond to the flavors ψ and χ. In Fig. 2 the double fermion loop on the rhs. of the plot couples to the chemical potential with the factor e (µ ψ +µχ) N T . It is obvious, that the weights in the partition sum (13) are real and positive also for arbitrary values of the chemical potentials µ ψ and µ χ . Thus in the dual formulation the sign problem is gone completely also for the case of finite chemical potential.
In addition to the plaquette expectation number and its susceptibility defined in (8), in the two flavor model with chemical potential we can also study the particle number densities n ψ , n χ and the corresponding susceptibilities χ ψ , χ χ . They are defined as Again we can apply the derivatives also to the dual form of the model to obtain the dual representation of these observables. Their dual form is particularly simple, i.e., the particle densities and susceptibilities are the first and second (connected) moments of the corresponding temporal winding numbers W (l) and W (l). These expressions illustrate an elegant feature of the dual representation: the net particle number can be identified as an integer on every single Note that for the second flavor (dashed lines) which is negatively charged, the gauge flux and the loop must run in the same direction for saturation. This is the reason why also loops from different flavors that run parallel to each other satisfy the gauge constraints. Finally we remark, that the vertical double loop on the rhs. of the plot is an example of a loop that couples to the chemical potential (both flavors have temporal winding number W (l) = W (l) = +1 in that configuration).
configuration. This is not the case for the conventional formulation, where the particle number cannot be defined as an integer on a single configuration. The fact that the particle number can be uniquely identified in the dual representation opens the door to canonical simulations, i.e., simulations at fixed particle number, which we briefly discuss at the end of Section 4 where we consider the two flavor case.

Update strategies and results for the one flavor model
We begin the discussion of the dual update strategies with the one flavor model with a topological term, since the one flavor case is simpler and the two flavor updates build on the steps developed here.

Steps of the update
The partition sum (7) is a sum over configurations of oriented loops and dimers. These have to obey the fermion constraints, i.e., every site of the lattice is either the endpoint of a dimer or is run through by a loop. The weight of each configuration, n I |p(n)| (2 √ ηη ) ( η/η ) p(n) , is computed from the plaquette occupation numbers p(n), which have to be such that the gauge constraints for the links along the fermion loops are obeyed (compare the example in Fig. 1). Let us begin with the update of the dimers in a given fixed configuration of fermion loops. It is easy to see from the example shown in Fig. 1 that different configurations of the dimers are compatible with a given configuration of loops. For example the two vertical dimers inside the loop at the top of the configuration can be replaced by a pair of horizontal dimers. Another possible modification is to shift the vertical line of dimers on the rhs. of the configuration upwards by one link (note that the lattice closes periodically). More generally one can identify closed contours where dimers alternate with empty links and shift all dimers along that contour by one link. This is a move that modifies only the dimers and leaves all fermion constraints intact without altering the fermion loops. Furthermore these moves of the dimers do not change the weights, since all plaquette occupation numbers p(n) remain the same and no Metropolis decision is needed for accepting a pure dimer change. It has to be pointed out that some dimers are frozen by the fermion loops surrounding them, e.g., the horizontal dimer in the top fermion loop, or the dimer inside the 2 × 3 fermion loop at the bottom.
The idea of identifying closed contours where dimers and empty links alternate naturally leads to a worm update. We start the worm from a randomly chosen lattice point which is not part of a fermion loop. At every site the worm makes a random decision along which link to continue next. These choices are restricted by the condition that the worm is not allowed to hit a fermion loop or itself (except when it reaches its starting point). For this bookkeeping the worm blocks the sites and links it already contains. The worm also may end up at a site where the only possible move would be to retrace the previous step. An example is the right end of the frozen horizontal dimer inside the top fermion loop in Fig. 1. In such a case the worm is deleted and no update is performed. Once the worm reaches its starting point it stops and the worm has identified a closed contour where links with dimers alternate with empty links. The update is completed by deleting all dimers along the contour and filling all previously empty links with a dimer. We remark that the dimer worm also includes the simple case of, e.g., rotating a pair of two neighboring horizontal dimers into a vertical pair.
Most important, however, is the fact that the worm can also update winding contours of alternating dimers and empty links, which is not possible with the local pair rotations alone. Only including the dimer worm makes the algorithm ergodic. The fact that the dimers give rise to configurations with topological properties can already be seen from the problem of filling an empty lattice with dimers: It is known (see, e.g., [12]) that in two dimensions the dimer configurations come in four different topological sectors which are not connected by local transformations. Thus it is necessary to include a global update such as the worm strategy discussed here.
Let us now come to discussing the update of the fermion loops and the plaquette occupation numbers attached to them. These updates can be done with three types of local steps on single plaquettes which we illustrate in Fig. 3. In each of these steps the plaquette occupation number of that plaquette changes as p(n) → p(n) ± 1, which leads to a changed weight in the dual representation (7). Consequently all the three types of steps have to be accepted in a Metropolis step with probability min {ρ ± , 1} where and the sign ± is used according to the change p(n) → p(n) ± 1.
A sweep of our algorithm runs through all plaquettes and attempts a change of the fermion variables on that plaquette together with the corresponding plaquette occupation number p(n). Depending on the situation the algorithm finds on the plaquette it works on, the algorithm suggests one of the types of ). For both, the insertion of the fermion loop, as well as for the removal of the fermion loop we have two possibilities which are offered with equal a-priori probability. The corresponding plaquette occupation number p(n) changes to p(n) ± 1. We remark that in all plots of Fig. 3 the changes of the plaquette occupation number p(n) are always denoted relative to the empty configuration, e.g., 0 → ±1. In the general case, e.g., when the plaquette is inside another loop, the change is p(n) → p(n) ± 1. Example of an update step where we either join two fermion loops or split a single fermion loop into two separate loops (which of these two situations applies depends on the global connectivity properties of the fermion loops). The corresponding plaquette occupation number p(n) changes to p(n) ± 1 in these steps.
moves illustrated in Fig. 3, and it is easy to see that these three types of steps together with the dimer worm described above give rise to an ergodic update. The first type of step ( Fig. 3.a., lhs. plot) is for the situation when the plaquette is occupied by a pair of parallel dimers. In this case the algorithm offers to replace the two dimers by an elementary fermion loop around the plaquette, where both possible orientations are offered with equal a-priori probability. The plaquette occupation number changes from p(n) to p(n) ± 1 and the step is accepted with the Metropolis probability min{ρ ± , 1} with ρ ± given in (16). The inverse step is illustrated in the rhs. plot of Fig. 3.a.: here the algorithm attempts to remove an elementary fermion loop around a single plaquette and to replace it with a pair of dimers, where again the horizontal and the vertical possibility are suggested with equal a-priori probability and accepted with probability min{ρ ± , 1} depending on the orientation of the loop that is replaced by the pair of dimers.
The second type of update step for fermion loops and plaquette occupation numbers is illustrated in the lhs. plot of Fig. 3.b. Here one adds a detour to a loop by removing a dimer that sits on the same plaquette as the segment of fermion loop. The plaquette occupation number is changed, p(n) → p(n)±1 according to the orientation of the new path of the fermion loop of that plaquette (in Fig. 3.b. we show the situation for an example relative to the empty configuration p(n) = 0). The move is again accepted with probability min{ρ ± , 1}. The rhs. plot of Fig. 3.b. shows the inverse process, i.e., a detour around a plaquette is taken out of a loop and the fermion constraint is satisfied by placing a dimer.
Finally, in Fig. 3.c. we illustrate the third possible case the algorithm attempts to update at a plaquette: Two anti-parallel segments of fermion loops on two opposing links of a plaquette. In this case the algorithm offers to delete the two segments and to insert them on the other two links of the plaquette. This leads to a configuration where two loops have been joined, or a loop has been separated into two loops (depending on the connectivity properties of the loop(s) the two original loop segments belong to). Again the plaquette occupation number changes as p(n) → p(n) ± 1, giving rise to a Metropolis acceptance probability of min{ρ ± , 1}.
In all other cases, e.g., at plaquettes with a piece of loop around only one corner plus the endpoint of a dimer on the opposite corner, or a plaquette with two pieces of loop running in the same direction on opposite links of the plaquette, the algorithm does not propose a change.
We stress that the moves for changing the fermion loops together with the plaquette occupation numbers we discussed here cannot introduce a non-vanishing net-winding number of loops. The local moves on the plaquettes can only create pairs of fermion loops that wind (spatially or temporally), where one loop in the pair winds forward and the other one backward (in space or time). Since forward and backward loops appear in pairs, no net winding number is introduced. However, this is exactly what the Gauss law enforces: a non-zero winding would correspond to a system that is not charge neutral. It is easy to see that the dual formulation implements the Gauss law in a simple geometrical way: a single loop that winds cannot be saturated with occupying plaquettes. Having understood that winding loops are not compatible with the constraints, it is trivial to see that the three types of steps shown in Fig. 3 give rise to an ergodic update for the fermion loops and the gauge fields in the one flavor case. In the case of two flavors of opposite charge discussed in the next section, we can have pairs of winding fermion loops and there additional strategies (i.e., worms for double fermion loops) will be needed.
For the simulation of the system with the topological term it is advantageous to add another update: a global shift of all plaquette occupation numbers. In this step we offer the change p(n) → p(n) ± 1 ∀n, where the sign ± is randomly chosen with equal probability. Obviously such a gobal shift of the plaquette occupation numbers leaves all constraints intact, since the additional flux on neighboring plaquettes cancels along the link they share. The step is accepted with Metropolis probability min{ρ g,± , 1}, where We stress that this step is not necessary for ergodicity, since such a global shift can also be generated by a loop growing to the size of the whole lattice, such that it wraps around the periodic boundaries and closes back onto itself. However, this is a rare event, and the explicit global change of the plaquette occupation number speeds up the propagation in configuration space.

Results and checks
Having developed the dual representation and a corresponding new update strategy it is paramount to test these. Here we use two types of tests. For the case of θ = 0, i.e., when no topological term is coupled, we can use a standard simulation, since the case of θ = 0 is free of the complex action problem also in the conventional formulation. In the standard simulation the fermions are integrated out and the corresponding fermion determinant is used as an additional weight factor (together with e −S G [U ] ) for the Monte Carlo sampling of the gauge configurations U . The fermion determinant is real and positive, since the eigenvalues of the staggered Dirac operator come in complex conjugate pairs. For the 2-dimensional case of the Schwinger model the discretized Dirac operator is a relatively small matrix, such that for generating the reference data we could simply compute the fermion determinant with a standard library. The second test we implemented was an exact summation over all possible fermion configurations on a small lattice, i.e., fillings of the lattice with dimers and fermion loops, and a numerical sum over all plaquette occupation numbers compatible with a given configuration of the fermion loops. The details for obtaining this exact result on small volumes are described in the appendix. The exact summation is possible also for θ = 0 and allows us to test the dual simulation also in the case where the complex action problem is present.
Before we come to presenting the tests and the numerical results, let us briefly discuss the parameters used in our simulations. We approach the continuum limit as in [3] by simultaneously sending to infinity the inverse gauge coupling β and the lattice volume N S N T at a fixed ratio For the tests shown here a typical value of the constant would be R = 0.1 and we considered lattice sizes up to N S = N T = 24. The values of θ were chosen in a range between −3π and +3π. Note that the geometrical definition of the topological charge gives an integer charge Q[U ] only in the continuum limit, such that also the 2π-periodicity in θ emerges only in that limit and it is useful to monitor a larger interval of θ-values than just the principal branch of [−π, π] (actually since observables are either even or odd in θ, already the interval [0, π] is sufficient in the continuum limit). The dual simulation of the one flavor model typically uses between 10 6 and 10 7 measurements which are decorrelated by a mix of 20 sweeps of local updates for loops and plaquettes, 10 dimer worms and 10 global updates for the gauge plaquettes. At each point in parameter space we use 5 × 10 4 cycles of these mixed updates for equilibration. The errors we show are statistical errors determined with single elimination jackknife. We begin the discussion of our tests with a comparison of all three approaches possible at θ = 0, i.e., conventional simulation, exact summation and dual simulation. In Fig. 4 we show the results of the three calculations as a function of β for lattice size 4 × 4 (larger lattices are too big for the exact summation). It is obvious that the dual simulation (circles), the conventional simulation (squares) and the results of the exact summation (full curve) agree perfectly, supporting the correctness of the dual representation at θ = 0 and the numerical implementation of the dual simulation.
For non-zero vacuum angle θ we can no longer use the conventional approach, but tested our dual simulation program against the results of the exact summation on a small lattice at finite θ the results of this test are presented in Fig. 5 where we show the plaquette expectation value U p (top plot) and the topological charge density q as a function of θ. The data from the dual simulation (circles) and the results of the exact summation (full curve) agree perfectly, confirming our update strategy for the dual representation also at non-zero θ.
We remark at this point that the fact that the comparison to the exact summation is possible only for small volume is not necessarily a disadvantage: on a small lattice configurations that wind in some way We compare the results from the dual simulation (circles) with the data from a conventional simulation (squares) and the results from an exact summation (full curve) on a 4 × 4 lattice. We find perfect agreement of the three results thus confirming the correctness of the dual representation and the update strategy.
around the periodic boundary conditions play a larger role and if not taken into account correctly by the update algorithm would lead to stronger deviations between the Monte Carlo data and the results from exact summation. In our case these configurations are fermion loops that wind around the compactified directions (in the one flavor case these can only be pairs of forward and backward winding loops), as well as sheets of plaquette occupation numbers that cover the entire torus (simply changing all plaquette occupation numbers from p(n) to p(n)±1 in a admissible configuration creates such a sheet which again gives rise to an admissible configuration). Furthermore, as already mentioned, also dimer configurations are known [12] to come in topologically distinct sectors that necessarily imply a global strategy such as our dimer worm for an ergodic update. The fact that we find perfect agreement on the 4 × 4 lattice indicates that also configurations winding around the periodic boundaries are taken into account correctly by the fermion/gauge algorithm and the worm for the dimers.
Having checked the correctness of the algorithm and the dual representation let us briefly analyze the physics of the vacuum angle that emerges from the formulation of the model with massless staggered fermions and the field theoretical definition of the topological charge. The study of the θ-dependence in scalar QED 2 [3] suggests that also here a 2π-periodicity of the dependence of observables on the vacuum angle θ should emerge only in the continuum limit. This is indeed manifest in Fig. 6 where we show the plaquette expectation value U p (lhs. plots) and the topological charge density q (rhs.) as a function of θ. From top to bottom we approach the continuum limit by increasing β and the volume. It is obvious that both observables start to develop the expected 2π-periodicity, as can, e.g., be seen by observing that the distance between two maxima (or minima) approaches the value 2π.
However, there is an important difference between scalar QED 2 and the massless Schwinger model. In the latter case one can use a chiral rotation and the anomaly of the fermion determinant to show that if one of the fermion masses vanishes (note that for staggered fermions we have to compare to the two-flavor continuum results), the physics of the continuum model is independent of θ. Thus it is interesting to analyze whether, and in which form this independence of θ emerges when approaching the continuum limit with the formulation with massless staggered fermions and the field theoretical definition of the topological charge 1 . When inspecting the results for U p (lhs. plots in Fig. 6) one finds that they resemble a parabola overlaid with oscillations, while the results for q resemble an inclined straight line overlaid with oscillations. We first note that in both cases the non-oscillating part becomes more flat towards the continuum limit: for U p the curvature of the parabola decreases when increasing β and the same holds for the slope of the data for q . Furthermore one clearly sees that also the amplitudes of the oscillations decrease quickly toward the continuum limit: the amplitude of the oscillation of U p (lhs. plots) between the minimum at θ = 0 and the first maximum roughly has the values 6 × 10 −2 , 8x8, β = 6.4 12x12, β = 14.4 Figure 6: Plaquette expectation value (lhs. plot) and topological charge density (rhs.) towards the continuum limit. We show the quantities as a function of θ for three different values β = 1. 6, 6.4 and 14.4 (top to bottom) and increase the volume such that R = β/N S N T remains constant at R = 0.1. We observe the emergence of the 2π-periodicity expected in the continuum limit.
the correct implementation of the anomaly by the lattice formulation with massles staggered fermions [14] (although other observables might be sensitive to the order of continuum-and chiral limit). Finally we address a related question in our exploratory study, namely the possibility of a phase transition at θ = π when approaching the thermodynamical limit at fixed inverse coupling β. For pure U(1) gauge theory in 2 dimensions, as well as for scalar QED 2 a first order phase transition emerges at θ = π, which can, e.g., be seen in a jump in q and a diverging topological susceptibility. These findings were established also with the geometrical definition of the topological charge (see, e.g., [3] and the two plots on the rhs. of Fig. 7): a first order jump emergences in q at θ = π, and the maximum in the topological susceptibility grows proportional to the volume.
In the lhs. plots of Fig. 7 we show our results for the topological charge density q (top plot) and the topological susceptibility χ t (bottom) as a function of θ for the Schwinger model. We work at a fixed β = 10.0 and increase the volume from 8 × 8 up to 24 × 24. It is obvious that both observables saturate as the volume increases and do not develop a phase transition, different from the results for the same observables in scalar QED 2 [3]. Also this result for the absence of a transition supports the correct approach of a θ-independent continuum limit.   Figure 7: Topological charge density q (lhs., top plot) and the corresponding susceptibility χ q (lhs., bottom) as a function of the vacuum angle θ. At fixed β = 10.0 we compare the results for different lattice sizes to study the volume dependence. For comparison the rhs. plots show the same quantities for scalar QED 2 , where indeed a first order behavior emerges at θ = π (plots from [3]).

Update strategies and results for the two flavor model with chemical potential
Having established the validity of our algorithmic approach in the somewhat simpler one flavor case with topological charge, and having discussed some of the observables, we now come to the two flavor case.
Here the focus is on finite density, i.e., finite chemical potential, and thus we now set the vacuum angle to θ = 0 (although its inclusion in the two flavor model is straightforward). The algorithmic challenge in the two flavor case is to efficiently update the winding fermion loops that carry the dependence on the chemical potentials µ ψ and µ χ .

Steps of the update
The update of the two flavor model builds on the steps already developed for the one flavor model, since the set of configurations in the dual two flavor partition sum (13) contains products of admissible fermion configurations of the one flavor partition sum (7) for each of the two flavors. Thus as the first part of our algorithm for both flavors we reuse the sweeps of the single flavor updates as discussed in the previous section. We alternately perform the local loop deformations illustrated in Fig. 3 and the dimer worms for each of the flavors, combined with offering a global shift of all plaquette occupation numbers. These updates access all configurations that have zero net winding number of fermion loops around both, spatial as well as temporal direction for each of the flavors. However, as the example in Fig. 2 shows, in the two flavor case also configurations are admissible where the net winding number around space and time has the same non-zero value for fermion loops of both flavors. In the example of Fig. 2 this is the temporal (= vertical) fermion double loop on the right hand side of the configuration. We will now present a worm algorithm that is capable of inserting and removing such fermion double loops. We stress at this point, that the fermion double loop in Fig. 2 is a very special configuration with non-zero net winding number, but it is easy to see that the local loop deformations of the previous section generate all possible instances of configurations in the same winding class. In particular, the two strands of the double loop can become split when, e.g., the update steps shown in Fig. 3.b. are applied to one of the two flavors.
The main challenge of inserting a winding loop is to find a path where the loop can be placed such that the fermion constraints are obeyed. To achieve such an insertion we first run a worm that identifies a closed contour where dimers of both colors alternate 2 . Once such a loop is identified, we replace it by a fermion double loop, where both possible orientations are chosen with equal probability. If the fermion double loop has zero temporal winding the change is always accepted. Otherwise we accept it with probability min {1, e W (µ ψ +µχ)N T }, where W is the temporal winding number of the double loop. This step is illustrated in the lhs. plot of Fig. 8, together with the inverse step, where a fermion double loop (again identified with a worm) is replaced by a closed contour of dimers (rhs.). Together with the single flavor updates of Section 3 this already constitutes an ergodic algorithm for the two flavor model. However, one may increase the efficiency of the update by also including steps that locally deform double fermion loops (see the steps in Fig. 9). In principle these steps lead to changes that are possible also via two other mechanisms: A complete removal of the double loop and a subsequent reinsertion after a dimer change, or as a sequence of two (or more) single flavor steps with a change of plaquette occupation numbers. The former possibility often (depending on the loop) has a low probability and the second one becomes strongly suppressed at small β. Thus we include steps for the fermion double loops that are similar to the local steps in the single flavor update, but do not require a worm search or the activation of plaquettes in an intermediate step. In Fig. 9.a. we show the step where a double dimer is replaced by a detour of the fermion double loop (lhs.) and the corresponding inverse step (rhs.). Fig. 9.b. shows a step where we join two double loops or split them (depending on the overall connectivity properties of the double fermion loop).

Results and checks
For the actual simulation of the two flavor model we use the following combined sweeps: One combined update sweep consists of one sweep of the local updates for each of the two flavors, as well as worm update sweeps for dimers of both flavors as discussed in Section 3.1. This is combined with one global plaquette shift (compare Eq. (17)), one sweep of local neutral loop updates as described in Fig. 9, as well as one winding loop update sweep as shown in Fig. 8. For each parameter set we equilibrate the system by performing 10 5 to 10 6 of these combined sweeps and for the evaluation of observables we use between 10 6 and 4 × 10 7 configurations separated by 6 to 20 combined sweeps for decorrelation. 2 A second possibility would be to find a loop where we have dimers of both colors on the same links alternating with empty links, but since the single flavor dimer loops which we run switch between these two cases, the second possibility needs not be included explicitly. As for the case of the one flavor model with topological term we begin the discussion of our simulation results with a comparison to the results of an exact summation of the partition sum on a small lattice (see the appendix for details of the summation). In the lhs. plot of Fig. 10 we show the particle number density n as a function of the chemical potential µ ψ = µ χ ≡ µ for a 4×4 lattice. The symbols represent the results of our dual simulation for β = 0, 1, 2 and 4, and the full curves are the corresponding reference data from the exact summation. Obviously, for all values of β the dual simulation results agree very well with the curves from the exact summation, which indicates that the dualization and the dual simulation were implemented correctly. The curves for n saturate at 1 as expected for fermions.
In the rhs. plot of Fig. 10 we show our results for the particle number density n as a function of µ for lattice size 8 × 8, again comparing the results for β = 0, 1, 2 and 4. Here we do not have reference data from an exact summation and we simply connect the data points that represent the results from the dual simulation to guide the eye. Comparison of the two figures shows that, as expected, the curves for n become steeper near the inflection point when increasing the volume. The limit of infinite lattice size corresponds to the limit of infinite spatial volume at zero temperature, where one expects that n behaves like a step function and our numerical results show that trend.
Having established the correctness of the dual simulation and thus the fact that the complex action problem is overcome in principle by the simulation in terms of worldlines/worldsheets, we also need to address efficiency issues of the approach. During the explorative studies presented here we partly encountered very long autocorrelations in our simulations, in particular for β > 2 at medium to large values of the chemical potential (µ > 0.5). The reason is the topological nature of the particle number in the dual representation. Increasing the chemical potential increases the particle number and thus the net winding number of the fermionic loops around the compact time direction. At sufficiently high µ this can lead to loops that wind several times around time and then close around the spatial direction. Note that the net winding number has to be the same for both flavors, and in such a configuration the winding loops can already occupy a sizable fraction of the lattice, and for large µ when n approaches 1 each site of the lattice belongs to a winding loop. Clearly such configurations are stabilzed by topology and despite the fact that we use a worm strategy for updating the double loops we found that it is very hard to break up such high-winding configurations.
This topologically stabilized autocorrelation is clearly visible in the time series of the particle number density which we show in Fig. 11. The particle number density n is plotted as a function of the Monte Carlo time t M C measured in units of 10 combined sweeps. The three plots of Fig. 11 are for µ = 0.5, 1.0 and 1.3 (top to bottom) and in each plot we show the time series for several different runs. The long correlations are obvious, in particular for the smaller two values of µ. The plots illustrate that although the complex action problem is solved by mapping the system to a representation with only real and positive contributions, the fermionic nature of the system together with the topological representation of the particle number in terms of the winding number make the Monte Carlo simulation of the system very stiff. In dimensions D ≥ 3 this problem is expected to be much milder, since then the 1-dimensional In particular for small µ we observe oscillations with low frequency corresponding to long autocorrelation times. These are due to topologically stabilized double fermion loops (see the discussion in the text).
winding loops do not separate the lattice into domains that cannot be intersected with another loop. We close this section with remarking that the dual formulation also allows for simulations of the canonical ensemble. Having seen that the system is highly constrained and that changing the particle number is very hard in some parameter regions, canonical simulations, i.e., simulations at fixed particle number, which in the dual representation corresponds to fixed temporal winding number, might be a powerful alternative approach. In order to implement a fixed winding number W one starts the simulation with an initial configuration that has this winding number, e.g., by placing W double fermion loops that wind once. In the subsequent update steps one implements a rejection step for the double fermion line worm whenever it tries to cross the last time-slice of the lattice, such that the worm cannot change the total winding number. Together with the local steps this updates all configurations with fixed winding number W and thus gives rise to a simulation of the canonical ensemble with particle number W for both flavors. First tests in a bosonic model [15] show that in some situations a dual canonical simulation clearly outperforms a grand canonical simulation and canonical simulations in the dual formulation might solve the problem with the long autocorrelations.

Summary and concluding remarks
In this paper we present the first dual simulation of the massless Schwinger model at finite vacuum angle θ and non-zero chemical potential. In the dual representation the partition sum has only real and positive contributions, such that the complex action problem is solved. However, in terms of the new variables, the system is highly constrained. In particular the new fermionic variables, the dimers and fermion loops have to obey the Pauli principle, which in the dual formulation requires that each site of the lattice is either run through by a loop or is the endpoint of a dimer. These constraints require new update strategies and in this paper we explore and evaluate different steps that combine into a suitable algorithm. We remark that in the massive case one also has monomers that can be used for filling the lattice, and since they saturate the fermionic constraints on a single site alone the constraints for the fermions are less rigid in the massive case (however, negative sign contributions re-appear [8]).
More specifically we here study the one flavor model with a topological term (i.e., θ > 0) at zero density and the two flavor model with finite chemical potential at θ = 0. The update for the one flavor model combines local loop deformations, a worm strategy for the update of the dimers and a global shift of the plaquette occupation numbers. The results of the dual simulation are cross checked against a conventional simulation at θ = 0, and against an exact summation on small lattices. We verify that the algorithm is ergodic and reproduces the reference data with high precision, thus establishing the correctness of the approach. In a first small physics application of the new method we analyze how the chosen formulation implements the θ-independence of physics in the continuum limit and analyze the difference to scalar QED 2 for the behavior at θ = π/2.
Generalizing the update strategy of the one flavor model, we construct an ergodic algorithm for the two flavor model, by adding a worm update for the winding double loops that couple to the chemical potential, as well as additional local steps for deforming double loops. Again we compare the results from the dual simulation with the outcome of an exact summation on small lattices and confirm ergodicity and correctness of the dual approach in the two flavor case. When exploring the parameter space at finite µ we found that for some parameter sets the system becomes quite stiff due to topological stabilization: increasing the chemical potential enforces high winding numbers for the fermion loops such that a large fraction of lattice sites is run through by loops that wind around the compact time (and space) directions. Obviously this gives rise to configurations that are hard to change -even with worm strategies. Improving the sampling of loop configurations in that parameter regime is a challenging task that goes beyond the exploratory first study presented here. A possible future strategy for overcoming the stiffness problem of the grand canonical approach at large µ are canonical worldline techniques, which we briefly discuss. configurations with a non-zero net winding number are excluded, since they cannot be saturated with occupied plaquettes (compare the discussion of Gauss' law in Section 2).
In most cases a given loop configuration is compatible with several dimer configurations. In the example in the top row of Fig. 12 there are 4 ways for placing the dimers, in the bottom example there are three possibilities. Since the dimers do not affect the weigths in the dual representation (7) the number of dimer configurations compatible with a given set of loops enters as a degeneracy factor. This factor is 4 in the example in the top row of Fig. 12 and 3 in the bottom.
The next step is to analyze the plaquette occupation numbers that are compatible with the given loop configurations. For this the possible orientations of the loops have to be taken into account. To illustrate the procedure we start with the top row example in Fig. 12 and assume that the loop is oriented in the mathematically negative sense. A possible configuration of plaquette occupation numbers then is given be p(n) = 0 for the 13 plaquettes outside the loop and p(n) = +1 for the 3 plaquettes inside the loop (lattice size is 4 × 4 = 16). However, more generally one can set p(n) = k for the 13 plaquettes outside the loop and p(n) = k +1 for the plaquettes inside, where k is an arbitrary integer. The complete contribution then is obtained by summing over all possible k ∈ Z. In case the orientation of the loop is chosen in the mathematically positive sense (both orientations have to be summed), the admissible plaquette configurations are p(n) = k for the plaquettes outside the loop and p(n) = k − 1 inside. Thus the full contribution of the diagram in the top row of Fig. 12 to the partition sum (7) is given by The asymptotic behavior I k (x) = I −k (x) ∼ (e x/2k) k / √ 2πk of the modified Bessel functions for k → ∞ (see, e.g., [16]) guarantees the (fast) convergence of the sum over k.
For the example in the bottom plot of Fig. 12 we have to sum over all four possible combinations of orientations of the two loops. An easy generalization of the first example then gives the corresponding Similar to the two examples discussed here, one can generate the sums over all possible loop configurations in a simple computer program, together with the corresponding degeneracy factors from the number of dimer configurations compatible with a given loop configuration.
The exact summation program can be generalized to the two flavor case in a straightforward way. The fermion constraints have to hold for both flavors independently, such that one obtains the set of admissible two flavor fermion configurations as the product of the sets of admissible fermion configurations for both flavors. The corresponding plaquette occupation numbers can be determined similarly to the one flavor case, with the small modification that for saturating the gauge constraints for the second flavor the plaquette flux or the orientation of a loop flux from the first flavor have to run in the same direction of the flux of the second flavor since it has opposite charge. For that reason, in the two flavor case also configurations where loops of the two flavors wind the same number of times around the periodic boundaries are possible. In case of a temporal winding they contribute to the dependence on the chemical potentials, i.e., their contributions are weighted with e N T (µ ψ +µχ) W , where W is the temporal winding number (W must be the same for both flavors).