Complex Langevin analysis of 2D U(1) gauge theory on a torus with a $\theta$ term

Monte Carlo simulation of gauge theories with a $\theta$ term is known to be extremely difficult due to the sign problem. Recently there has been major progress in solving this problem based on the idea of complexifying dynamical variables. Here we consider the complex Langevin method (CLM), which is a promising approach for its low computational cost. The drawback of this method, however, is the existence of a condition that has to be met in order for the results to be correct. As a first step, we apply the method to 2D U(1) gauge theory on a torus with a $\theta$ term, which can be solved analytically. We find that a naive implementation of the method fails because of the topological nature of the $\theta$ term. In order to circumvent this problem, we simulate the same theory on a punctured torus, which is equivalent to the original model in the infinite volume limit for $ |\theta|<\pi$. Rather surprisingly, we find that the CLM works and reproduces the exact results for a punctured torus even at large $\theta$, where the link variables near the puncture become very far from being unitary.


Introduction
The θ term provides an interesting avenue of research in quantum field theories. Due to its topological nature, its effects on physics should be genuinely nonperturbative, if present at all. In particular, it does not affect the equation of motion, which implies that θ is a parameter that does not exist in the corresponding classical theory. For instance, the θ term in QCD is given by S θ = −iθQ, where Q is the topological charge defined by which takes integer values on a compact space. This term is renormalizable by powercounting and hence it is a term which is perfectly sensible to add in the action. However, it breaks parity and time-reversal symmetries, and hence the CP symmetry. This leads to a non-vanishing electric dipole moment of a neutron, which is severely restricted by experiments. The upper bound on θ thus obtained is |θ| 10 −10 [1], which is extremely small although there is no reason for it theoretically. This is a naturalness problem known as the strong CP problem.
A popular solution to this problem is the Peccei-Quinn mechanism [2,3,4,5], which introduces axions as a pseudo Nambu-Goldstone boson of a hypothetical global U(1) PQ symmetry. In this mechanism, the potential for the axions induced by QCD chooses a CP invariant vacuum automatically. Recently, gauge theories with a θ term have attracted attention also from the viewpoint of the 't Hooft anomaly matching condition [6,7,8] and the gauge-gravity correspondence [9,10,11,12]. In particular, there is an interesting prediction for a phase transition at θ = π, which claims that either spontaneous CP breaking or deconfinement should occur there [6,7]. In order to investigate gauge theories from first principles in the presence of a θ term motivated either by the physics related to axions or by the recent predictions, one needs to perform nonperturbative calculations based on Monte Carlo methods. However, this is known to be extremely difficult because the θ term appears as a purely imaginary term in the Euclidean action S. The Boltzmann weight e −S becomes complex, and one cannot interpret it as the probability distribution as one does in Monte Carlo methods.
One can still use the reweighting method by treating the phase of the complex weight as a part of the observable. In the case at hand, this amounts to obtaining the histogram of the topological charge at θ = 0 and taking an average over the topological sectors characterized by the integer Q with the weight e iθQ . Various results obtained in this way are nicely reviewed in Ref. [13]. Clearly, the calculation becomes extremely difficult due to huge cancellations between topological sectors when topological sectors with |Q| ≫ π/|θ| make significant contributions to the partition function, which occurs either for |θ| ∼ π or for smaller |θ| with sufficiently large volume. This is the so-called sign problem which occurs in general for systems with a complex action.
There are actually a few promising methods that can handle systems which suffer from this problem such as the complex Langevin method (CLM) [14,15,16,17,18,19], the generalized Lefschetz thimble method [20,21,22,23,24,25], the path optimization method [26,27,28,29] and the tensor network method [30,31,32,33,34]. Each method has its pros and cons, however. In this work, we focus on the CLM, which can be applied to various physically interesting models with large system size in a straightforward manner. (See, for instance, Refs. [35,36,37].) The only drawback of the method is that it can give wrong results depending on the system, the parameter region, and even on the choice of the dynamical variables. Recently, the reason for this behavior has been understood theoretically [16,17,18,19,38,39]. In particular, Ref. [19] proposed a practical criterion for correct convergence, which made the CLM a method of choice for complex action systems as long as the criterion is met. There are indeed many successful applications to lattice quantum field theory [40,41,42,43,44,45,46,47,48,49,50,51,35] and matrix models [52,53,54,55,36,47,56,57] with complex actions.
As a first step towards its application to 4D non-Abelian gauge theories, here we apply it to 2D U(1) lattice gauge theory with a θ term, which suffers from a severe sign problem despite its simplicity. 1 In fact, the model can be solved analytically with finite lattice spacing and finite volume on an arbitrary manifold [60,61,62], which makes it a useful testing ground for new methods [60,63,64,65] aiming at solving the sign problem. By using the reweighting method [64], for instance, one can only reach θ ∼ 2.2 with a 16 × 16 lattice, and in particular, it seems almost impossible to approach θ = π by this method. Note also that the region of θ that can be explored by this method shrinks to zero as one increases the lattice size.
We find that a naive implementation of the CLM fails. The reason for this is that the configurations that appear when the topology change occurs during the Langevin process necessarily result in a large drift term. Due to this fact, the criterion [19] for the validity of the method based on the histogram of the drift term cannot be satisfied. If one tries to suppress the appearance of the problematic configurations by approaching the continuum limit, the criterion can be satisfied, but the topology change does not occur during the Langevin process, hence the ergodicity is lost.
In order to cure this problem, we introduce a puncture on the torus, which makes the base manifold noncompact. We have obtained exact results for this punctured model as well. Even in the continuum limit, the topological charge is no longer restricted to integer values and the 2π periodicity in θ does not hold. However, if we take the infinite volume limit with |θ| < π, one cannot distinguish the model from the original non-punctured model as far as the observables that make sense in that limit are concerned. Note that in that limit, the topological charge can take arbitrarily large values and therefore it does not really matter whether it is an integer or not.
On the other hand, the situation of the complex Langevin simulation changes drastically for the punctured model. The topology change occurs freely and the appearance of the problematic configurations can be suppressed by simply approaching the continuum limit. Thus the criterion for the validity of the CLM is met without losing the ergodicity, and we are able to reproduce the exact results for the punctured model.
The most striking aspect of our results is that the CLM works even if the link variables close to the puncture become very far from being unitary. This can happen because the direct effect of the θ term on the complex Langevin dynamics is actually concentrated on these link variables. While the link variables are allowed to be non-unitary in the CLM in general in order to include the effects of the complex action, all the previous work suggested that the condition for the validity cannot be satisfied unless the non-unitarity is sufficiently suppressed. Precisely for this reason, the gauge cooling [66,18,19] was invented as a crucial technique in applying the CLM to gauge theories. In fact, we also use the gauge cooling in our simulation, but the link variables close to the puncture nevertheless become far from being unitary when θ or the physical volume gets large. Yet the criterion for the drift term is not violated and the exact results are perfectly reproduced.
The rest of this paper is organized as follows. In Section 2, we make a brief review on 2D U(1) gauge theory on a torus with a θ term and discuss how to put it on a lattice. In Section 3, we apply the CLM to this theory and show that a naive implementation fails. In Section 4, we introduce a puncture on the torus in order to circumvent the problem encountered in Section 3, and show the equivalence of the punctured model and the original non-punctured model in the infinite volume limit for |θ| < π. In Section 5, we apply the CLM to the punctured model and show that it works perfectly even at large θ, where the link variables near the puncture become very far from being unitary. Section 6 is devoted to a summary and discussions. In Appendix A, we review how one can solve the theory analytically for various boundary conditions, and obtain exact results for the punctured and non-punctured models, in particular. In Appendix B, we show our results of the CLM for the punctured model using a different lattice definition of the topological charge.
Note added. When this paper was about to be completed, we encountered a preprint [67] on the arXiv, which investigates the same theory numerically by the density of states method. The exact results are reproduced for |θ| < π up to L = 24, which goes far beyond the reweighting method [64]. It is crucial to use an open boundary condition, which is similar in spirit to introducing a puncture in our work although the purpose is quite different.

2D U(1) lattice gauge theory with a θ term
In this section, we review 2D U(1) gauge theory with a θ term and discuss how to define it on a lattice.
In the continuum 2D U(1) gauge theory on a Euclidean space, the action for the gauge field A µ (x) (µ = 1, 2) is given by where g is the gauge coupling constant and F µν is the field strength defined as We add a θ term S θ = −iθQ (2.3) in the action, where Q is the topological charge defined by which takes integer values if the space is compact. We put this theory on a 2D torus, which is discretized into an L × L periodic lattice with the lattice spacing a. On the lattice, we define the link variables U n,µ ∈ U(1), where n labels the lattice site as x µ = an µ . We also define the plaquette which is a gauge invariant object. Here we write U −1 n,µ instead of U † n,µ , which will be important later in applying the CLM, where we complexify the dynamical variables respecting holomorphicity.
The lattice counterpart of the field strength (2.2) can be defined as where we take the principal value for the complex log; namely log z = log |z| + i arg z with −π < arg z ≤ π. Since the plaquette can then be written in terms of F n,µν as the lattice counterpart of the gauge action (2.1) can be defined as which approaches in the continuum limit up to an irrelevant constant with the identification (2.10) In the present 2D U(1) theory, the topological charge can be defined as which gives an integer value even at finite a. This can be proved easily by noting that n P n = 1 since each link variable appears twice in this product with opposite directions. We call this definition (2.11) the "log definition". As an alternative definition, we consider which approaches (2.4) in the continuum limit recalling (2.7). Note, however, that the topological charge defined on the lattice in this way can take non-integer values in general before taking the continuum limit. We call this definition (2.12) the "sine definition". Thus the lattice theory is given by where S g is given by (2.8) and S θ is given by (2.3) with Q defined either by (2.11) or by (2.12). Since this theory is superrenormalizable, we can take the continuum limit a → 0 with fixed g, which is set to unity throughout this paper without loss of generality. In this unit, the physical volume of the 2D torus is given by where we have used (2.10).

Applying the CLM to the 2D U(1) gauge theory
Since the θ term is purely imaginary in general, it makes Monte Carlo studies of gauge theories extremely difficult due to the sign problem. We overcome this problem by using the complex Langevin method (CLM) [14,15,16,17,18,19], which extends the idea of stochastic quantization to systems with complex actions. In this section, we discuss how to apply the CLM to 2D U(1) gauge theory with a θ term and show some results, which suggest that a naive implementation of the method fails.

the complex Langevin equation
The first step of the CLM is to complexify the dynamical variables. In the present case of U(1) gauge theory, we extend the link variables U n,µ ∈ U(1) to U n,µ ∈ C \ {0}, which corresponds to extending the gauge field A µ (x) ∈ R to A µ (x) ∈ C in the continuum theory. Then we consider a fictitious time evolution U n,µ (t) of the link variables governed by the complex Langevin equation where η n,µ (t) is a real Gaussian noise normalized by η n,µ (s)η k,ν (t) = 2δ n,k δ µ,ν δ s,t . The term D n,µ S is the drift term defined by first for the unitary link variables U n,µ (t), and then it is defined for the complexified link variables U n,µ (t) by analytic continuation in order to respect holomorphicity. Using the action (2.13), we obtain D n,µ S = D n,µ S g + D n,µ S θ , where the first term is given as The second term D n,µ S θ depends on the definition of the topological charge. If one uses the log definition (2.11), Eq. (3.2) for the θ term becomes a δ-function, which vanishes identically except for configurations with P n = −1 for some n, reflecting the topological nature of the definition. Such configurations are precisely the ones that appear when the topology change occurs within the configuration space of U n,µ . It is not straightforward to extend such a term to a holomorphic function of U n,µ .
On the other hand, if one uses the sine definition (2.12), the drift term becomes which may be viewed as an approximation of the δ-function mentioned above. Moreover, it can be readily extended to a holomorphic function of U n,µ . For this reason, we use the sine definition for the non-punctured model. The criterion [19] for the validity of the CLM states that the histogram of the drift term should fall off exponentially or faster. There are two cases in which this criterion cannot be met. The first case occurs when the configuration comes close to the poles of the drift terms (3.3), (3.4), which correspond to configurations with P n = 0 for some n. If this happens during the Langevin process, there is a possibility of violating the criterion. This problem is called the singular-drift problem [68,69], which was found first in simple models [70,71]. In the present model, the same problem is caused also by approaching configurations with |P n | = ∞ for some n, which are related to the poles by the parity transformation.
The second case occurs when the dynamical variables make large excursions in the imaginary directions [16]. This problem is called the excursion problem. In the present model, this corresponds to the situation in which the link variables have absolute values |U n,µ | far from unity.
Both the singular-drift problem and the excursion problem can occur because the link variables U n,µ are not restricted to be unitary in the CLM. In order to avoid these problems, it is important to perform the gauge cooling, which we explain in the next section.

gauge cooling
The idea of gauge cooling [66] is to reduce the non-unitarity of link variables as much as possible by making gauge transformations corresponding to the complexified Lie group after each step (3.1) of the Langevin process. This procedure can be added without affecting the argument for justifying the CLM as demonstrated explicitly in Refs. [18,19]. Recently, the mechanism of the gauge cooling for stabilizing the complex Langevin simulation has been investigated [72].
The deviation of the link variables from U(1) can be defined by the unitarity norm The gauge cooling reduces this quantity by a complexified gauge transformation, which is determined as follows.
First we consider an infinitesimal gauge transformation where ǫ n ∈ R. The change of the unitarity norm due to the transformation is given by where G n is defined as Therefore, we find that the unitarity norm is reduced most efficiently by choosing ǫ n ∝ −G n . Using this result, we consider a finite gauge transformation U n,µ → g n U n,µ g −1 n+μ ; g n = e −αGn , (3.9) which makes the unitarity norm depending on α in (3.9). We search for an optimal α that minimizes (3.10). Note here that it is typically a small number since the gauge cooling is performed after each step of the Langevin process. We therefore expand Eq. (3.10) with respect to α up to the second order and obtain the value of α that minimizes it as We repeat this procedure until the unitarity norm changes by a fraction less than 10 −5 .

adaptive stepsize
When we solve the complex Langevin equation in its discretized version (3.1), it occasionally happens that the drift term becomes extremely large, in particular during the thermalization process. This causes a large discretization error, which either makes the thermalization slow or destabilizes the simulation. We can avoid this problem by using a small stepsize ∆t, but the computational cost for a fixed Langevin time increases proportionally to (∆t) −1 and the calculation becomes easily unfeasible. The adaptive stepsize [73] is a useful technique, which amounts to reducing the stepsize only when the drift term becomes large. In our simulation, we measure the magnitude of the drift term defined as at each step, and choose the Langevin stepsize ∆t in (3.1) as where ∆t 0 is the default stepsize, and v 0 is the threshold for the magnitude of drift term.
In the present work, the default stepsize is set to ∆t 0 = 10 −5 , and the threshold is set to v 0 = 2β, considering a bound u ≤ 2β for θ = 0, where the CLM reduces to the real Langevin method. The measurement of the observables should be made with the same interval in terms of the Langevin time but not in terms of the number of steps.

results with the naive implementation
In this section, we present our results obtained by the CLM, which is implemented naively using the non-punctured model explained above as opposed to the punctured model, which we use later. As for the definition of the topological charge, we adopt the sine definition (2.12) for the reason given in Section 3.1.
We have performed simulations at various θ for (β, L) = (3, 10), (12,20) corresponding to a fixed physical volume V phys ≡ L 2 /β = 10 2 /3. Below we show our results only for θ = π, where the sign problem becomes severest, but the situation is the same for all values of θ.
In Fig. 1(Right), we plot the histogram of Re Q sin obtained by the CLM for (β, L) = (12, 20) with θ = π, which has a sharp peak at Re Q sin ∼ 0. In the same figure, we also plot the exact result for (β, L) = (12, 20) with θ = 0 for comparison, which exhibits a few sharp peaks at integer values within the range −2 Re Q sin 2. From these two plots, we conclude that the transitions between different topological sectors are highly suppressed in the simulation, which causes a problem with the ergodicity.
This occurs also at θ = 0 for large β, and it is called the "topology freezing problem" in the literature. In fact, the results one obtains by simulations suffering from this problem correspond to the expectation values restricted to the topological sector specified by the initial configuration. This is true for both θ = 0 and θ = 0. In this case, however, the effect of the θ term cancels between the numerator and the denominator of the expectation values, and the calculation essentially reduces to that of the real Langevin method at θ = 0.
For (β, L) = (3, 10) with θ = π, on the other hand, the histogram of Re Q sin obtained by the CLM has broad peaks that overlap with each other, which looks similar to the exact result for (β, L) = (3, 10) with θ = 0. This implies that the topology freezing problem does not occur for (β, L) = (3, 10). See also Fig. 3.
Below we define the observables we investigate in this paper. First, we define the average plaquette by Hereafter, V denotes the number of plaquettes in the action, which is V = L 2 for the nonpunctured model and V = L 2 − 1 for the punctured model we define in Section 4.1. The topological charge density is defined by 15) which is zero at θ = 0 and purely imaginary for θ = 0. Finally, the topological susceptibility  is defined by which is real for all θ. In fact, the topological susceptibility χ is related to the topological charge density (3.15) through Note, however, that this relation can be violated if the CLM fails to calculate the expectation values correctly. In Fig. 2, we show the results obtained by the CLM for the non-punctured model. We also plot the exact results for comparison, which are derived in Appendix A.4. In the left column, we present our results for (β, L) = (3, 10), which suffer from the incorrect convergence, whereas in the right column, we present our results for (β, L) = (12,20), which suffer from the topology freezing problem. In either case, our results do not reproduce the exact results as anticipated. Note that our results at θ = 0 agree with the exact results for (β, L) = (3, 10) but not for (β, L) = (12,20). This is because the topology freezing problem occurs for large β even at θ = 0, where the sign problem is absent.
Thus we find that the CLM with the naive implementation fails for both (β, L) = (3, 10) and (β, L) = (12, 20) for different reasons. For (β, L) = (3, 10), the topology change occurs but the criterion for correct convergence is not satisfied due to the large drifts. For (β, L) = (12,20), the criterion for correct convergence is satisfied, but the ergodicity is violated due to the topology freezing problem. We have searched for a parameter region in which neither of the problems occur, but we could not find one. In fact, we will see in the next section that these problems are related to each other at least in the present model.

the appearance of large drifts and the topology change
In this section, we provide more in-depth discussions on the relationship between the appearance of large drifts and the topology change in the non-punctured model. Let us first recall that the drift terms are given by (3.3) and (3.4), which depend on P n . When β is large, the gauge action S g favors configurations with P n ∼ 1 for all n, which implies that the drift terms are small.
On the other hand, the notion of topological sectors can be defined by the real part of (2.11), which takes integer values, even for complexified configurations that are generated in the CLM. In order for a transition between different topological sectors to occur, one of the plaquettes has to cross the branch cut; namely the phase of the plaquette has to jump from −π to π or vice versa. When this occurs, large drift terms can appear as can be seen We have also confirmed that the large drift term appears for the link variables composing the plaquette that crosses the branch cut.
In order to understand this observation better, we focus on a particular link variable U k,1 , and consider the corresponding drift term, which depends on the plaquettes P k and P k−2 sharing the link. For simplicity, we set P k−2 = 1 and consider the drift term v as a function of where we have defined a complex parameter φ by φ = −i log P k . A large drift appears when |Im φ| → ∞. In Fig. 4(Left), we plot the drift term as a flow diagram for β = θ = 1.
Considering that the contribution of the drift term v to the change of φ at a Langevin step is given by ∆φ = −v∆t, we actually plot (−v) in the complex φ plane.
In what follows we assume that β > θ/2π. Then we find from Eq. (3.18) that there are two fixed points corresponding to v = 0. One is φ = 0 and the other is φ = i log[(θ/2π + β)/(θ/2π − β)], which is close to ±π for β ≫ θ/2π. As one can see from Fig. 4(Left), the fixed point φ = 0 is attractive, which confirms that P k tends to become unity when β is large. The other fixed point φ ∼ ±π is repulsive, and the magnitude |v| grows exponentially as one flows away in the imaginary direction; See Fig.4(Right). As we mentioned above, when the transition between topological sectors occurs, one of the plaquettes crosses the branch cut, which corresponds to Re φ = ±π in the flow diagram. When this happens, the configuration can flow in the imaginary direction, which causes a large drift.

Introducing a puncture on the 2D torus
Since the problem we encounter in the previous section occurs due to the topological nature of the θ term, a simple remedy would be to change the topology of the base manifold to a noncompact one. Here we consider introducing a puncture on the 2D torus. Once we introduce a puncture, the drift term D n,µ S θ with the log definition of the topological charge has nonzero contributions for the link variables surrounding the puncture, which enable us to include the effect of the θ term correctly in the CLM as we will see in Section 5. Therefore, for the rest of this paper, we basically use the log definition to simplify our discussions. Unlike the non-punctured model, the topological charge is no more restricted to integer values, and it can be changed freely.
Since the puncture affects the theory only locally, its effect is expected to die out in the infinite volume limit for |θ| < π as we demonstrate explicitly in this section using the exact results. Thus unless we are interested in a theory with a finite volume, the punctured model is as good as the original model, the difference simply being a different choice of "boundary conditions". In fact, we will see that the non-punctured model has slow convergence to the infinite volume limit for θ ∼ π, which is not the case in the punctured model.

defining the punctured model on the lattice
There are various ways to introduce a puncture on the periodic lattice. Here we consider removing a plaquette as a simple choice. More precisely, we define the punctured model by removing one plaquette, let say P K , from the sum appearing in (2.8) and (2.11) when we define the action (2.13).
As an alternative method, we have also tried introducing a slit at a particular link, which amounts to duplicating the corresponding link variable and including each of them in the plaquettes that share the link. The results turn out to be qualitatively the same as the ones obtained by removing a plaquette. There are, of course, many others, but in any case, one can obtain exact results for a finite lattice as we explain in Appendix A, and using them, one can demonstrate explicitly that the punctured model is equivalent to the original non-punctured model in the infinite volume limit for |θ| < π.

equivalence in the infinite volume limit
In this section, we show the equivalence of the non-punctured model and the punctured model in the infinite volume limit. Here we use the log definition of the topological charge, but a similar statement holds as far as the same definition is used for the two models. 2 The partition function for the non-punctured model is given by (See Appendix A.2 for derivation.) for finite V = L 2 , where the function I(n, θ, β) is defined by Let us take the infinite volume limit V → ∞, in which the sum over n in (4.1) is dominated by the term that gives the largest absolute value |I(n, θ, β)|. This corresponds to the n that minimizes | θ 2π − n|. Thus in the infinite volume limit, the free energy is obtained as whereθ is defined byθ = θ − 2πk with the integer k chosen so that −π <θ ≤ π.  for finite V = L 2 − 1, which implies that the free energy is actually V independent. Hence all the observables that can be derived from it has no finite size effects. Note also that this model does not have the 2π periodicity in θ. By comparing (4.3) and (4.5), one can see that the two models are equivalent in the infinite volume limit for |θ| < π.
The observables defined in Section 3.4 can be calculated for the two models using (4.1) and (4.4) by numerical integration (See Appendix A.4 for the details.). In Fig. 5, we plot the average plaquette (Top) defined by (3.14), the imaginary part of the topological charge density (Middle) defined by (3.15) and the topological susceptibility (Bottom) defined by (3.16) for L = 10 (Left) and L = 20 (Right), respectively, with the same β = 12. The results for the two models tend to agree as L increases for |θ| < π.
We can evaluate the free energy (4.5) for the punctured model more explicitly for large β, which is relevant in the continuum limit. By integrating over φ in Eq. (4.2) as From this, we can obtain various observables for the punctured model as for finite V , which explains the θ dependence observed in Fig. 5. From Fig. 5, we also find that the results for the non-punctured model have sizable finite volume effects, in particular around θ ∼ π, which is absent in the punctured model. While the volume independence of the punctured model may well be peculiar to the present 2D gauge theory case, the advantage of the punctured model compared with the non-punctured model from the viewpoint of finite volume effects may hold more generally.

Application of the CLM to the punctured model
In this section, we apply the CLM to the punctured model using the log definition Q log of the topological charge. Our results reproduce the exact results discussed in the previous section as long as we are close enough to the continuum limit. We also show that the topology freezing problem is circumvented without causing large drifts thanks to the puncture.

the drift terms for the punctured model
We have discussed the drift terms in the non-punctured model in Section 3.1. For the punctured model, we only have to modify the drift terms for the four link variables surrounding the puncture; i.e., U K,1 , U K+2,1 , U K,2 and U K+1,2 . Thus we obtain where we have ignored the issue of δ-function discussed in Section 3.1. This is justified if all the plaquettes in the action never cross the branch cut; i.e., |Im log P n | ≤ π − ǫ for ∀n = K with a strictly positive ǫ during the Langevin simulation. We will see that this assumption is justified at sufficiently large β in Section 5.3. Note that the drift term from the θ term appears only for the link variables surrounding the puncture, and it is actually a constant independent of the configuration. While these properties are peculiar to the log definition Q log , similar properties hold also for the sine definition Q sin at large β, where all the plaquettes P n approach unity except for P K , which corresponds to the puncture. We discuss the case with the sine definition in Appendix B, where we see that the obtained results are qualitatively the same as those obtained with the log definition.

the θ dependence of the partition function
As we have seen in Section 4.2, the punctured model is equivalent to the non-punctured model in the infinite volume limit for |θ| < π, beyond which the equivalence ceases to hold. In particular, the punctured model does not have the 2π periodicity in θ, which exists in the non-punctured model. In order to understand this point better, we discuss the θ dependence of the partition function in this section. Let us first note that the partition function for arbitrary θ is related to the topological charge distribution ρ(q) for θ = 0 through Fourier transformation as Therefore, the absence of the 2π periodicity in θ in the punctured model is directly related to its property that the topological charge can take non-integer values even if we use the log definition Q log . Going beyond the fundamental region −π < θ ≤ π simply amounts to probing the fine structure of the topological charge distribution ρ(q), which is irrelevant in the infinite volume limit. By making an inverse Fourier transform, we can obtain the topological charge distribution ρ(q) for θ = 0 as We calculate this quantity for the punctured model by the CLM for θ = 0. In Fig. 6 Note that the calculation actually reduces to that of the real Langevin method due to the absence of the sign problem for θ = 0. We therefore have no concerns about the criterion for correct convergence here. While the sign problem is absent for θ = 0, the topology freezing problem can still be an issue for large β. The agreement we see for (β, L) = (12,20) confirms that this problem is resolved in the punctured model at least for θ = 0.

validity of the CLM
In this section, we discuss the validity of the CLM for the punctured model. Fig. 7(Left) shows the histogram of the drift term for (β, L) = (3, 10) and (β, L) = (12, 20) with θ = π, which are the parameters used in Section 3.4 for the non-punctured model. We find that the criterion is satisfied for (β, L) = (12, 20) but not for (β, L) = (3, 10), similarly to the situation in the non-punctured model. The difference from the non-punctured model is seen, however, in Fig. 7(Right), where we show the histogram of Re Q log obtained by the CLM for (β, L) = (12, 20) with θ = π. (The result for (β, L) = (3, 10) looks quite similar to this plot.) It is widely distributed within the range −3 Re Q log 3, which is in sharp contrast to the plot in Fig. 1(Right) for the same (β, L) = (12,20) in the case of the non-punctured model. In fact, it turns out to be close 3 to the exact result obtained for the same (β, L) = (12, 20) with θ = 0, which is plotted in the same figure. Thus we find that the topology freezing problem at large β is circumvented in the punctured model and yet the CLM remains valid.
Next we discuss the reason why the punctured model can avoid the topology freezing problem without causing large drifts. The difference from the non-punctured model is that one of the plaquettes, P K , is removed from the action. Note that the topological charge Q log for the punctured model is given by where the first term is nothing but the topological charge defined for the non-punctured model, whose real part takes integer values. The second term has a real part which lies within the interval [− 1 2 , 1 2 ). Therefore it makes sense to define the "topology change" in the punctured model as the situation in which the real part of the first term changes by ±1. As we discussed in Section 3.5 for the non-punctured model, one of the plaquettes should inevitably cross the branch cut in order for the topology change to occur in the above sense. When β is large, this process is highly suppressed for all the plaquettes that are included in the action. In the non-punctured model, the topology freezing problem occurs precisely for this reason. However, in the punctured model, the particular plaquette P K is removed from the action, and therefore it can change freely even for large β. This is demonstrated in Fig. 8, where we plot the probability distribution of the phase of the plaquette P K as well as that of the other plaquettes P n (n = K) for (β, L) = (3, 10) (Left) and (β, L) = (12, 20) (Right). We find that the phase of the removed plaquette P K is almost uniformly distributed for both (β, L). On the other hand, the distribution of the phase of the other plaquettes depends on (β, L). It has a compact support for (β, L) = (12, 20) but not for (β, L) = (3, 10). In the former case, there is no distribution at the branch cut, which implies that the branch cut crossing of the plaquettes P n (n = K) does not occur at all. In the latter case, there is a small but finite distribution at the branch cut, which means that the value of β is not large enough to suppress the branch cut crossing of the plaquettes P n (n = K) completely. This is consistent with the fact that the histogram of the drift term has fast fall-off for (β, L) = (12, 20) but not for (β, L) = (3, 10) considering the discussion given in Section 3.5. While the flow diagram in Fig. 4(Left) is obtained for the sine definition of the topological charge, it looks similar for the log definition, which simply corresponds to setting θ = 0 in (3.18). Therefore, large drifts can appear when one of the plaquettes P n (n = K) crosses the branch cut, which indeed occurs for (β, L) = (3, 10) also for the punctured model. For (β, L) = (12,20), on the other hand, the topology change is made possible by allowing the removed plaquette P K to cross the branch cut freely, but all the plaquettes that are included in the action are forced to stay close to unity because of large β. This justifies our assumption that the issue of δ-function can be neglected in deriving the drift terms (5.1) and (5.2). Since the plaquette P K does not appear in the drift terms, it does not cause large drifts even if it crosses the branch cut. This makes it possible for the punctured model to avoid the topology freezing problem without causing large drifts. Let us next discuss how the unitarity norm (3.5) behaves in our complex Langevin simulations. As we can see from (5.1) and (5.2), the link variables surrounding the puncture have a drift term in the imaginary direction coming from the θ term. At each Langevin step, two of the link variables are multiplied by e θ∆t/2π and the other two are multiplied by e −θ∆t/2π so that the removed plaquette is multiplied by e 2θ∆t/π due to this drift term. Therefore, there is a danger that the magnitude of these four link variables increases or decreases exponentially and hence the unitarity norm (3.5) grows exponentially with the Langevin time.
In Fig. 9, we plot the history of the unitarity norm (3.5) for various θ with (β, L) = (5, 16). Similar results are obtained for other (β, L). (Here and for the rest of this subsection, we restrict ourselves to the parameter sets, for which the histogram of the drift term has fast fall-off.) Indeed we observe an exponential growth at early stage, but the unitarity norm actually saturates to a constant depending on θ at sufficiently long Langevin time. This saturation occurs since the non-unitarity of the four link variables surrounding the puncture propagates to all the other link variables on the lattice due to the interaction caused by the gauge action S g , which tries to make each plaquette except the removed one close to unity. We find that thermalization of various observables can be achieved only after the saturation of the unitarity norm.
In fact, we find that the unitarity norm is not distributed uniformly on the lattice due to the existence of the puncture, as is also expected from the above discussion. In order to see this, we define the "local unitarity norm" by which is an average of the unitarity norm for the four link variables surrounding each plaquette P n . The unitarity norm defined by (3.5) is simply an average of N (n) over all the plaquettes including the removed one; namely N = 1 L 2 n N (n). In Fig. 10(Left), we plot this quantity N (n) against n = (n 1 , n 2 ) for (β, L) = (12,20) with θ = π, where the puncture is located at n = K = (10, 10). We observe a sharp peak at the puncture, which goes up to N (K) ∼ 6 × 10 3 . The plaquettes adjacent to the puncture have a local unitarity norm ∼ 1.5 × 10 3 . This implies that the unitarity norm is mostly dominated by the four link variables surrounding the puncture.
The local unitarity norm N (K) at the puncture depends not only on θ but also on β and L. In Fig. 10(Right), we plot this value against x = |θ|V phys = |θ|L 2 /β for various θ, β and L. All the data can be fitted to a single curve N (K) = x p exp(c 1 x + c 2 ), which reveals  an exponential behavior at large x.
What actually matters for the validity of the CLM is not so much the local unitarity norm N (n) as the absolute value of each plaquette P n , which we plot in Fig. 11 against n = (n 1 , n 2 ) for the same parameters as in Fig. 10(Left). The absolute value of P K corresponding to the removed plaquette is close to ( N (K)) 4 ∼ 3.6 × 10 7 , which implies that |U K,1 |, |U K+1,2 |, |U −1 K+2,1 | and |U −1 K,2 | are close to N (K). Except for this removed plaquette, the absolute value of the plaquette deviates only slightly from unity due to large β.
In fact, this deviation of |P n | from unity for n = K has a physical meaning since Im Q log = − 1 2π n =K log |P n | as one can see from (5.5). From the exact result (4.9) obtained at large β, we find that |P n | ∼ e −θ/(2πβ) for n = K, which is ∼ 0.96 for θ = π and β = 12 in agreement with the value observed in Fig. 11. If we flip the sign of θ, which corresponds to the parity transformation, we find that |P n | → |P n | −1 for all n.
Note also that P K does not appear in the drift term, which implies that its absolute value can become large without causing large drifts. We have confirmed that the criterion for correct convergence is satisfied for sufficiently large β, and indeed the exact results for various observables can be reproduced correctly as we will see in the next section. This remains to be the case even for large θ and/or large V phys , where the unitarity norm becomes large. 4 Thus the present model provides a counterexample to the common wisdom that the CLM fails when the unitarity norm becomes large.

results for the observables
In this section, we calculate the observables for the punctured model by the CLM and compare our results with the exact results derived in Appendix A.4. Let us recall that, in the definitions (3.14), (3.15) and (3.16), V denotes the number of plaquettes in the action, which is V = L 2 − 1 for the punctured model. In contrast, we define the physical volume V phys by Eq. (2.14) not only for the non-punctured model but also for the punctured model, which simplifies the relationship between β and L for fixed V phys .
In Fig. 12, we show our results for the average plaquette w (Top), the topological charge (Middle) and the topological susceptibility (Bottom) against θ for (β, L) = (3,10) and (12,20) in the left and right columns, respectively, which correspond to a fixed physical volume V phys ≡ L 2 /β = 10 2 /3. The exact results obtained for the same model with the same parameter sets are also shown for comparison. We find from our results for the average  plaquette that the exact results are reproduced for (β, L) = (12,20), but there is slight deviation for (β, L) = (3,10). This is consistent with our observation in Section 5.3 that the condition for correct convergence is met for (β, L) = (12, 20) but not for (β, L) = (3, 10). For the topological charge, we find that our results reproduce the exact results not only for (β, L) = (12, 20) but also for (β, L) = (3, 10). The same holds for the topological susceptibility. We consider that the agreement observed here for (β, L) = (3, 10) is accidental, though, since the condition for correct convergence is not satisfied. The fact that the results of the CLM for the punctured model with (β, L) = (3, 10) is not as bad as those for the non-punctured model with the same (β, L) shown in Fig. 2(Left) can be understood by considering that the effect of the θ term is included correctly by the drift terms for the link variables composing the removed plaquette, but it is only the infrequent branch cut crossing of the other plaquettes that spoils the validity of the CLM.

Summary and discussions
In this paper, we have made an attempt to apply the CLM to gauge theories with a θ term. As a first step, we applied the CLM to the 2D U(1) case, which is exactly solvable on a finite lattice with various boundary conditions. We find that a naive implementation of the method fails due to the topological nature of the θ term.
While the gauge configurations are complexified in the CLM, one can still define the notion of topological sectors by Re Q log ∈ Z. When a transition between different topological sectors occurs, one of the plaquettes has to cross the branch cut inevitably, which causes the appearance of large drift terms. This indeed happens at small β, where we find that the criterion for correct convergence of the CLM is not satisfied. Increasing β makes all the plaquettes close to unity. The large drift terms do not appear in this case, and the criterion for correct convergence of the CLM is satisfied. However, the topology change does not occur during the simulation and the ergodicity is violated. This is analogous to the topology freezing problem, which is known to occur for θ = 0. The results obtained in this case correspond to the expectation values for an ensemble restricted to a particular topological sector specified by the initial configuration.
In order to avoid this problem, we have considered the punctured model, which can be obtained by removing one plaquette from the action, both from the gauge action and from the θ term. While the quantity Re Q log is no more restricted to integer values, we can still formally classify the complexified configurations into "topological sectors" by adding back the contribution of the removed plaquette to Re Q log . Even for large β, the removed plaquette can cross the branch cut easily, which results in frequent transitions between different "topological sectors". Note also that, as far as β is sufficiently large, all the other plaquettes are close to unity, and hence large drift terms do not appear. Thus the criterion for correct convergence of the CLM can be satisfied by simply approaching the continuum limit without causing the topology freezing problem. Indeed our results obtained by the CLM for the punctured model reproduce the exact results even at large θ.
In the case of the punctured model, the drift term from the θ term appears only for the link variables composing the removed plaquette, and it is given by ±i θ 2π , which causes multiplication by a constant factor e ∓∆t θ 2π to these link variables at each Langevin step. The local unitarity norm of these link variables grows exponentially at early Langevin times, but it saturates at some point to some constant, which increases exponentially for large |θ|V phys . We have seen that the CLM works perfectly even in this situation as far as β is sufficiently large. This provides a counterexample to the common wisdom that the CLM fails when the unitarity norm becomes large. Thus our results also give us new insights into the method itself.
The punctured model is actually equivalent to the non-punctured model in the infinite volume limit for |θ| < π. In that limit, the topological charge can take arbitrarily large values, so the discretization of Q to integers is no more important. This equivalence has been confirmed explicitly by obtaining exact results for the punctured model. In fact, the exact results also reveal the absence of finite volume effects in the punctured model as opposed to the non-punctured model, which exhibits sizable finite volume effects around θ ∼ π. It is conceivable that the smearing of the topological charge somehow results in the reduction of finite volume effects. If so, a similar conclusion should hold more generally.
We are currently working on the application of the CLM to the 4D gauge theory with a θ term. (See Ref. [74] for an earlier attempt.) Some preliminary results obtained in the SU(2) case look promising. We plan to investigate, in particular, the interesting phase structure around θ = π predicted in Ref. [6].

A Derivation of the exact result
In 2D lattice gauge theory, we can obtain the partition function explicitly on any manifold at finite lattice spacing and finite volume [62], from which various observables can be obtained. In this section, we review the derivation using the so-called K-functional [61].

A.1 the K-functional
Let us consider a lattice gauge theory with a θ term on a 2D lattice manifold M. Here we take the gauge group to be U(N), which is a generalization of U(1) considered so far. Note that the topology of the gauge field becomes trivial for SU(N) in 2D gauge theories.
As a building block for evaluating the partition function, we define the K-functional K A for the region A ⊂ M defined by [61] where the integral goes over the link variables inside A leaving out those on the boundary C. (See Fig. 13.) The action S A in Eq. (A.1) is given by where the sum goes over the plaquettes P i included in the region A. Here we use the log definition (2.11) of the topological charge, but the results for the sine definition (2.12) can be obtained in a similar manner as we mention at the end of Section A.4. The K-functional depends on the link variables on the boundary C = ∂A, but due to the gauge invariance, it actually depends only on which is a consecutive product of link variables along the loop C. The choice of the starting point of the loop C does not matter since a different choice simply corresponds to making a gauge transformation of Γ, which leaves the K-functional invariant. We can calculate the K-functional for any A by gluing the K-functional for a single plaquette P , which is nothing but K(P ) = exp Tr β 2 P + P −1 + θ 2π log P . Note here that (A.4) is a function of the group element P ∈ U(N), which is invariant under It is known that any function having this property can be expressed by the so-called character expansion K(P ) = r λ r χ r (P ) , (A. 6) which is analogous to the Fourier expansion for periodic functions. Here χ r (P ) is the group character, which is defined by the trace of P for an irreducible representation r, and it satisfies the orthogonality relation Using this relation, the coefficient λ r in the expansion (A.6) can be readily obtained as As an example, let us obtain the K-functional K 2×1 for a 2 × 1 rectangle by gluing two neighboring plaquettes P 1 = U 1 Ω and P 2 = Ω −1 U 2 as shown in Fig. 14. The group elements U 1 and U 2 are the products of three link variables, and Ω represents the link variable shared Figure 14: The K-functional K 2×1 for a 2 × 1 rectangle is obtained by considering the K-functional for the two plaquettes P 1 = U 1 Ω and P 2 = Ω −1 U 2 , which are glued together by integrating out the shared link variable Ω.
by P 1 and P 2 . Integrating out the shared link variable Ω, we get where d r = χ r (1) is the dimension of the representation r and we have used a formula Iterating this procedure, we obtain the K-functional for any simply connected region A as where |A| is the number of plaquettes in A, and Γ is defined by (A.3).
In the U(1) case, the representation can be labeled by the charge n ∈ Z, and the dimension of the representation is d n = 1 for ∀n ∈ Z. Since the character for the plaquette P = e iφ is given by χ n (P ) = e inφ , the K-functional for a single plaquette (A.6) reduces to where the coefficient λ n is a function of θ and β given explicitly as λ n ≡ I(n, θ, β) using (A.4) with P = e iφ . This function reduces to the modified Bessel function of the first kind for θ = 0. The character expansion in the U(N) case is more complicated, so we only show the end results referring the reader, for instance, to the appendix of Ref. [75] for the details. The representation of the U(N) group is labeled by N integers ρ = (ρ 1 , ρ 2 , · · · , ρ N ) ∈ Z N (A.14) satisfying ρ i ≥ ρ i+1 , and the dimension of the representation ρ can be calculated by The coefficient λ ρ in (A.6) that corresponds to the representation ρ is expressed as a determinant λ ρ = det M(ρ, θ, β) , (A. 16) where the matrix M(ρ, θ, β) is given as which may be viewed as a generalization of (A.13).

A.2 partition function for the non-punctured model
Let us evaluate the partition function for the 2D U(N) lattice gauge theory on a torus. For that, we first consider the K-functional K L 1 ×L 2 for a rectangle composed of V = L 1 × L 2 plaquettes, which can be expressed as (A.11). As is shown in Fig. 15, we identify the top and bottom sides represented by U −1 and U, respectively, and identify the left and right sides represented by W −1 and W , respectively. Integrating out the group elements U and W , we obtain the partition function for the non-punctured model as As one can see from (A.13), the integral I(n, θ, β) has a property I(n, θ + 2πk, β) = I(n − k, θ, β) for ∀k ∈ Z, which guarantees the 2π periodicity of (A.20) in θ. Let us consider taking the V → ∞ and β → ∞ limits simultaneously with fixed V phys ≡ V /β, which corresponds to the continuum limit. In this limit, the integral (A.13) can be evaluated as Plugging this into (A.20), we obtain omitting the divergent constant factor.

A.3 partition function for the punctured model
Let us extend the calculation in the previous section to the punctured model. First, we calculate the K-functional for a rectangle with a hole shown in Fig. 16, which we divide into two regions A 1 and A 2 by cutting along two segments Ω 1 and Ω 2 . The outer and inner boundaries of the rectangle are divided into two segments (U 1 , U 2 ) and (ω 1 , ω 2 ), respectively. Then, the K-functional for each region is given, respectively, as By gluing the two regions A 1 and A 2 together at Ω 1 and Ω 2 , we obtain the K-functional for the rectangle with a hole as Let us introduce the group elements U and W for the outer boundary as we did in Fig. 15 so that U 1 U 2 = UW U −1 W −1 , and define ω = ω 1 ω 2 for the inner boundary. Integrating out the group elements U and W , we obtain the K-functional for the punctured torus as Finally, we integrate out the link variables surrounding the puncture to get the partition Figure 16: The K-functional for a rectangle with a hole is obtained by gluing the two regions A 1 and A 2 . From this, the K-functional for the punctured torus is obtained similarly to what we did in Fig. 15. Integrating out the link variables surrounding the puncture, we obtain the partition function for the 2D U(N) gauge theory on a punctured torus.
function for the punctured model as where r = 0 corresponds to the trivial representation, which has d 0 = 1.
In the U(1) case, the partition function reduces to which does not have the 2π periodicity in θ. Let us consider taking the V → ∞ and β → ∞ limits simultaneously with fixed V phys ≡ V /β, which corresponds to the continuum limit. Similarly to the case of the non-punctured model discussed in Section A.2, we obtain omitting the divergent constant factor. This coincides with (A.22) in the V phys → ∞ limit for |θ| < π. Note, however, that the equivalence between the punctured and non-punctured models does not hold for finite V phys .

A.4 evaluation of the observables
We can evaluate the expectation values of various observables defined in Section 3.4 from the partition function derived above, namely (A.20) for the non-punctured model and (A.28) for the punctured model. Since the latter case is easier due to the absence of an infinite sum, we only discuss the former case in what follows. The average plaquette w defined by (3.14) is given as where we have defined A(n, θ, β) = ∂ ∂β log I(n, θ, β) Similarly, the topological charge density defined by (3.15) can be obtained from Note that I(n, θ, β) and the functions (A.31), (A.33) and (A.35) derived from it are all real-valued, and we can calculate them by numerical integration with sufficient precision. Also, when we evaluate the infinite sum in the expressions (A.30), (A.32) and (A.34), we have to truncate it at some n. Note here that |I(n, θ, β)| vanishes quickly as |θ/2π − n| increases. We can therefore evaluate the infinite sum with sufficient precision by keeping only a few terms when the lattice volume V is sufficiently large.

(B.2)
At large β, all the plaquettes except P K , namely the one that is removed, approach unity. The drift term from the θ term therefore vanishes for all the link variables except for those surrounding the puncture, which have constant drifts ±i θ 2π . Thus in the continuum limit, the drift terms for the sine definition agree with those for the log definition given by (5.1) and (5.2). This connection makes it easier to understand why we can safely ignore the issue of δ-function in the drift term for the log definition described in Section 3.1. It is therefore expected that the results of the CLM for the sine definition are essentially the same as those for the log definition for large β. In Fig. 17, we show our results for the punctured model with the sine definition for the same (β, L) as those in Fig. 7 with the log definition. For (β, L) = (12,20), we find that the histogram of the magnitude u of the drift term falls off rapidly, and that the histogram of Re Q sin obtained by the CLM is widely distributed within the range −3 Re Q sin 3. Hence the topology freezing problem is circumvented without causing large drifts similarly to the situation with the log definition.
On the other hand, for (β, L) = (3, 10), we find that the histogram of the magnitude u of the drift term falls off fast and that the condition for the validity of the CLM is satisfied unlike the case of the log definition. As a result, all the observables are in complete agreement with the exact results for all values of θ even with (β, L) = (3, 10). This can be seen from Fig. 18, where we show our results for the punctured model with the sine definition for the same values of (β, L) as the ones used in Fig. 12. For (β, L) = (1.92, 8) corresponding to the same V phys ≡ L 2 /β, however, we actually find that the histogram has a power-law tail similarly to the case of the log definition. Therefore, the difference between the two definitions is merely a small shift in the validity region of the CLM.
We also show the exact results for the punctured model with the log and sine definitions, which tend to agree as β is increased with fixed V phys ≡ L 2 /β, which corresponds to the continuum limit.