Dual simulation of the 2d U(1) gauge Higgs model at topological angle $\theta = \pi\,$: Critical endpoint behavior

We simulate the 2d U(1) gauge Higgs model on the lattice with a topological angle $\theta$. The corresponding complex action problem is overcome by using a dual representation based on the Villain action appropriately endowed with a $\theta$-term. The Villain action is interpreted as a non-compact gauge theory whose center symmetry is gauged and has the advantage that the topological term is correctly quantized so that $2\pi$ periodicity in $\theta$ is intact. Because of this the $\theta = \pi$ theory has an exact $Z_2$ charge-conjugation symmetry $C$, which is spontaneously broken when the mass-squared of the scalars is large and positive. Lowering the mass squared the symmetry becomes restored in a second order phase transition. Simulating the system at $\theta = \pi$ in its dual form we determine the corresponding critical endpoint as a function of the mass parameter. Using a finite size scaling analysis we determine the critical exponents and show that the transition is in the 2d Ising universality class, as expected.


Introduction
Perhaps one of the most interesting aspects of quantum field theories is that they allow for topological θ-terms -terms which couple to the global properties of the theory, rather than local degrees of freedom. As such they are invisible to perturbation theory, which is local in character. Still, the θ-terms can significantly change the behavior of the theory, in particular since such terms typically violate some symmetries, so their inclusion has huge phenomenological consequences.
On the other hand the effect of such terms is difficult to study. Detailed analytic studies of their effects either require higher symmetry (e.g., supersymmetry or integrability) or a regime which is weakly coupled and amenable to computations. On the other hand such terms are imaginary even in the Euclidean formulation, such that first-principle numerical simulations are hindered by a complex-action problem.
The system we analyze in this paper is the U(1) gauge Higgs model in two dimensions with a topological term. Such models are tightly connected with various interesting systems. In condensed matter physics the U(1) gauge Higgs model with two scalar fields is an effective theory of SU (2) anti-ferromagnetic spin chains, falling into two separate universality classes: integer spin with θ = 0 and half-integer spin for θ = π. Moreover, in the large spin limit the antiferromagnetic spin chain becomes a O(3) nonlinear sigma model 1 with θ = 2πS, where S is the spin quantum number of the underlying system. Finally, such models, being asymptotically free and having classical instanton solutions, are popular also in the particle physics community as toy models of QCD which shares some similar properties. The work we present here is a first step to studying such systems from first principles.
Two values of the topological angle, θ = 0 and θ = π, are special because in those cases the system is invariant under charge conjugation C, which is a global Z 2 symmetry. In the spin systems this symmetry corresponds to a translational symmetry mod 2. Such a symmetry may break spontaneously or may be preserved in the ground-state. For θ = π a scenario that breaks C spontaneously is ensured when the scalar fields are sufficiently massive, in which case the theory is a free U (1) gauge theory and is exactly solvable 2 . On the other hand we expect that as the mass-squared of the scalars decreases towards negative infinity the system will undergo a transition into a trivially gapped phase [4] (i.e., a phase with a gap, and no ground state degeneracy or topological order), which is conjectured to be a 2d Ising transition. This case is in distinction to the case with multiple scalar flavors in which case there is a 't Hooft anomaly between a global flavor symmetry and/or spacetime symmetries and charge conjugation C which forbids a trivially gapped phase. See [1][2][3][4][5][6][7][8] for related discussions.
In the presence of a topological term Monte Carlo simulations are usually plagued by the complex action problem. A possible way to solve the complex action problem is to exactly map the lattice model to a dual form in terms of worldlines and worldsheets, which is the approach we follow here (compare [9][10][11][12][13][14][15][16][17][18][19][20] for use of these techniques in abelian gauge Higgs systems related to the model studied here). More specifically we use the dual formulation that is based on the Villain action [21] appropriately endowed with a θ-term (for an earlier study of the same system with the Wilson action see [18]).
We will arrive at the Villain action from considering a non-compact U (1) lattice gauge theory coupled to matter. Such a theory has a non-compact Z center symmetry, and does not allow for a θ-term (total fluxes on the closed manifold must be identically zero). By gauging the center symmetry, we must introduce a plaquette (a 2-form) Z-valued gauge field, which allows the global Z symmetry to be promoted to a gauge symmetry. Furthermore this allows for a natural introduction of a properly quantized θ-term. We show that by choosing an appropriate gauge the action reduces to the Villain action which we use. However, note that there are other gauges which may be of interest in simulating fermionic theories whose dual descriptions are still plagued by a sign-problem.
Using the Villain action has the advantage that the dual form (a short derivation of the dual form for the Villain action is given in Section 2.2) implements the charge conjugation symmetry that emerges at θ = π as an exact global Z 2 symmetry for the discrete dual variables -we discuss this dual form of the symmetry in detail in Section 2.3.
We study the topological charge and the topological susceptibility numerically as a function of the mass parameter of the Higgs field, keeping the quartic coupling fixed. We find that the system undergoes the expected second order transition at a critical value of the mass parameter with the topological charge being the corresponding order parameter. We use finite size scaling techniques to show that the critical exponents are those of the 2d Ising model. 2 The model and its dual representation 2.1 Continuum description and phase diagram A detailed discussion of the continuum expectations for the phase diagram can be found in [4], and here we only highlight some important points for clarity and completeness.
The problem can be formulated in the continuum with the action where D µ = ∂ µ +iA µ is the U (1) covariant derivative, and F 12 = ∂ 1 A 2 −∂ 2 A 1 is the field strength tensor. The theory is considered on a two-torus T 2 . In that case the flux of F 12 is quantized in integer units of 2π. For that reason the partition function, is periodic with respect to the shifts θ → θ + 2π. We can analyze the system semi-classically in the regime where |m 2 | e 2 . If m 2 is positive and large we can integrate out the φ-field and end up with a pure gauge theory in 1+1 dimensions, which can be exactly solved. The spectrum is given by where L is the spatial volume. To one-loop the effective coupling e 2 ef f is related to the bare coupling e 2 as which is valid as long as m 2 e 2 . Notice that the spectrum is invariant under the change θ → θ + 2π, k → k + 1. Particularly interesting is the θ = π case, where the ground state is degenerate with equal energies for k = 0 and k = 1, i.e., E 0 = E 1 . The topological susceptibility χ t for both of these ground states is entirely determined by the coupling e 2 ef f , On the other hand we could also take m 2 −e 2 . In this case the system is in the Higgs phase, with the modulus of the φ-field weakly fluctuating. Let us set φ = e iϕ v. We then have a term v 2 (∂ µ ϕ + A µ ) 2 (6) in the effective action. This term is minimized by A µ = −∂ µ ϕ, i.e., a pure gauge. The vacuum is unique and charge conjugation invariant (C invariant). Furthermore, the θ-dependence is determined by instanton configurations whose size is cut off by the mass of the φ fluctuations which is proportional to |m 2 |, allowing for instanton calculations. This yields the θ-dependence of the ground state energy given by where S 0 is the instanton action and K is a mass-dimension 2 pre-factor. As we change the parameter m 2 from large and positive to large and negative at θ = π, we must encounter a Z 2 2d Ising phase transition. The phase diagram as a function of m 2 and θ is expected to be as in Fig. 1. The bulk of the current work is to confirm this picture by first-principle simulations.
2.2 Lattice representation: Gauging the center symmetry and the Villain action As we have mentioned several times already, our dualization starts from the Villain action. In order to be useful for our purposes we first have to properly define the θ-term for the Villain action. Further, since we are specifically interested in the physics at θ = π, and in particular the realization of the charge conjugation symmetry at this point, it is vital that the topological charge is properly quantized. To see how to do this, we will start from a slightly unconventional starting point and consider a non-compact gaugetheory, which in the presence of matter has a Z center symmetry. This symmetry is somewhat special, as it acts on line operators (the Polyakov loops) and is dubbed a 1-form symmetry (as distinct from the usual 0form symmetry). This terminology follows [22], and is rooted in the fact that 1-form symmetries act on line operators (i.e., integrals of 1-forms) as opposed to point operators, which can be thought of as 0-forms. To turn the theory into a compact gauge theory this symmetry has to be gauged with a plaquette (or 2-form) Z-valued gauge field B x . We will see that by doing this the theory has a natural θ-term -a sum over the plaquette based fields B x . This theory is endowed with a 1-form discrete gauge symmetry in addition to the usual 0-form gauge symmetry. By fixing a gauge, the formulation becomes equivalent to the Villain action, endowed with a θ-term.
In the conventional lattice representation the dynamical degrees of freedom for the 2-dimensional gauge-Higgs model are the complex valued field variables φ x ∈ C assigned to the sites x and the U(1) valued link variables U x,µ ∈ U(1) assigned to the links (x, µ) of the lattice. Both variables obey periodic boundary conditions on the lattice of size N S ×N T which we denote by Λ. The corresponding partition sum is then given by Here However, here we will take an approach different from the conventional lattice discretization, which, as we shall soon see, will result in the Villain action. To start let us consider a non-compact gauge theory on a lattice. Such a theory has R-valued link variables A x,µ , assigned to the links (x, µ).
We define the discretized field strength F x as In the second step we have, for further use, also defined a lattice derivative d, acting on link fields. The F x are assigned to the plaquettes of the lattice which we label by their lower left corner x. Now we define the partition function as where we have discretized the gauge field action and defined the measure for the non-compact gauge fields as the product measure over all links. Z H is the same Higgs field partition sum as before, but in the argument we made explicit the representation of the gauge links as U x,µ = e iAx,µ . Note that as it stands the partition function is divergent and hence ill defined. The divergence comes about because of a global non-compact center symmetry, as well as the non-compact gauge group. The center symmetry can be thought of as a symmetry which shifts the link fields, A x,µ → A x,µ + λ x,µ , where λ x,µ ∈ 2π × Z, with the constraint that (dλ) x = 0 ∀x. Such transformations are global Z-valued 1-form symmetry transformations as they shift Polyakov loops by a constant, (x,µ)∈l A x,µ → (x,µ)∈l A x,µ + (x,µ)∈l λ x,µ , where l denotes a closed non-contractible loop on our lattice. Notice that the shift of the loop is an integer multiple of 2π independent of its shape. In the absence of matter the link shifts λ x,µ can also be chosen as sets of real numbers with (dλ) x = 0 ∀x.
The theory so far is a non-compact lattice gauge theory, and as such it does not have a θ-term. It is also ill-defined because of the infinity associated with the integration measure.
To address these two issues, we will now gauge the 1-form center symmetry 3 . This promotes the set of λ x,µ to be an arbitrary set of integer multiples of 2π, without the constraints (dλ) x = 0. Clearly the field strength is not gauge invariant because now F x transforms as F x → F x + (dλ) x . To make it gauge invariant we must introduce a plaquette (2-form) gauge field, i.e., a 2πZ-valued field B x on the plaquettes x, and replace F x is a gauge symmetry that makes F x + B x invariant. Note that this gauge symmetry is in addition to the usual gauge symmetry variation which sends A x,µ → A x,µ +φ x+μ −φ x . One should not be surprised that we managed to formulate a theory which has more gauge symmetries than the usual formulation of the compact-gauge theories, as gauge symmetries are redundancies, and are not a measurable feature of the theory. Indeed one could always just fix the gauge, eliminating the redundancy altogether, and the physics should remain unchanged.
Notice that now also Wilson loops in the form (x,µ)∈l A x,µ , where the loop l closes trivially on the torus are not gauge invariant under the 1-form gauge symmetry described above. This means they are not observables in our theory. Instead a Wilson loop l defined as e ik (x,µ)∈l Ax,µ with k ∈ Z is indeed invariant and thus an observable. This indicates that by gauging the Z center symmetry of the non-compact gauge theory, we have obtained a compact gauge theory, with the gauge action given by We note that the theory can be endowed with a natural θ-term given by Note that this is a natural θ-term since it is related to d 2 xF x in the continuum which was replaced by x (F x + B x ) in our discretization. On the lattice F x is a total derivative, so it drops out when summed over x, such that x (F x + B x ) = x B x , which is the form of S θ on the rhs. of (12). Clearly e −S θ is gauge invariant and moreover, since B x ∈ 2πZ, the partition function with a topological term based on (12) is periodic under θ → θ + 2π. Note that at the special values θ = 0 and θ = π mod 2π the system also has a charge conjugation symmetry C, implemented as φ This charge conjugation symmetry and our correct implementation on the lattice will play a key role for the critical endpoint behavior discussed in this paper.
Finally we address the question of ill-definitnes due to an infinity which stems from the integration range (−∞, ∞) of the non-compact link fields A x,µ . Because we can always use a gauge symmetry where λ x,µ = 2π m x,µ with m x,µ ∈ Z, we can use the gauge invariance to fix the gauge such that A x,µ ∈ (−π, π]. To impose this gauge let us introduce a function Θ −π,π (x) which equals 1 for −π < x ≤ π and vanishes otherwise. By inserting unity in the form the partition function turns into where {k} = x,µ kx,µ∈Z and {B} = x Bx∈2πZ . We can now per- (15) Now we note that the infinite sum {k} is just a multiplicative factor and can be dropped. After dropping this term, and setting B x = 2πn x , where n x ∈ Z, we obtain the Villain representation [21] of the partition function, now containing the correctly discretized 2π periodic θ-term that implements charge conjugation as an exact Z 2 symmetry.
It is important to stress that the Villain action is just a particular gauge in our discretization framework. We could have made other choices. For example we could make a choice which eliminates the field B x from all plaquettes except one, labelled by x 0 . What makes this gauge attractive is that we could then use it to propose "instanton" updates simply by proposing a change in B x 0 . This eliminates the topological charge sampling problem, because it reduces it to just one variable: the value of B x 0 . While here we will not make use of this gauge, it could potentially be of use when simulating theories with fermions, in which case dual descriptions are still plagued with the sign-problem.
We can now summarize the lattice discretization by providing the final form of the Boltzmann factor and specifying also the Higgs field action and Higgs field partition sum. The gauge field Boltzmann factor in the Villain form with the topological term is given by where in the second expression for B G [U ] we have inserted F x in the exponent coupling to θ, which, as we have already discussed, is an equivalent form because of x F x = 0. The latter expression will be used in the appendix for transforming B G [U ] to its dual form. The partition function Z H [U ] for the Higgs field in a background U = {e iAx,µ } reads where the action S H [φ, U ] for the Higgs field is given by λ denotes the quartic coupling and the mass parameter M is defined as The path integral measures are a product of U(1) Haar measures for the gauge fields and a product of integrals over C for the Higgs field, The charge conjugation transformation which we have addressed corre- It is easy to see that this transformation is a symmetry of the action S H [φ, U ] for the Higgs field (18) and of the path integral measure. For the field strength the transformation gives rise to F x → −F x such that n x → −n x compensates this in the quadratic part of the exponent of the gauge field Boltzmann factor (16), while the linear part of the exponent is symmetric only at θ = π (and of course at the trivial point θ = 0). Analyzing the spontaneous breaking of this symmetry for the θ = π case with Monte Carlo simulation techniques is the goal of this paper.
It is obvious that a non-zero vacuum angle in the representation (16) gives rise to a complex Boltzmann factor and thus to a so-called complex action problem. This complex action problem can be solved completely by mapping the partition sum exactly to a dual form where the new degrees of freedom are worldlines for matter fields and world-sheets for the gauge fields.

Dual representation
This section summarizes the exact transformation of the model in Villain formulation to the corresponding dual formulation in terms of worldlines and plaquette occupation numbers. We begin the presentation of the mapping to the dual form for the Villain action by quoting from the literature the dual representation of the partition function Z H [U ] for the Higgs field in a gauge field background. This form is derived by expanding the Boltzmann factors for the individual nearest neighbor terms and subsequently integrating out the original field variables φ x (see, e.g., the appendix of [11] for the derivation of the form used here). In the dual form the new degrees of freedom are flux variables j x,µ ∈ Z assigned to the links of the lattice and auxiliary variables h x,µ ∈ N 0 also living on the links. The partition sum Z H [U ] is a sum {j,h} ≡ x,µ jx,µ∈Z hx,µ∈N 0 over all possible configurations of these variables. It is given by Each configuration of the flux variables j x,µ and the auxiliary variables h x,µ comes with a weight factor W H [j, h] given by The first part of the weight factor, i.e., the product over the links collects the combinatorial factor from the expansion of the nearest neighbor Boltzmann factors. This is followed by a product over all sites where the integrals I (f x ) emerge when integrating out the moduli of the complex fields φ x . These integrals depend on the non-negative combinations f x of the flux and auxiliary variables given in (23). For the Monte Carlo simulation the integrals I (f x ) are numerically pre-calculated and stored for a sufficient number of values f x for the chosen parameters M and λ.
Integrating over the complex phases of the fields φ x generates constraints for the flux variables j x,µ . These constraints appear as a product of Kronecker deltas in (20) where we use the notation δ(n) = δ n,0 . The constraints enforce vanishing divergence of the flux j x,µ at every site x, i.e., Along the closed loops of the j-flux the U(1) variables on the links (x, µ) of the loop are activated according to the flux j x,µ . This gives rise to the last product in (20) where for negative values of j , for links that are run through in negative direction the complex conjugate link variable is taken into account. The link variables along the loops now have to be integrated over with the gauge field Boltzmann factor B G [U ]. As we show in the appendix, one may use Fourier transformation with respect to the field strength F x to write the Boltzmann factor of the Villain action in the form At every plaquette of the lattice we sum over a plaquette occupation number p x ∈ Z, and the corresponding sum over all configurations will be denoted as {p} ≡ x px∈Z . The configurations of the plaquette occupation numbers come with a Gaussian weight factor at every site that depends on both, the inverse gauge coupling β and the vacuum angle θ. For this weight factor we introduce the abbreviation (V = N S N T denotes the number of sites) The gauge field dependence of the Boltzmann factor B G [U ] is collected in the last term of (24), which is a product over the Fourier factors e iFxpx at all plaquettes. Inserting the explicit expression (9) we can reorganize this product, where in the last step we have shifted the nearest neighbor displacement from the A x,µ to the plaquette occupation numbers p x by accordingly relabelling the variable x in the product over all sites. In this form we can combine this factor with the gauge field dependence of Z H [U ] in (20) and integrate over the link variables. We find where in the last step we have used π −π dA/2π e iA n = δ(n). Thus the integration over the U(1) link variables gives rise to a second set of constraints which connect the j-flux and the plaquette occupation numbers p x such that at each link the sum of j-flux and flux from the plaquettes containing that link has to vanish.
Putting things together we find the final form of the dual representation The partition function is a sum over all configurations of the link-based flux variables j x,µ ∈ Z, the link-based auxiliary variables h x,µ ∈ N 0 and the plaquette occupation numbers p x ∈ Z. The j x,µ must obey the zero divergence constraint that forces admissible configurations to have the form of closed loops of j-flux. A second set of constraints links these fluxes to the plaquette occupation numbers such that the combination of j-flux and flux from the plaquettes vanishes at each link of the lattice. As a consequence admissible configurations of the plaquette variables are either closed surfaces formed by occupied plaquettes or open surfaces with j-flux at the boundaries.
The configurations come with a weight factor W H [j, h] for the j-flux and the unconstrained auxiliary variables h x,µ which together describe the Higgs field dynamics in the dual language. This weight factor is given explicitly in (21). The second weight factor W G [p] for the plaquette occupation numbers governs the gauge dynamics and is given in (25). The weight has the form of a Gaussian for the plaquette occupation number p x at each site. The width is given by the square root of the inverse gauge coupling and the mean of the Gaussian is at θ/2π. Obviously all weights are positive at arbitrary values of θ and thus the dual representation solves the complex action problem.

Charge conjugation symmetry at θ = π and observables
Let us now discuss an important aspect of the dual representation for the Villain action, i.e., the implementation of the charge conjugation symmetry at θ = π as global Z 2 symmetry for the discrete dual variables. At θ = π the product of Gaussians in the weight leaves the weight factor W G [p] invariant. In order to show that this gives rise to a full symmetry we need to discuss the interaction with the flux variables j x,µ , encoded in the constraints δ j x,1 + p x − p x−2 δ j x,2 − p x + p x−1 appearing in (28). First we note that a uniform constant shift of all p x , such as the additive term −1 on the rhs. of (29) leaves the constraints intact, since in each of the Kronecker deltas the difference of two plaquette occupation numbers appears, such that a constant shift cancels. Flipping the sign of the p x , which is the second ingredient of the transformation (29), can be compensated in the constraints by changing the signs of all j-fluxes, The final step for completely establishing the symmetry is to note that the flux variables j x,µ enter in the weight factors W H [j, a] in Eq. (21) only via their absolute values |j x,µ |, such that also the weight factor W H [j, a] remains unchanged. Thus the combination of (29) and (30) implements a global Z 2 symmetry for the discrete dual variables. The fact that the symmetry is Z 2 is obvious, since applying the transformation (29), (30) twice gives the identity transformation (p x = p x ∀x, j x,µ = j x,µ ∀x, µ).
In this paper we want to study the Z 2 symmetry and thus need a corresponding order parameter. A suitable order parameter is given by the topological charge density 4 q. Its vacuum expectation value q and the corresponding topological susceptibility χ t are defined as As an additional observable we also study the gauge field action density s G defined as In order to determine the dual form of these observables we perform the derivatives with respect to θ using the dual representation for Z and obtain It is easy to see that for θ = π the topological charge density q, which in the dual representation (33) is given by q | θ=π = (V 2πβ) −1 x (p x + 1/2), changes sign under the Z 2 transformation (29), (30), and thus is a suitable order parameter for the breaking of that symmetry. Note that on a finite lattice the Z 2 symmetry at θ = π cannot be broken spontaneously such that q vanishes exactly. Thus in our numerical analysis at θ = π we will partly evaluate |q| , which also can be used to identify and analyze the phase transition.
The situation which we find for the U(1) gauge-Higgs model at θ = π thus is reminiscent of a ferromagnet. The corresponding order parameter, i.e., the magnetization here is replaced by the topological charge density q. The role of the external magnetic field is played by i.e., the topological angle shifted such that ∆θ = 0 corresponds to the symmetrical point where explicit symmetry breaking is absent. Finally we need to identify the control parameter that drives the system through the phase transition. For the ferromagnet this parameter is the inverse temperature, which here is replaced by the parameters M and λ. In our analysis we will keep λ fixed and use the mass parameter M = 4 + m 2 as the control parameter to drive the system at the symmetrical value ∆θ = 0 (i.e., θ = π) through the expected phase transition.

Numerical simulation 3.1 Details of the simulation
For the Monte Carlo updates of abelian gauge Higgs models in the dual representation different strategies can be followed. One type of strategy offers only trial configurations that already obey the constraints. This strategy can, e.g., be implemented by offering the following local changes where x 0 is some site of the lattice and the increment ∆ is randomly chosen in {−1, +1}. The proposed changes (37) are accepted with the usual Metropolis probability that can easily be computed from the weights of the dual representation and the new configuration generated by (37) obeys all constraints. Together with a conventional Metropolis update of the auxiliary variables h x,µ , which are not subject to constraints, the update (37) already constitutes an ergodic algorithm. However, for faster decorrelation it is useful to also add a global update of all plaquette occupation numbers, where again ∆ is an increment randomly chosen in {−1, +1}. The other strategy for the update of the dual gauge Higgs model is based on the idea of the worm algorithm [23]. In an initial step one locally violates the constraints by changing some of the dual variables. Subsequent update steps randomly propagate the defect on the lattice until the starting point is reached and the defect is healed. All intermediate propagation steps are accepted with suitable Metropolis decisions and after the worm has closed, one again has a configuration that obeys all constraints and the total transition probability fulfills detailed balance. For applying the strategy to dual abelian gauge theories the structures generated in the propagation steps are surfaces of plaquette occupation numbers bounded by j-flux. This socalled surface worm strategy was discussed in detail in [11] and in some parameter regions clearly outperforms the local updates in four dimensions. However, here we simulate the considerably less demanding 2-dimensional case where the updates (37) and (38) are sufficient for very accurate results.
In our numerical simulation we combine sweeps of the local updates (37) and the global plaquette occupation number updates (38). More specifically, for each parameter set we perform 10 5 mixed sweeps for equilibration. The statistics ranges between 10 4 and 10 7 measurements, depending on the parameter region, i.e., the distance to the critical point. We typically use 20 to 50 (for the larger lattices) mixed sweeps in between measurements for decorrelation. The errors we show for our results are statistical errors determined from a jackknife analysis.
We carefully tested our dual simulation algorithms using several test cases. For pure gauge theory at arbitrary θ one can compute the partition sum in closed form: In that case the constraints for the remaining dual variables p x imply that only constant configurations ∀x : p x = j ∈ Z are admissible and the partition function and observables are obtained as fast converging sums over j ∈ Z (in fact Z is given by the θ 3 elliptic function). Our simulation algorithm reproduces this analytic result with very high statistics. A second case where we tested the dual simulation is the full theory at θ = 0. In that case the complex action problem is absent and reference simulations with the conventional representation are possible. Also for this test we found excellent agreement of results from the dual simulation with those from the conventional representation.
In all our numerical simulations we kept the quartic coupling fixed at a value of λ = 0.5. For the inverse gauge coupling β we used β = 3.0, and β = 5.0.

Results for bulk observables
We begin the presentation of our results with a plot of the topological charge density q and the gauge action density s G as a function of the symmetry breaking parameter ∆θ. We study the behavior for two different values of the mass parameter M = 4+m 2 , given by M = 2.0 (lhs. plots) and M = 3.5 (rhs.). The other parameters are β = 3.0 and λ = 0.5 and in Fig. 2 we show our results for three different volumes.
It is obvious that for a mass parameter value of M = 2.0 both observables q and s G show only a very weak dependence on the symmetry breaking parameter ∆θ = θ − π and essentially no scaling with the volume. This clearly is different for M = 3.5, where in the thermodynamic limit we see the emergence of a first order jump in q , which has also been observed in both, pure U(1) gauge theory in 2-dimensions [24], as well as in the simulation of the U(1) gauge Higgs system with the Wilson gauge action [18]. Also the gauge action density is sensitive to ∆θ at M = 3.5 and in the thermodynamic limit develops a cusp at the symmetrical point ∆θ = 0, i.e., at θ = π.
From the fact that at M = 3.5 we observe a θ-dependence of the observables with an emerging first order behavior at θ = π and only very weak dependence on θ for M = 2.0, we expect that the first order line at θ = π must have an endpoint at some critical value M c somewhere between M = 2.0 and M = 3.5.
In order to locate and analyze the phase transition we study |q| and χ t as a function of the control parameter M = 4 + m 2 . In Fig. 3 Figure 3: The topological charge density |q| and the topological susceptibility shifted by a constant χ t + β 4π at θ = π for different volumes. We compare the results for β = 3.0 (top plots) and β = 5.0 (bottom) at λ = 0.5 and plot the observables as a function of our control parameter, i.e., the mass parameter M = 4 + m 2 . For comparison on the top of the plots we also show the horizontal axis labelled with m 2 , which at the point of the transition is negative. lattice (N S = N T = L here) is expected to be χ t,max ∼ L γ/ν = L 7/4 for the conjectured Ising transition. Despite the fact that for a precise determination of the maxima a much finer spacing of the data points would be needed, it is already plausible that the maxima indeed scale as conjectured. In the next section we now determine the critical mass parameter M c as well as some of the critical exponents using precise finite site scaling techniques.

Definition of the critical exponents and the strategy for their determination from finite size scaling
In this subsection we define the critical exponents we want to determine and briefly summarize the setting of our finite scaling analysis. This is of course text book material (see, e.g., [25,26]), but needs to be slightly adapted to the observables of the U(1) gauge Higgs model studied here.
In our system the mass parameter M is the quantity that drives the system through the phase transition at the critical value M c . Consequently we define the reduced mass m r , We already discussed that the topological charge density q has to be identified with the magnetization, and the topological susceptibility χ t with the magnetic susceptibility. Near the transition at m = 0 these observables behave as where ξ is the correlation length. The terminology for the corresponding critical exponents β, γ and ν follows the conventions for ferromagnets. Using the third equation in (40) to express |m r | in terms of ξ and the fact that the correlation length is limited by L (we use lattices with size N S = N T = L), we can scale out the leading term in |q| and χ t to find the finite size scaling form where x = L 1 ν m r is the scaling variable and F q (x) and F χ (x) are the scaling functions that are regular at m r = 0 and have no further L dependence in the scaling region around m r = 0.
For the determination of the critical exponents we implement the Ferrenberg-Landau strategy [27]. We introduce the fourth order Binder cumulant [28] The maximum slope of this cumulant scales with L 1 ν , such that the peak of its mass-derivative scales as [28] dU The first step of our determination of the critical exponents (see the next subsection for the discussion of the actual implementation) is the determination of ν based on (43), with a numerically computed derivative dU/dM of the Binder cumulant based on a multiple-histogram technique (see below). In addition we also determine ν by studying observables that have the same scaling behavior as U , for instance the logarithmic derivatives of moments of the topological charge. In particular we study the derivatives which have maxima that also scale with L 1 ν . Thus these derivatives can again be fit as described in (43) and allow for an independent determination of ν.
Once the exponent ν has been calculated with satisfactory precision we can determine the critical mass M c . We first estimate the pseudo-critical mass values M pc (L) by determining the location of the peaks in observables evaluated at lattice size L. Suitable observables for this analysis are the derivative dU/dM (43), the logarithmic derivatives (44) as well as the location of the maximum slope in |q| and the peak in the susceptibility χ t . The critical mass M c was then determined in a fit of the data for the pseudo-critical mass values M pc (L) with the scaling ansatz Once ν and M c are determined, we know the zero of the reduced mass m defined in (39) and the scaling variable x = L 1 ν m r . Thus we can compute the remaining critical exponents from the finite size scaling relations (41). β is determined from the scaling of |q| (mr,L) , and γ from the scaling of χ t (mr,L) .
To reliable determine the precise locations M pc (L) of the peaks in observables we use the multiple-histogram method [29]. This approach allows one to perform only a few simulations in the region around the peaks and to compute the observables at a continuum of values M with high precision. An error estimate for the necessary interpolation can be obtained by jackknife resampling of the data.

Numerical results for the critical exponents
For the determination of the critical exponents with the finite size scaling strategy outlined in the previous subsection we use L × L lattices with sizes of L = 16, 32, 48, 64, 80, 96 and 112 and an inverse gauge coupling of β = 3.0. As already outlined, for each lattice size L we perform simulations at 4 different values of the mass M within the critical region and interpolate with fine interpolation steps using the multiple histogram method [29]. For an error estimation we divide our measurements into blocks and perform a jackknife resampling of the data blocks. The critical exponents and the critical mass are then extracted following the strategy outlined in the previous subsection. For the simulations to determine the critical exponents we use ensembles with 4 × 10 6 configurations after 5 × 10 4 initial equilibration sweeps and 50 decorrelation sweeps between the measurements.
The first step in the analysis is the determination of the critical exponent ν, which can be performed independent of other critical exponents and without the knowledge of the critical mass M c . As discussed we determine ν from the maxima of the mass-derivatives of U , ln |q| and ln q 2 (compare Eqs. (43) and (44)). In the lhs. plot of Fig. 4 the symbols indicate the maxima of the mass derivatives of the three observables plotted against the lattice size. In our double logarithmic plot the maxima fall onto straight lines that are essentially parallel and have a slope given by 1/ν. From a fit of the maxima with the scaling behavior (43) we determine ν.
By performing a series of fits where we vary the minimal size L taken into account in the fit range we can see that up to L = 32 the Binder cumulant is sensitive to additional finite size effects and we thus omit the corresponding data for the smallest two lattice sizes from the fit, while for the other observables all measurements could be used. The values of ν extracted from the fits of the different observables were combined in an average leading to a final result of ν = 1.003 (11). The statistical error was propagated through the entire analysis chain starting with the initial jackknife resampling of the data.
With the knowledge of the critical exponent ν we can now determine the critical mass M c from the scaling relation (45) for the pseudo-critical values M pc (L) obtained as the maxima of suitable observables. As outlined, we use the mass derivatives of U , |q| , ln |q| and ln q 2 as well as the location of the peaks in χ t . The symbols in the rhs. plot of Fig. 4 show the corresponding values M pc (L) as a function of L −ν , and the straight lines the fits with (45). We indeed observe the expected linear behavior and the fits for the different observables extrapolate for L → ∞ to the corresponding estimates for M c . 16   In both plots we annotate the final results for the parameters determined from the fits.
Again we varied the fit range to analyze higher order finite size corrections and discarded the L = 16 data for all observables and in addition removed the L = 32 data for the Binder cumulant. The extrapolated results for the critical mass were then averaged with a final result of M c = 2.989 (2).
After estimating ν and M c we can finally extract the exponents β and γ from the scaling of |q| (mr,L) and χ t (mr,L) , as given in (41). In the double logarithmic plots in Fig. 5 we show the results for |q| (mr,L) and χ t (mr,L) at M c , i.e., at m r = 0 as a function of L. Fits according to the scaling relations (41) appear as straight lines and their slopes are given by −β/ν and by γ/ν. Using the previously determined result for ν we obtain values of β = 0.126(7) and γ = 1.73 (7).
This concludes our determination of the critical exponents and the crit- ical mass parameter M . We summarize our results in Table 1, where we display also the corresponding critical values of the 2d Ising model. We observe good agreement within error bars, which provides strong numerical evidence that the critical endpoint of the 2d gauge Higgs model at θ = π is in the 2d Ising universality class.

Summary and future prospects
In this paper we have presented the results of our lattice analysis of the critical endpoint for the 2d U(1) gauge Higgs model at topological angle θ = π, where charge conjugation is an exact symmetry. Charge conjugation is spontaneously broken when the scalars with a gauge charge are heavy. As the scalar mass squared is lowered and becomes sufficiently negative, the symmetry can become unbroken with a transition that is conjectured to be in the 2d Ising universality class. For studying this scenario numerically, the lattice formulation has to deal with two main challenges, the implementation of the topological charge as an integer in order to make charge conjugation an exact symmetry at θ = π, and the complex action problem at non-zero topological angle. The former problem is overcome here by using the Villain action, while the latter challenge is solved by using a dual representation in terms of worldlines and plaquette occupation numbers for the simulation.
Using the topological charge density q, and the corresponding topological susceptibility χ t as our main observables we show that for sufficiently large mass parameter M = 4 + m 2 the system undergoes a first order transition at θ = π. As the mass parameter is lowered this first order behavior disappears in a critical endpoint M c . As a function of M the observables q and χ t indicate a second order transition at M c . Using finite size scaling analysis we determine the corresponding critical exponents and show that the transition is in the 2d Ising universality class as conjectured.
The work presented here is only the first step in exploring the rich physics of 2d quantum field theories at finite topological θ-angle. The methods presented here can be adopted for studying theories with more scalar matter fields. If the fields have the same mass, the global symmetry group is enhanced, and there exist nontrivial 't Hooft anomaly matching conditions at θ = π, which guarantee that there is no trivially gapped phase (see [4][5][6][7][8]). Another interesting case is to promote the charge-1 scalar field to a charge-2 scalar field. Such a theory, in addition to the C-symmetry has a Z 2 center symmetry, which leads to a mixed 't Hooft anomaly. This model, which can be studied on the lattice with our methods, will exhibit a very different behavior 5 from the charge-1 model. Similar reasoning can also be used to formulate compact U (1) gauge theories in higher dimensions. By doing this, the theory of compact electrodynamics can be formulated which allows for an exclusion of dynamical monopoles as well as complete control over their contributions to a partition sum (charges, actions, etc.). One application of such a formulation is that a theory with only fixed monopole charges (typically multiples of 2 or 4) can be formulated, or indeed even without any monopoles. Such theories are the effective theories of anti-ferromagnets, and have been connected with deconfined quantum criticality 6 and emergent symmetries at the phase transition between the monopole phase, and the Higgs phase.
Finally our approach can also be adopted for studying various non-linear sigma models, such as SU (N )/U (1) N −1 , which have N − 1 θ-angles and can be formulated in terms of U (1) N −1 scalar gauge theories. For N = 2 this is the famous O(3) or CP (1) nonlinear sigma model, which at θ = π is an effective theory of half-integer spin chains, while for N = 3 the model corresponds to an effective theory of SU (3) spin chains [30][31][32]. Being asymptotically free, the model also has two θ-angles, and a rich phase diagram as a function of the two angles (see [7,32]) which can be explored with the methods presented here 7 .

Appendix
In order to keep this paper self-contained, in this appendix we derive the Fourier representation of the Boltzmann factor B G [U ] which we use in (24). In its conventional form the Boltzmann factor with Villain action is the product over factors