Simulating Yang-Mills theories with a complex coupling
Jan M. Pawlowski, Manuel Scherzer, Christian Schmidt, Felix P.G. Ziegler, Felix Ziesché
SSimulating Yang-Mills theories with a complex coupling
Jan M. Pawlowski,
1, 2
Manuel Scherzer, Christian Schmidt, Felix P.G. Ziegler, and Felix Ziesch´e Institut f¨ur Theoretische Physik, Universit¨at Heidelberg,Philosophenweg 16, D-69120 Heidelberg, Germany ExtreMe Matter Institute EMMI, GSI, Planckstraße 1, D-64291 Darmstadt, Germany Fakult¨at f¨ur Physik, Universit¨at Bielefeld, Postfach 100131, D-33501 Bielefeld, Germany. CP3-Origins and Danish IAS, Department of Mathematics and Computer Science,University of Southern Denmark, Campusvej 55, 5230 Odense M, Denmark (Dated: January 12, 2021)We propose a novel simulation strategy for Yang-Mills theories with a complex coupling, basedon the Lefschetz thimble decomposition. We envisage, that the approach developed in the presentwork, can also be adapted to QCD at finite density, and real time simulations.Simulations with Lefschetz thimbles offer a potential solution to sign problems in Monte Carlocalculations within many different models with complex actions. We discuss the structure of Gen-eralized Lefschetz thimbles for pure Yang-Mills theories with a complex gauge coupling β and showhow to incorporate the gauge orbits. We propose to simulate such theories on the union of thetangential manifolds to the relevant Lefschetz thimbles attached to the critical manifolds of theYang-Mills action. We demonstrate our algorithm on a (1+1)-dimensional U(1) model and discusshow, starting from the main thimble result, successive subleading thimbles can be taken into accountvia a reweighting approach. While we face a residual sign problem, our novel approach performsexponentially better than the standard reweighting approach. I. INTRODUCTION
The notorious sign problem hampers numerical sim-ulations of many interesting physical systems, rangingfrom high energy to condensed matter systems. A signproblem is faced in numerical calculations of statisticalmodels whenever the action becomes genuinely complex.Hence, standard Monte Carlo methods and in particularimportance sampling drastically lose their efficiency withincreasing lattice volume.Examples of theories with a sign problem include realtime calculations in lattice-regularized quantum field the-ories, i.e. lattice QCD in Minkowski space-time, latticeQCD with a nonzero vacuum angle θ , and lattice QCDwith a nonzero baryon chemical potential µ B . For thelatter two cases, many methods have been developedthat potentially circumvent or solve this problem in thecontinuum limit. These methods include reweighting[1, 2], Taylor expansions [3–5], analytic continuation frompurely imaginary chemical potentials [6, 7], canonicalpartition functions [8, 9], strong coupling/dual methods[10–13], the density of states method [14–16], and com-plex Langevin dynamics [17–20], see [21] for a review ofrecent developments. However, all these methods so farface severe limitations that restrict their applicability inthe continuum limit.In the past decade deformations of the original inte-gration manifold into the complex domain have been in-troduced, based on complex saddle points of the action,the Lefschetz thimbles [22, 23]. If these deformationsare chosen well, all physical expectation values obtainedfrom an oscillatory integral remain unchanged but thesign problem is drastically alleviated. Thus, a numericalevaluation of the theory is accessible. By definition, theimaginary part of the action (phase of the probability density) is stationary on the thimble. A first Lefschetzthimble algorithm in the context of the QCD finite den-sity sign problem was introduced in [24], for a recent re-view see [25]. Despite its great potential for beating thesign problem, simulations with Lefschetz thimbles haveto overcome the following intricacies:( i ) A parametrization of the thimble is a priori un-known and has to be obtained as the numerical so-lution of a flow-equation. The parametrization andthe necessary Jacobian of the variable transforma-tion are numerically demanding.( ii ) The curvature of the thimble manifold introducesa so-called residual sign problem through the Jaco-bian.( iii ) In most cases, one has to include relevant contri-butions from a large number of thimbles. Theirrelative weight may give rise to a further, residual,sign problem. In any case, the probability densitybecomes multi-modal.These shortcomings have been addressed by formulatingparticularly efficient methods for the calculation of theJacobian [26, 27], by optimizing the deformation of theintegration domain either based on a model ansatz [28,29] or by means of machine learning [30–32] and finally byusing tempering [33, 34] and other advanced strategies tofoster transitions between thimbles and take into accountcontributions from multiple thimbles [35–37].In the present work we analyze the structure of theLefschetz thimble decomposition of pure gauge theorieswith gauge groups U( N ) and SU( N ) and complex cou-pling β . We propose to sample on the tangential man-ifold attached to the thimble of the main critical point, i.e. the critical point with the smallest action value. As a r X i v : . [ h e p - l a t ] J a n R e < / ( P + P - ) > Re( β ) analyticReweighting (N=120000)32 Tangent Spaces (N=3000) FIG. 1. Comparison of the novel Takagi simulation methodwith standard reweighting on an 8 × β , while keeping Im ( β ) = 1constant. Both reweighting and our approach took about thesame amount of computing time. expected, we find a full hierarchy of critical points (sad-dle points) that have to be considered depending on thecoupling parameter β . We include successive subleadingsaddle point contributions via reweighting. The reweight-ing procedure is set up by using linear mappings from themain tangential manifold to the tangential manifold at-tached to the thimble of the target critical point.Our approach has the following advantages over theflow-based generalized Lefschetz thimble approach [38]:( A
1) During the sampling procedure, we neither have toflow our configuration nor have to calculate a Jaco-bian since it is constant on the tangential manifold.( A
2) As contributions from subleading thimbles aretaken into account by reweighting [36], there is noergodicity problem due to potential barriers be-tween thimbles. The reweighting procedure doesalso not introduce any overlap problem since criti-cal points are mapped onto critical points.The disadvantages are:( D
1) Sampling on the tangential manifold rather thanon the thimble itself does not completely resolvethe sign problem, however, there is no residual signproblem from the Jacobian [26].( D
2) The critical points need to be known.The advantages are clearly demonstrated for the bench-mark case of an U(1) gauge theory, as summarized inFig. 1. There we compare our results for the real partof the plaquette with that produced by the standardreweighting procedure. For the same numerical costs theerror bars of our results are reduced by orders of magni-tude. This will be discussed in more detail in this work,which is organized as follows: In Sec. II we derive nec-essary formulas for our setup. Emphasis is devoted todiscussing the critical manifolds of the Yang-Mills actionas well as the Lefschetz thimbles. Moreover we explainour update and reweighting procedures. In Sec. III weapply our algorithm to the case of a 2-dimensional U(1)gauge theory. In Sec. IV we discuss prospect of our ap-proach. Finally we conclude in Sec. V.
II. FORMULATIONA. Overview
We consider the standard discretization of the Yang-Mills action [39] S = β (cid:88) x (cid:88) µ<ν (cid:26) − N (cid:0) Tr P µ,ν ( x ) + Tr P − µ,ν ( x ) (cid:1)(cid:27) , (1)where P µ,ν ( x ) = U µ ( x ) U ν ( x + ˆ µ ) U − µ ( x + ˆ ν ) U − ν ( x ) de-notes the elementary plaquette in the ( µ, ν )-plane at lat-tice site x ∈ Λ where Λ ⊂ Z d . The summation is donesuch that each plaquette is counted with only one ori-entation. The link variables U µ ( x ) are elements of thegauge group, which we consider to be a Lie group and inparticular U( N ) or SU( N ).A sign problem is introduced by choosing general com-plex couplings β .We aim to update a given gauge field configurationsuch that the imaginary part of the action Im( S ) variesslowly and hence the remaining sign problem is mild.In order to achieve this goal the link variables are com-plexified, i.e. they are allowed to take values within thelarger group GL( N, C ) or SL( N, C ), respectively. By gen-eralized Picard-Lefschetz theory, there is a smooth mid-dle dimensional manifold connected to each complex crit-ical manifold (simply connected union of points, where ∇ S = 0) of the action on which Im( S ) stays constant.These manifolds are called Lefschetz thimbles. As al-ready stated, updates that stay on these thimbles arecomputational very demanding and are invariably globalupdates [24, 38, 40, 41]. To reduce the computational de-mand, it has been argued that one does not have to stayexactly on the thimble in order to reduce the sign prob-lem significantly [29, 42]. Any deformation of the originalintegration domain which has the correct asymptotic be-havior will be sufficient. Here we construct an integral de-formation, which is the union of all tangential manifoldsto the critical points. As for the pure gauge theory allcritical manifolds are known, this manifold is straightfor-ward to parametrize. After discussing the critical man-ifolds in Sec. II B and II C we discuss properties of Lef-schetz thimbles and tangential manifolds in Sec. II D andII E. Our algorithmic approach is presented in Sec. II F,II H and II G. B. The critical points of the Yang-Mills action
A Lefschetz thimble is generally defined to be the unionof flow lines generated by the steepest descent equationd U d t = − (cid:18) δSδU (cid:19) ∗ , (2)which end in a non-degenerate critical point of the action.The degeneracy of critical points due to gauge symme-try necessitates the application of Witten’s concept ofgeneralized Lefschetz thimbles [43, 3.3].The critical manifolds can be described in terms of pla-quette variables. This is seen by examining the gradientof the action ∂ x,κ,a S = − iβ N × (3)Tr (cid:34)(cid:32)(cid:88) κ<ν P κ,ν ( x ) − P − κ,ν ( x ) + P κ, − ν ( x ) − P − κ, − ν ( x ) − (cid:88) µ<κ P µ,κ ( x ) − P − µ,κ ( x ) + P − µ,κ ( x ) − P − − µ,κ ( x ) (cid:33) t a (cid:35) . Here t a are the generator matrices of the related Lie-Algebras su ( N ) , u ( N ). In our notation those are Her-mitian, and for N > t a t b ] = δ ab . The derivative with respectto the gauge link U in the direction of t a is defined as ∂ a f ( U ) := ∂∂ω f ( e iωt a U ) | ω =0 . Negative signs in the sub-script of the plaquette variables refer to reversed direc-tions in their orientation, e.g. P κ, − ν ( x ) = U κ ( x ) U − ν ( x + ˆ κ − ˆ ν ) U − κ ( x − ˆ ν ) U ν ( x − ˆ ν ) . A necessary condition for a critical configuration is a van-ishing gradient of the action.In the following we derive relations that constrainpossible plaquette values from a critical configuration.Eq. (3) vanishes ∀ a , if the matrix in round brackets isproportional to for plaquette values in SL( N, C ). Forplaquette values in GL( N, C ), the matrix has to be zero.For a proof, see the App. D. This criticality conditionyields relations for adjacent plaquettes sharing one link.Note that in d dimension one link is shared by 2( d − d = 2 case: For plaquette values inGL( N, C ), we can directly read off the relations P , ( x ) = P − , − ( x ) or P , ( x ) = − P , − ( x ) , and P , ( x ) = P − − , ( x ) or P , ( x ) = − P − , ( x ) . (4)If we assume that our critical configuration consists of commuting (Abelian) link variables, i.e. if link variablesare diagonal, the above relation simplifies to P , ( x ) = P , ( x − ˆ ν ) or P , ( x ) = − P − , ( x − ˆ ν ) , (5)with ν ∈ { , } . It follows that each critical configura-tion exhibits at most two distinct plaquettes values. Ifplaquettes take values in SL( N, C ) we get, e.g. , for a linkin ˆ1-direction( P , ( x ) − P − , ( x ) − P − , − ( x ) + P , − ( x )) = α , (6)for an arbitrary α ∈ C . For the Abelian case, this equa-tion reformulates again to a relation between adjacentplaquettes. Restricting this further to the original groupSU( N ), i.e. the maximal torus of SU( N ), we find(Im P , ( x ) − Im P , ( x − ˆ ν )) ii =(Im P , ( x ) − Im P , ( x − ˆ ν )) jj ∀ i, j, ν . (7)For d >
2, we obtain the same constraint on the imagi-nary parts of diagonal entries of adjacent plaquettes, butnaturally there are more adjacent plaquettes.From the periodic boundary condition we can derivefurther constraints. Under the assumption, that the rel-evant critical configurations are Abelian, the product ofall plaquettes in an arbitrary 2-dimensional hyperplaneΛ µν must be one, i.e. (cid:89) x ∈ Λ µν P µν ( x ) = 1 . (8) C. Locality and critical manifolds
We are left with a fairly large number of criticalconfigurations, which we will boil down to a set of basisconfigurations, getting the others by transpositions andsymmetry relations. The ultimate aim is to obtain ahomotopic covering of the original integration domain[U( N )] dV or [SU( N )] dV .As a guiding inspiration, we start with a discussionof the one-plaquette model, i.e. we just have only oneplaquette degree of freedom, as Eq. (1) is a sum of localplaquette terms. The action is defined as S = − β N Tr (cid:2) P + P − (cid:3) , (9)where an irrelevant constant has been omitted. The Liederivatives with respect to P are given by ∂ a S = − iβ N Tr (cid:2) ( P − P − ) t a (cid:3) . (10)Observations are:(i) Eq. (10) vanishes for self-inverse plaquettes P inU( N ). Therefore, the critical points consists onlyof matrices whose eigenvalues are +1 and − P ∈ SU( N ) Eq. (7) implies that the imagi-nary part of all eigenvalues must be identical. Forvanishing imaginary part we obtain the self-inverseelements of SU( N ). The constraint det[ P ] = 1 im-plies an even number of ( − N ( − ) . For non-vanishing imaginary part, thecenter elements of SU( N ) are solutions. For N ≥ ± w ≡ e − Re ( S ) may be exponentially suppressed. For the one-plaquette model we find the following hierarchy of criticalpoints: For P ∈ U( N ), we have S = − βN (cid:16) N − N ( − ) (cid:17) . (11)The importance of a critical point thus shrinks with thenumber of ( − S = − β , P = S = − β , P ∈ (cid:110) e i π k , k = 1 , (cid:111) S = − β , P ∈ { diag(1 , − , − , diag( − , , − , diag( − , − , } . (12)Critical points from the one-plaquette model can beused to construct certain critical configurations for thefull lattice theory. We pick critical plaquette values fromthe one-plaquette model and distribute them in accor-dance with Eq. (4) and (6) over the lattice. Each criticalconfiguration obtained in this way is one representativeof a critical manifold, which consists of its gauge orbitand additional zero modes of the action. For this sim-ple procedure we obtain an additional constraint fromEq. (8): The product of eigenvalues at every positionover every two-dimensional hyperplane must be one. Weare therefore limited to configurations, where the num-ber of ( − N ) case) in every hyperplane or in principal theiroverall product is one including additional roots of unity.Note however, that not all critical configurations canbe found by the above prescription. Recall that Eq. (4)restricts the plaquette values in a hyperplane at any givenposition to only two possible values. We can constructa further critical configuration by setting one (or more)plaquette values in that hyperplane to e i(cid:15) , while choosing e i ( π − (cid:15) ) for all remaining plaquettes in the hyperplane.Possible (cid:15) values are constrained by Eq. (8), and hence kπ + ( V − k ) (cid:15) = 2 πl , (13) where k is the number of plaquettes with the specific di-agonal entry e i ( π − (cid:15) ) . The 2 πl factor stems from the 2 π periodicity. For d = 2, N = 1, l is the actual topo-logical charge. It is constant on the thimble, since theanti-holomorphic gradient flow can be seen as an analyticcontinuation of the classical gradient flow, which leavesthe topological charge invariant [44]. Note especially for k (cid:54) = V /
2, that (cid:15) has to be a real number, so the criticalconfigurations are all in the original group space. Picard-Lefschetz theory tells us now, that if the tangent spaceof the thimble is not normal at this point to the originalgroup manifold, then the intersection number is non-zero(in our case one). We will see, that for Im β (cid:54) = 0, this isthe case and they all contribute. For SU( N ), we have toadd the constrain that the determinant equals one, effec-tively reducing the number of critical configurations. D. The Takagi decomposition and generalizedLefschetz thimbles
Next, we construct the tangent spaces at each criticalmanifold described Sec. II C. To that end we solve theTakagi equation H ∗ ξ ∗ = λξ with λ ∈ R . (14)Here H denotes the Hessian of the action evaluated onthe critical manifold. Modes ξ associated with positive λ ,point in the direction of the thimble. In the following werefer to such a mode as Takagi vector . Correspondingly,a mode ξ associated with a negative λ points in directionof the anti-thimble. It is called anti-Takagi vector (seee.g. [38]) The λ = 0 vectors do not change the action andare therefore referred to as zero-modes. They result forinstance from the gauge degrees of freedom. The Hessianof the action can be written as ∂ x,κ,a ∂ y,η,b S = β N c Tr (cid:88) U κ ( x ) ,U η ( y ) ∈ P ( P + P − )+ (cid:88) U − κ ( x ) ,U − η ( y ) ∈ P ( P + P − ) − (cid:88) U κ ( x ) ,U − η ( y ) ∈ P ( P + P − ) − (cid:88) U − κ ( x ) ,U η ( y ) ∈ P ( P + P − ) t b t a , (15)where the plaquettes appear in different orientations de-pending on the position of the referred link. Since byconstruction we have considered only those representa-tives of our critical manifolds which have diagonal links.They are Abelian and we can permute our variables tobring all generators to the right.For critical configurations, where the plaquettes areelements of the original group, the Hessian splits into areal matrix with a complex prefactor H = βM , whoseeigenvectors v and eigenvalues α can be computed. For α (cid:54) = 0, the solutions to Eq. (14) take the form ξ ( ± ) = (cid:115) ± sign( α ) β ∗ | β | v , (16)where the ξ (+) ( ξ ( − ) ) indicate the thimble (anti-thimble)directions. The α = 0 eigenvectors correspond to the λ = 0 solutions of Eq. (14) and can have an arbitrarycomplex prefactor.We observe that for most of our critical configurations.We will come to the exceptions in Sec. III B. The Hessiandoes not change under field transformations in the direc-tion of these zero modes. Consequently the critical man-ifold { U crit , µ ( x ) } is independent of those. Therefore, wecan deduce that the projection of the subspace spannedby its zero modes in the Lie algebra is the critical mani-fold itself. U crit µ ( x ) = U crit , µ ( x ) exp i (cid:88) k,a ˜ b k v x,µ,ak ( α = 0) t a (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ b k ∈ C , x ∈ Λ , µ = 1 , . . . , d . (17)The k index enumerates the different zero-eigenvectors v of matrix M . For ˜ b k ∈ R , this is a compact manifold.In complexified space with ˜ b k ∈ C , we have non-compactimaginary directions. Nevertheless, the manifold is stillcritical. To obtain the generalized Lefschetz thimble [43],we start with a compact submanifold (a cycle) of realdimension α = 0) = n (0) . Therefrom its tangent spaceis spanned. The compact submanifold is commonly called gauge orbit . At every point of this cycle, we use the α (cid:54) = 0 Takagi vectors to span the rest of the tangentspace. Since the latter is invariant under zero-modes apoint on this tangent space can be directly written interms of U µ ( x ) = U crit , µ ( x ) exp i (cid:88) k,a b k c k v x,µ,ak t a , with b k ∈ R and c k = (cid:40) , α k = 0 (cid:113) sign( α k ) β ∗ | β | , α k (cid:54) = 0 . (18)Here the summation index k runs over all eigenvectors.The real dimension n of the thimble thus splits into zeroand non-zero mode directions as n = n (0) + n (+) = dV N g ,where N g is the dimension of the Lie algebra. This con-struction can be generalized to more complicated Hes-sians, see App. A. FIG. 2. Upper plot: Schematic picture of the generalized Lef-schetz thimble, spanned on a compact submanifold as pro-posed by Witten [43, 3.3]. Lower plot: The thimble will ex-hibit infinitely many Riemann surfaces if we choose a non-compact submanifold from the critical cylinder. In both plotsthe critical cylinder is depicted in orange and the thimble inblue.
If we do not conform to that construction, e.g. , tiltthe real vectors by a complex factor, we loose homo-topy to the generalized Lefschetz thimble. The resultingmanifold is non-compact, see Fig. 2. However, the non-compact directions refer to zero-modes, i.e. they leavethe action invariant. This corresponds to multiplying thepartition sum with an infinite volume factor, which dropsout for all physical observables, when calculating expec-tation values. The prerequisite for these observables is,that they are independent from the zero-modes, i.e. thatthey are gauge-invariant observables.The global minimum among our critical manifolds isgiven by { P µν ( x ) = ∀ x, µ < ν } . Since this is a mini-mum in the non-complexified gauge theory, the matrix M is positive semidefinite. As the gauge field is constant,the complex prefactor is identical to all Takagi-modes.We are free to choose the same complex prefactor alsofor the real zero-modes. Since the eigenvectors of M span the whole R n space and since we chose the samecomplex prefactor for all directions, we can make a ba-sis transformation to the unit basis { e , . . . , e n } . Movingin one of these directions e i corresponds to a change ofa single (local) color degree of freedom. This basis canthus be used to construct a local update algorithm onthe generalized main Lefschetz thimble, see Sec. II H.As we are interested in a single homotopic covering ofour original compact integration domain, we need to limitthe tangent spaces. It turns out that the construction ofa continous manifold that is homotopic to the original in-tegration domain is the main conceptual difficulty in thisapproach. For the U(1) one-plaquette model, we have -1 0 1 2 3 4 -4 -3 -2 -1 0 1 I m ( P ) Re(P)ThimblesTMsU(1) -2-1.5-1-0.5 0 0.5 1 1.5 2 -3 -2 -1 0 1 2 3 φ I φ R ThimblesTMsU(1)
FIG. 3. Lefschetz thimbles and tangential manifolds (TMs) bounded by their intersection for the one-plaquette model at β = 1 + 3 i in group and angular (algebra) representation. only two critical points and tangent spaces. They give,when glued together at their intersections, a manifold ho-motopic to U(1), see Fig. 3. We introduce boundaries onthe main tangential manifold by identifying intersectionswith the tangential spaces of the other subleading criticalmanifolds. We will call these limited tangent spaces, tan-gential manifolds (TM) from now on. For the full latticetheory we discuss several possible choices of boundariesthroughout this work. E. Hierarchy of critical manifolds
The choice of critical manifolds introduces a naturalhierarchy which is reflected by the values of the action.Their importance decreases according to their weight fac-tor e − S with increasing action. Given our choice of criti-cal configurations with diagonal plaquettes from the orig-inal integration domain (U( N ) or SU( N )), we can easilyexpress the action in terms of the plaquette eigenvalues.We find S = β (cid:88) x (cid:88) µ<ν (cid:34) − N N (cid:88) k =1 cos( φ ( k ) µ,ν ( x )) (cid:35) , (19)where the φ ( k ) x is the angle of the k -th eigenvalue of theplaquette P µν ( x ). These critical action values are minimaof the attached thimbles, since the real part of the actionnaturally increases if one moves away from the criticalmanifolds for Re ( β ) > φ ( k ) x = 0 ∀ x, k ) defines the globalminimum of the action in this many TMs scenario. Con-sequently with increasing Re ( β ), certain TMs becomeexponentially suppressed and we obtain a pronouncedhierarchy with the main TM as leading order. For purelyimaginary values of β , every thimble contributes equally. F. Update algorithm on a TM with a Takagi basis
Next, we discuss possible sampling algorithms that arerestricted to a single TM. Following our strategy to spanthe tangent space by the Takagi vectors and real zero-modes, we readily have a parametrization at hand. Ac-cording to Eq. (18) we can express each configurationon the tangent space by real coordinates b k , specifying avector in the Lie-algebra. Note that for the zero-modes( α = 0) we have no complex prefactor. Using these co-ordinates, one can think of applying various different up-date procedures on the tangent space, starting from asingle random walk Monte Carlo (crude Monte Carlo),over a Hybrid Monte Carlo to a trained flow-based neu-ral network [45, 46]. The important steps are:1. Choose a critical configuration, to specify the tan-gent space on which to carry out the updates. Thisconfiguration may also serve as starting configura-tion.2. Calculate the real Hessian M and determine itseigenpairs ( α i , v i ). This has to be done only oncefor a given critical configuration, independent ofthe coupling β .3. Propose a new configuration by drawing a set ofreal coordinates b k . The proposed configurationis than specified according to Eq. (18). Performan accept/reject step based on the real part of theaction difference between the old and new config-uration. If the proposed configuration is outsidethe boundaries of the TM, we assign to it a zeroprobability and reject the proposed configuration,recording the old configuration like in the Metropo-lis algorithm.4. Finally take the remaining sign problem into ac-count by reweighting with the imaginary part ofthe action. G. Leading and subleading thimbles
So far, we have restricted the sampling to a singleTM attached to a specific critical configuration. Ulti-mately, we are interested in sampling a compact man-ifold that is homotopic to the original integration do-main. In the standard thimble decomposition the com-bination of various Lefschetz thimble leads, by defini-tion, to a multi-modal probability distribution. Sam-pling such distributions by a Monte Carlo procedure isdifficult. As a solution, a tempered sampling procedurewas proposed [33, 34]. Another possibility are indepen-dent Monte Carlo processes on each thimble. However,in this case the relative weights between thimbles needto be known. One possibility to infer these values is byusing prior knowledge of a physical observable for nor-malization [35].Here we construct a homotopic manifold by piecewisedefinition, where we use the TMs as building blocks. Asour construction deviates from the thimble decomposi-tion especially close to the boundaries, we do not haveinfinite action barriers between the patches. Hence, asampling procedure that proposes configurations acrossboundaries would be in principle possible. However, wefound that it is most convenient to sample the main tan-gent space with a single Monte Carlo chain and take allremaining patches into account via reweighting. We ex-emplify this procedure for a system of two TMs, τ , τ .Calculating expectation values over this extended regionrequires the relative weight Z /Z . Here Z denotes thepartition function corresponding to the subleading tan-gential manifold and Z refers to the corresponding quan-tity on the main tangential manifold. It holds (cid:104)O(cid:105) τ ∪ τ = (cid:82) τ d U O [ U ] e − S [ U ] + (cid:82) τ d U O [ U ] e − S [ U ] (cid:82) τ d U e − S [ U ] + (cid:82) τ d U e − S [ U ] = (cid:104)O(cid:105) τ + ( Z /Z ) (cid:104)O(cid:105) τ Z /Z ) . (20)Following the method proposed in [36], we introduce amapping f : τ −→ τ . (21)It maps configurations from one of the two patches tothe other. With this mapping we can express the ratio Z /Z as Z Z = (cid:82) τ d U e − S [ f ( U )]+ S [ U ] det[d f ] e − S [ U ] (cid:82) τ d U e − S [ U ] = (cid:10) e − S ◦ f + S det[d f ] (cid:11) τ . (22)It remains to find a suitable f . Since we consider onlytangent spaces, f is linear and can be constructed as abasis transformation from the Takagi basis of τ to τ .The Jacobian det[d f ] is therefore a constant factor. As the eigenbases of the real Hessian M are orthogonal andcan be chosen to have determinant one, the Jacobian de-pends consequently only on the the different sets of com-plex prefactors c k of the Takagi vectors. These dependin turn only on the sign of the eigenvalues α k . Thereforewe find for the Jacobiandet[d f ] = (cid:16)(cid:113) − β ∗ | β | (cid:17) n ( − ) τ (cid:16)(cid:113) + β ∗ | β | (cid:17) n (+) τ − n (+) τ , (23)where n (+) and n ( − ) are the number of positive/negativeeigenvalues of the real Hessian M (see Sec. II D) at therespective patch. Here we assumed that τ is the patchattached to the main critical manifold having only posi-tive eigenvalues being the global minimum. If the patches τ and τ are of different size, one may introduce an inde-pendent scale parameter to the mapping f , which is thanalso reflected as factor in the Jacobian. As the main tan-gential manifold is usually the largest such a factor is notnecessarily needed in practice. It may however be usedfor optimization purposes. H. Alternative updates on the main tangent space
As outlined in Sec. II D we can sample the main tan-gent space by just setting each complex factor c k to (cid:112) β ∗ / | β | . In practice this tilts every link, since in thiscase all eigenvalues are positive or zero. Having identicalcomplex prefactors, we can perform a basis transforma-tion of our coordinates to the unit basis. Therein each co-ordinate corresponds to a single gauge degree of freedom.This allows for applying a local heat bath or Metropolisalgorithm. The situation discussed here is shown at thebottom of Fig. 2. The unbounded critical manifold inimaginary direction has to be dealt with.We limit the plaquette values in accordance with theintersection points of the tangent spaces in the oneplaquette-model, creating TMs (see Fig. 3). Similarlyas before, we define the region outside the boundary tohave zero probability.But since link variables can still diverge in imaginarydirection we have to apply a second limit or just recordvariables, which are unaffected by the zero modes. Inour theory, these are the plaquettes variables. An alter-native is to adopt gauge cooling [47] to make this problemmilder. A severe limitation of this method stems from thefact that only observables invariant under all zero-modescan be measured. Taking pure Yang-Mills theory it isnot possible to measure the Polyakov loop. It is invari-ant under gauge transformations but not under a globalzero mode. The latter is represented by changing all linksin the same direction and amount leaving the plaquettesand therefore the action invariant. III. APPLICATION TO A 2-DIMENSIONALU (1) -GAUGE THEORY
In this section we apply the above outlined method totwo-dimensional pure U(1) gauge theory with periodicboundary conditions. Recently, complementary studieson this theory with a sign problem have appeared in theliterature. In [48] the complex Langevin method is em-ployed. The complex action here is caused by a non-zerovacuum angle. On the other hand in [49] the theory witha complex gauge coupling is investigated by means of thepath-optimization method.
A. The effective degrees of freedom
We formulate the theory in terms of its effective de-grees of freedom.Eq. (8) can be verified easily by noting that every link ap-pears twice in all plaquettes. Hence the links cancel eachother when being multiplied altogether. Consequently,one plaquette can be expressed in terms of all others.The action of the two-dimendional theory is rewritten asfollows S = − β (cid:88) ( x,t ) (cid:54) =(0 , (cid:0) P , ( x, t ) + P − , ( x, t ) (cid:1) (24)+ (cid:89) ( x,t ) (cid:54) =(0 , P − , ( x, t ) + (cid:89) ( x,t ) (cid:54) =(0 , P , ( x, t ) , neglecting constant terms. The last term is called toron term. One can reformulate the full theory in termsof these plaquettes variables having a reduced partitionsum, which gives the same expectation values for observ-ables that are invariant under zero modes Z = (cid:90) (cid:89) ( x,t ) (cid:54) =(0 , d θ ( x, t ) × exp β/ (cid:88) ( x,t ) (cid:54) =(0 , (cid:16) e iθ ( x,t ) + e − iθ ( x,t ) (cid:17) (25) × exp (cid:104) β/ (cid:16) e − i (cid:80) ( x,t ) (cid:54) =(0 , θ ( x,t ) + e i (cid:80) ( x,t ) (cid:54) =(0 , θ ( x,t ) (cid:17)(cid:105) . For a full derivation see App. B. The periodic boundaryconditions are represented by the toron term we placedat position (0 , -0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 0.5 1 1.5 2 2.5 3 3.5 4 R e < / ( P + P - ) > Re( β ) Approx ord 0Approx ord 2Approx ord 4Approx ord 6Approx ord 8BesselLOBesselNLORew(high stat) FIG. 4. Thimble hierarchy depending on Re ( β ) in the ap-proximation on a 4 × β ) = 1. the integral factorizes in plaquettes variables yielding Z = (cid:34)(cid:90) U(1) d P e β/ P + P − ) (cid:35) V = [ I ( β )] V . (26)The plaquette expectation value is then <
12 ( P + P − ) > = I ( β ) I ( β ) , (27)which is exactly the the same as for the one-plaquettemodel. There is no volume dependence. Since the dif-ference between periodic and open boundary conditionsvanishes in the infinite volume limit, this is the expectedvalue.We construct the approximation scheme by successivelyincluding TMs as integration domains. Thereto we splitthe integral according to the critical points P = ± Z = (cid:20)(cid:90) τ d P e β/ P + P − ) + (cid:90) τ d P e β/ P + P − ) (cid:21) V =: [ Z + Z ] V = V (cid:88) k =0 (cid:18) Vk (cid:19) Z V − k Z k . (28)We can map this to the lattice by k denoting the num-ber of plaquettes being − (cid:18) Vk (cid:19) is the number of such combinations on the lattice.Eq. (28) allows to calculate approximate values for thecomparison with numerical simulations (see Fig. 4).Complementary, there exists a formal solution for thefor the full lattice theory with periodic boundary con-ditions involving no approximations. We can write thepartition sum as Z = (cid:90) d U exp ( − S [ U ]) = + ∞ (cid:88) n = −∞ [ I n ( β )] V , (29) ord1ord2ord3ord4 - - - - ( β ) r e l a t i v e w e i gh t Z R , j / Z R , FIG. 5. Comparison of the ratio of the partition functions Z R,j /Z R, = (cid:104) exp( − ( S R,j − S R, )) (cid:105) i.e. the relative weightsbetween the lowest four subleading and the leading order TMs(basis configurations). The index j labels the subleading or-ders and the subscript 0 indicates that the expectation valueis calculated on the configurations measured on the leadingorder tangential surface, as in [36]. Here, the lattice is of size4 ×
4, Im ( β ) = 1 and spherical boundaries are used. Thedata reflects very well the hierarchy described in Sec. III B.As a comparison the functions exp( − S R [ P crit ]) are added tothe plot in light colors. being a series in modified Bessel functions I n ( β ), where V is the number of plaquettes [50, 51]. The leading orderof this series corresponds to our approximation. The fol-lowing orders take finite volume effects into account andyield the exact result provided that the series convergesfor the given value of β . B. Critical manifolds and their hierarchy
For an even number k in Eq. (28) the critical config-urations in the approximation and in the original latticetheory coincide. In contrast, for odd k , this does nothold due to the periodic boundary conditions. However,it is possible to construct critical manifolds being (arbi-trarily) close to the corresponding configurations in theapproximation. For odd k < V / − k plaque-ttes to e i ( π − (cid:15) ) and the rest to e i(cid:15) such that Eq. (13) issatisfied. This leads to the choice (cid:15) = πV − k , (30)which we refer to together with the critical configurationsfor even k as basis configurations from now on. Further-more as discussed in Sec. II C, Eq. (13) gives rise to amultiplet of configurations being related to the basis con-figurations by a symmetry. To see this, note that these (cid:15) configuration have topological charge (cid:100) k (cid:101) . Now, wecan change the topological sector by adding/ subtracting πV − k from (cid:15) . We observe that if | (cid:15) | < π/
2, the Takagivectors do not change. Consequently we can apply thistransformation directly to measured plaquettes values bymultiplying them with e ± π/ ( V − k ) , where the sign is cho- sen whether we have a e i ( π − (cid:15) ) or e i(cid:15) plaquettes. Thistransformation does not only apply to odd k but also toeven k . By this procedure, we reach 2( V / − (cid:100) k (cid:101) ) criticalmanifolds in different topological sectors for odd k and2( V / − − k ) + 1 for even k , respectively.Second, we can go from the k -th TM to the ( V − k )-th TMby changing the plaquettes from the critical configurationaccording to P crit1 , ( x ) −→ − ( P crit1 , ( x )) − . (31)The corresponding real Hessian M → ( − M ) changes signand has the same eigenvectors and zero modes. Onlythe non-zero eigenvalues α → ( − α ) change sign. Conse-quently and using the fact that eigenvectors are uniqueup to a non-zero scalar multiplication we get the Tak-agi vectors for the ( V − k )-th TM by multiplying theTakagi vectors from the k -th TM by √− i . For k (cid:54) = V / , V / ±
1, all zero modes do not change the pla-quettes. Writing a plaquette configuration on the k -thtangent space as P crit1 , ( x ) e i ∆ ϕ ( x ) , we can write the map-ping to the configuration on the opposite ( V − k )-th tan-gent space as P crit1 , ( x ) e i ∆ ϕ ( x ) → − ( P crit1 , ( x )) − e − ∆ ϕ ( x ) , (32)having again a transformation for directly measuring theplaquette values. Having these, the same procedure toreach the different topological sectors can be applied tothese plaquette values.An exception is the case k = V /
2, which has an addi-tional zero mode, which can be parametrized by P ( x ) = e iϕ and P ( x ) = e i ( π − ϕ ) , (33)for V / ± i , which would be assigned to k = V / ± (cid:18) VV / (cid:19) .The critical manifolds form a hierarchy depending ontheir associated value of the action S = 2 kβ, k even , ≈ β (cid:16) k + π V − k ) (cid:17) , k odd . for the basis configurations. Depending on Re ( β ) thecritical manifolds differ in importance for the partitionsum since on thimbles and suitably bounded TMs thereal part of the action is minimal at the critical manifold.Consequently, if β is purely imaginary, every thimble orTM contributes equally. Otherwise one can obtain anapproximate result by taking only a few thimbles or TMsinto account as the others are exponentially suppressed.Fig. 5 illustrates the hierarchy of the critical manifolds0 R e ( S ) R e ( S ) -Plaquette( )-Plaquette FIG. 6. Real part of the action of a single plaquette on thetangent space of the thimble. At the top, we look at thesituation φ = 0, i.e. P = 1 and P = − (cid:15) = π/ considering subleading TMs up to order k = 4 restrictedto the basis configurations. The simulation setup usedfor the shown data is described in detail in Sec. III D andIII F. C. Ensuring homotopy
To ensure global homotopy, we need to make sure, thatthe TMs form a patchwork covering [U(1)] V . Therefore,we have to create boundaries, which match each otherexactly. Strictly speaking, this is not generally possible inhigher dimensions, since there is no theorem that tells us,that these tangent spaces have to intersect in this manneras it is the case for thimbles. But we can at least getclose to something alike minimizing the systematic errorintroduced by homotopy violations as much as possible.Our approach comes with thinking in effective degrees offreedom being plaquettes variables with a toron term asshown in Sec. III A. Looking at the different (subleading)tangent spaces, we have applied several schemes, eachbased on different criteria: 1. Real plaquette boundaries based on the one plaque-tte model: These are shown in Fig. 3. The plaque-tte variables take values in a region bounded by theintersection points of the two tangents in the one-plaquette model. In the full lattice theory, we stillfind these intersection points for the transitions inconfigurations where k is even. We therefore limitthe real parts of the plaquettes by these intersec-tions. For (cid:15) (cid:54) = 0 configurations the transition looksdifferent and we take it into account by a shift ofthe boundaries. Having transition points does notmean that we get full homotopy. We would needto identify intersections with real dimension 2 V − S R is a simple rising quadraticfunction [52, Lemma 2.2]. Since the TM is close toit in the vicinity of the critical manifold, we observethe similar rising behaviour up to some distanceillustrated for the local action in Fig. 6. We limitthe plaquettes variables by the local maxima of S R .Since configurations are exponentially suppressedby it and the fact that this distance goes beyondthe intersection points mentioned before allows alarger part of configuration space to be explored,while the systematic error stays small.4. Spherical boundaries: This is the most conservativechoice being independent from the coordinate sys-tem. We observe that wandering along the Takagivector with the highest eigenvalue α = 8 (this istrue for all even lattice sizes including and greaterthan 2 ×
2) for the main TM turns it into the ulti-mate subleading TM (whose critical configurationhas all plaquettes equal to −
1) intersecting withits counterpart ( α = −
8) from there. This can beunderstood in terms of the effective d.o.f. as a diag-onal connection in a hypercube. The intersectionmarks a corner in the cube containing the main1 R e < / ( P + P - ) > Re( β ) Approx ord 0BesselLOBesselNLOMain TM | < e - i S I > | Re( β ) Main TMReweighting FIG. 7. Result of the simulation on the main TM with real plaquette boundaries for constant Im ( β ) = 1 on a 4 × β ), the result approximates the BesselNLO result. This indicates that we simulate the actuallattice theory and not by accident a cartesian product of the one-plaquette model. On the right is the average sign comparedwith standard reweighting. TM. We now choose the radius of the inner sphereof the cube to make sure, we do not intersect withother TMs, where we can assign a radius in thesame way. We get r max = √ V π √ V − (cid:115) β | β | , (34)as the radius of the effective sphere (the distancefrom the critical manifold is calculated in terms ofthe non-zero eigenmodes).A disadvantage is that here the curse of dimension-ality hits directly into the calculations, since withrising dimension, the volume of the inner sphere be-comes negligible in comparison to the cube and weexplore only small portions of configuration space.One can counter that by scaling the radii payingthe price of overlapping TMs introducing anothersystematic error.5. Im(S) boundaries: The idea here is that the TMsgo into the other TMs by their intersections. Sothe imaginary part of the action of one TM has tochange continuously to the one of the other (Onthimbles Im S is constant and changes abruptly attheir singular intersection, which in our case is atinfinity.). Having a pronounced hierarchy allows usto set the boundaries using this feature by e.g. al-lowing one tangent space only to vary in an intervalof Im S and the other one on a subsequent interval.Problematic is the fact that we have intersectionsof one TM with multiple TMs in different orders.So here we possibly limit the explored space toomuch.All in all, we have found a combination of real plaque-ttes boudaries with (cid:15) shifts and action boundaries most promising. To that end we calculate the former and cor-rect them, if they extend futher than the action bound-aries, which is especially important for TMs with high (cid:15) . Otherwise we can fall into the regions with negativeaction, where the the TM is far away from the thimble,see Fig. 6. D. Algorithm for sampling on the main TM
In the following we explain the sampling method togenerate configurations on the main tangential surface.(1) Diagonalize the Hessian at the main critical pointwhere all plaquettes are +1.(2) Construct the parametrization of the main TMsurface. The eigenvectors from (1) with non-zeroeigenvalues correspond to thimble directions. Tiltthose as described in Eq. (16) above to obtain theTakagi basis { ξ i i = 1 , . . . , V } .(3) Define boundaries. Before running the simulationspecify the configuration space to be sampled bychoosing suitable boundaries of the main TM. Forpossible choices see Sec. III C.(4) Run the Monte Carlo simulation. A vector on theleading TM in the Takagi basis is given by ξ = φ i ξ i with φ i ∈ R . Beginning with a cold start at thecritical manifold, i.e. φ i = 0 ∀ i = 1 , . . . , V samplethe φ i via the Metropolis algorithm using propos-als ∆ φ i ∼ N (0 , σ ). Here, a sweep is defined by ap-plying a Metropolis accept-reject step for every di-rection i . The Metropolis updates are constructedsuch that the following conditions are satisfied(i) Configurations with S R < | < e - i S I > | Volume Main TMReweighting
FIG. 8. The average sign plotted for increasing volume atconstant β = 2 + 1 . i . (ii) A proposed configuration outside of the spec-ified boundaries is being rejected.The measured configurations are stored to be used forthe reweighting on the subleading TMs. This part of thealgorithm is described in detail in Sec. III F.Like already mentioned in Sec. II H, there is also a localupdate algorithm (see App. C). E. Numerical results on the main TM
Looking at the approximation in Fig. 4, we expect themain tangent results to roughly follow the zero order ap-proximation. For high β R , the other tangent spaces areexponentially suppressed allowing for convergence to thefull result. Noticing that the full result for this rangeis slightly above the one plaquette model due to finitevolume effects, we hope to see the same from the simula-tion, which is the case, see Fig. 7. This proves that we donot accidentially just simulate the one-plaquette modelby our procedure.Another important point, is that the sign problemshould be reduced in comparison with standard reweight-ing, which is also the case (see second plot in Fig. 7).So, we can expect for high enough β R (e.g. here largerthan 3) correct results for simulations at complex betawith a lesser sign problem. To extend this range, we needto take the subleading TMs into account, see Sec. III F.For standard phase reweighting, the average signshould exponentially decrease with increasing space-timevolume. Since our simulation works similarly just on atilted space, we see the same happening here. The dif-ference is that our average sign is higher than the onefor standard reweighting and the slope is less steep, seeFig. 8. Therefore, higher lattice volumes are more eas-ily accessible in our approach. Indeed, we needed to in-crease statistics for the reweighting simulations, while thenumber of samples for the different volumes in the TMsimulations remained the same. | < σ / | σ | > | Re( β )Takagi ord 0Takagi ord 10Takagi ord 20Takagi ord 32Reweighting FIG. 9. Absolute values of the expectation values of the phaseof the reweighting factor of Eq. (35) depending on the num-ber of included tangent spaces and β in comparison with astandard reweighting simulation. F. Reweighting of subleading TMs
For the reweighting, the observable from Eq. (20) be-comes < O > = (35) < e − iS I O + (cid:80) nk =1 m k det[d f k ] e S R − S ◦ f k O ◦ f k > τ < e − iS I + (cid:80) nk =1 m k det[d f k ] e S R − S ◦ f k > τ , where n denotes the number of subleading tangent spaceorders, one wants to incorporate, m k are the multiplici-ties and f k are linear transformations projecting the mainTM τ onto the subleading tangent space τ k . This is doneby aligning the Takagi basis using the real vectors v j tocalculate a transformation matrix γ ( k ) ij = ( v ( k ) j ) T v (0) i , (36)which we apply to the subleading Takagi basis z ( k ) i = (cid:88) j γ ( k ) ij c ( k ) j v ( k ) j , (37)where the c ( k ) j are the complex prefactors from Eq. (16).The z ( k ) i are now our new aligned basis vectors for τ k .Therefore we can directly project from the leading tosubleading tangent spaces. As already mentioned inSec. II G, depending on the other boundaries of the sub-leading tangent space, one can apply scaling factors tothe variables. We observed, that for β R , β I ≥
0, themain TM is always larger (or equally large) as the othersubleading TMs. Instead of rescaling, we use an indica-tor function χ τ k ( U ) for the boundaries: If the projected3 R e < / ( P + P - ) > Re( β ) BesselLOBesselNLOTakagi ord0Takagi ord5Takagi ord10Takagi ord15Takagi ord20Takagi ord25Takagi ord30Takagi ord31Takagi ord32 I m < / ( P + P - ) > Re( β ) BesselLOBesselNLOTakagi ord0Takagi ord5Takagi ord10Takagi ord15Takagi ord20Takagi ord25Takagi ord30Takagi ord31Takagi ord32 FIG. 10. Takagi simulation incorporating different orders of reweighted tangent spaces on an 8 × β ) = 1.The dashed line indicates the Bessel result, see Eq. (29). Each order k contains for practical reasons the V − k order and alltopological TMs associated with these. Shown are the real part and imaginary part of the plaquette expectation value. space f k ( τ ) = ˜ τ k ⊇ τ k , then we have (cid:90) τ k d U g ( U ) = (cid:90) ˜ τ k d U g ( U ) χ τ k ( U ) (38)with g ( U ) an arbitrary function on the subspace. So,we simply set the integrand to zero in the reweightingprocess, if the projected configuration is out of bounds.We use several symmetries to incorporate the differenttopological sectors for each order as well as a mappingfrom the k -th to the ( V − k )-th TM already discussed inSec. III B. In the end, we have to diagonalize only V /
V /
V / − V / ϕ = π of the additional zero mode (33).During the reweighting process, we need to check for ev-ery contribution, if it is in its predefined boundaries, sincethey also differ over topological sectors (see Sec. III C).We applied the procedure on an 8 × β ) = 1 and compared different orders of reweightingwith the BesselNLO result (see end of Sec. III A). Theresults are shown in Fig. 10 as well as the average sign inFig. 9 depending on how many orders of TMs one takesinto account. The boundaries chosen are real plaquetteboundaries based on the one plaquette model, since theyprevent overlapping of the TMs, while allowing a largespace to be explored. IV. DISCUSSION
As we have stated already in the introduction, the par-ticular choice of our deformation has two important ad-vantages:( A
1) Due to the flatness of the patches a parametriza-tion in terms of real coordinates and basis vectorsis easily constructed. It is thus straight forward to realize a sampling procedure on a particular patch.In contrast to the generalize Lefschetz thimble ap-proach, there is no need to solve a flow equationto propose a new configuration nor to evaluate aJacobian at each sampling point.( A
2) Since our coordinate systems originate from a crit-ical configuration on each patch, which is the con-figuration that receives the largest weight on thatpatch, a reweighting from one patch to another ispossible without any severe overlap problem.In the case of a pronounced hierarchy among criticalconfigurations, we have demonstrated in the case of the2d U(1) theory, that a well controlled approximationscheme emerges if successive contributions from sup-pressed patches are taken into account.On the other hand there are two – as we believe lesssevere – disadvantages:( D
1) Sampling on the flat tangential manifold ratherthan on the curved thimble, does not erase the signproblem completely. It would be interesting to an-alyze whether the optimal deformation of the orig-inal manifold, which minimizes the combined signproblem of the action on the integration domainand the one introduced by the Jacobian, is closerto the thimble decomposition or our decompositionfrom flat patches. We leave this for future investi-gations. We have demonstrated that the resultingsign problem in the case of the 2d U(1) gauge the-ory is, although sill exponential in volume, verymild compared to the standard reweighting proce-dure.( D
2) In order to construct our deformation, we need toknow all relevant critical configurations. Depend-ing on the space-time dimension, volume and gaugegroup, this can be a very large amount of criticalconfigurations.4It is one of the main results of this work that we haveidentified a large number of critical configurations, whichare characterized by specific distributions of plaquettesacross on the lattice. In particular, we have demonstratedthat we can chose these critical configuration from themaximal torus of the original gauge group. We haveshown further that there exists a large number of de-generate critical configurations and patches due to lat-tice symmetries. For the reweighting procedure we thusneed to take only one of these patches into account ifappropriate combinatorial multiplicity factors are used.The results look promising and systematic errors by non-homotopy vanish with growing lattice size.We want to emphasize here that the choice of a com-plex coupling β is not at all a pathological choice. Inthe limit of a purely imaginary coupling, the kernel ofthe discussed partition sum can be seen as the real timepropagator of the theory. For the evaluation of real timecorrelation functions in a thermal bath, the Schwinger-Keldysh formalism is usually applied. Our samplingstrategy might also be applied to the Schwinger-Keldyshcontour, even though in this case an additional sign prob-lem arises from the edges of the contour.Moreover, the critical configurations we have identifiedhere are not only critical in the case of the Yang Millsaction with complex coupling β . The same configurationsremain critical when we introduce fermionic matter fieldswith a chemical potential µ . In this case the effectiveaction might be written as S eff ( µ ; U ) = βS G ( U ) − Tr ln D ( µ ; U ) . (39)Hence, the action gradient is the sum of a contribu-tion from the gauge ( S G ( U )) and the fermionic part( S F ( µ ; U )) of the action. That the action gradient ofthe gauge part vanishes at our critical configurations hasbeen discussed in detail above. The fermionic contribu-tion to the gradient is given as ∂ x,ν,a S F = i Tr (cid:34) D − (cid:16) e µδ ν, U ν ( x ) − e − µδ ν, U † ν ( x + ˆ ν ) (cid:17) t a (cid:35) . (40)For those critical points that are not only chosen fromthe maximal torus of the gauge group but are also con-structed from center elements of the original gauge group,the fermionic contribution vanishes as well. As all linkvariables are proportional to the unit matrix, the (nonesparse) inverse of the fermion matrix D − contains N × N diagonal blocks which are also proportional to the unitmatrix. We conclude that the matrix multiplying thegenerator t a is proportional to the unit matrix and assuch the whole expression vanishes. We thus hope thatour strategy might also prove useful in this case, i.e. thefield theoretical description of dense matter, includingQCD at net-baryon number density. V. CONCLUSION
We have put forward here a novel nonperturbative lat-tice approach for Yang Mills theories with a complexgauge coupling β . The approach is based on a defor-mation of the original integration domain of the theoryinto complex space. Guided by the Lefschetz thimble de-composition of the partition sum we have chosen a newintegration domain. This is constructed piecewise frompatches of tangential manifolds to the relevant Lefschetzthimbles.For the numerical implementation we have chosen toset up a Monte Carlo procedure on the main tangentspace only. All further contributions from other patchesare taken into account by reweighting. We have testedour approach by applying it to the case of 2d U(1) gaugetheory. Here, it far out-performs our benchmark simula-tion with standard reweighting, see Fig. 1.Based on this observation we plan to apply our ap-proach to general U( N ) and SU( N ) gauge groups in 4d.While our approach has shown to feature an exponen-tially better performance than standard reweighting, itremains to be shown that the sign problem stays numeri-cally manageable also in these cases. We hope to addressthese questions in a forthcoming publication. Moreover,we envisage simulations of fermionic matter fields at fi-nite chemical potential as well as real-time lattice theo-ries with expectation values along the Schwinger-Keldyshcontour. ACKNOWLEDGEMENTS
The authors thank Andrei Alexandru, Benjamin J¨ager,Alexander Lindemeier and Ion-Olimpiu Stamatescu fordiscussions.C. Schmidt and F. Ziesch´e acknowledge support byDeutsche Forschungsgemeinschaft (DFG, German Re-search Foundation) through the Collaborative ResearchCentre CRC-TR 211 ’Strong-interaction matter underextreme conditions’ project number 315477589 and fromthe European Union’s Horizon 2020 research and in-novation program under the Marie Sk(cid:32)lodowska-Curiegrant agreement No H2020-MSCAITN-2018-813942 (Eu-roPLEx). This work is further supported by the ExteMeMatter Institute EMMI, the Bundesministerium f¨ur Bil-dung und Forschung (BMBF, German Federal Ministryof Education and Research) under grant 05P18VHFCAand by the DFG through the Collaborative ResearchCentre CRC 1225 (ISOQUANT) as well as by DFG underGermany’s Excellence Strategy EXC-2181/1-390900948(the Heidelberg Excellence Cluster STRUCTURES). M.Scherzer acknowledges support from DFG under grantSTA 283/16-2. F. P. G. Ziegler acknowledges supportby Heidelberg University where a part of this work wascarried out.5
Appendix A: General complex Hessians
If we can not write the Hessian as H = βM , with a realmatrix M as in Sec. II D, we consider real and imaginaryparts of the Takagi Eq. (14) separately (see e.g. [38])( H R − iH I )( v R − iv I ) = λ ( v R + iv I ) ⇔ (cid:18) H R − H I − H I − H R (cid:19) (cid:18) v R v I (cid:19) = λ (cid:18) v R v I (cid:19) . (A1)The matrix is by definition symmetric and we have onlyreal eigenvalues, positive λ for Takagi and negative λ for Anti-Takagi vectors. The zero modes span again thewhole critical manifold. Since for compact gauge groups,we are only interested in the real subspace, which arenaturally compact. Therefore, we impose v I = 0 and get H R v R = 0 and H I v R = 0 . (A2)This implies that the critical manifold is spanned by themutual zero modes of the real and imaginary part of theHessian. The strategy therefore is to calculate the zeromodes either parts and test, if these are also zero modesof the other part. The number of real zero modes and ofthe Takagis have to add up to the overall dimension. Appendix B: The toron formulation
We reformulate the theory as mentioned in Eq. (24).Then we gauge all time-like directions into one link. Thisreduces the degrees of freedom by N x ( N t − θ and expressthe links (in their algebra representation φ µ ( x, t )) interms of these. Thereto we write direction vectors( (cid:126)ϕ ( x, t ) , (cid:126)ϕ ( x ) , (cid:126)ϕ (0)) in link space indicating how thenew variables change the links in their respective direc-tion. This yields the following parametrization φ µ ( y, t ) = (cid:88) ( x,τ ) (cid:54) =(0 , θ ( x, τ )( (cid:126)ϕ ( x, τ )) φ µ ( y,t ) (B1)+ (cid:88) x θ ( x )( (cid:126)ϕ ( x )) φ µ ( y,t ) + θ (0)( (cid:126)ϕ (0)) φ µ ( y,t ) for each remaining link. The θ ( x, t ) shall denote plaque-tte variables, while the θ ( x ) and θ (0) denote zero modesat space slices or the zero time slice. We will integratethem out later. All space-like links can be replaced byplaquette variables (cid:126)ϕ ( x, t ) = ˆ φ ( x, t mod N t ) − ˆ φ ( x, t − N t ) ,t ∈ { , . . . , N t } , (B2) and a zero mode (cid:126)ϕ ( x ) = N t − (cid:88) t =0 ˆ φ ( x, t ) , for each space slice. We use the notation ˆ φ to denote aunit vector in link space corresponding to the variable φ .We have V variables, which are linear independent,since we can express each space-like link byˆ φ ( x, t ) = 1 N t (cid:16) (cid:126)ϕ ( x ) + ( t −
1) mod N t (cid:88) k =1 k (cid:126)ϕ ( x, k ) (B3) − ( − t ) mod N t (cid:88) l =1 l (cid:126)ϕ ( x, N t − l ) (cid:17) . We replace the remaining N x time-like links with plaque-tte variables (cid:126)ϕ ( x,
0) = ˆ φ ( x,
0) + ˆ φ ( x,
1) (B4) − ˆ φ ( x + 1 mod N x , − ˆ φ ( x, , for x ∈ { , . . . , N x − } and the zero mode (cid:126)ϕ (0) = N x − (cid:88) x =0 ˆ φ ( x, . Using the fact that we already have a basis transform forthe space-like links, we can express the remaining time-like links in terms of these variables in a similar fashion toEq. (B3). This proves that our variable transformationis invertible and therefore we have a non-zero Jacobian.By the parametrization (B1) and the toron action (24),we rewrite the partition sum to Z = (cid:90) d θ (0) (cid:89) x d θ ( x ) (cid:89) ( x,t ) (cid:54) =(0 , d θ ( x, t ) det (cid:34) ∂ (cid:126)φ∂(cid:126)θ (cid:35) × exp β/ (cid:88) ( x,t ) (cid:54) =(0 , (cid:16) e iθ ( x,t ) + e − iθ ( x,t ) (cid:17) (B5) × exp (cid:104) β/ (cid:16) e − i (cid:80) ( x,t ) (cid:54) =(0 , θ ( x,t ) + e i (cid:80) ( x,t ) (cid:54) =(0 , θ ( x,t ) (cid:17)(cid:105) . Since our transformation is linear, the Jacobian det (cid:104) ∂(cid:126)φ∂(cid:126)θ (cid:105) is constant and drops out when taking expectation val-ues. The same is true for the integrals over the zeromodes. They do not change the plaquettes and the ac-tion. Their contribution is given by a constant factor.Therefore they can be dropped without changing expec-tation values. This reduces our degrees of freedom by N x + 1 to V −
1. The remaining effective degrees of free-dom θ ( x, t ) are called toron variables.6 Appendix C: A local update procedure for a2-dimensional
U(1) gauge theory
We parametrize this tangent space by arclength.Therefore, we have ϕ µ ( x ) = (cid:115) β ∗ | β | λ µ ( x ) , (C1)where λ µ ( x ) ∈ R . We limit this space by the intersectionsseen in the one plaquette model (see Sec. II D). Thereforewe haveΛ( x ) = λ ( x ) + λ ( x + ˆ0) − λ ( x + ˆ1) − λ ( x ) ∈ [ − πR, + πR ] , with R = Re (cid:115) β ∗ | β | . (C2)We enforce this limit by setting the probability forproposed configurations outside this limits to zero, ef-fectively rejecting them in the Metropolis step. The al-gorithm is according to that1. Go through the lattice by a checker board patternand pick accordingly λ µ ( x ).2. Make a proposal ∆ λ µ ( x ) ∼ Gauss( µ = 0 , σ ). andcalculate the two adjacent plaquettes Λ( x ), Λ( x − ˆ ν )with ν (cid:54) = µ .3. If Λ( x ) and/or Λ( x − ˆ ν ) is not within the bound-aries, set e − S (cid:48) = 0. Otherwise calculate the changeof the action ∆ S ( P ( x ) , P ( x − ˆ ν ))4. Do a Metropolis Accept/Reject step and beginagain from 1.Since our links can wander off into the imaginary di-rection, we perform gauge cooling steps, see i.e. [47].Locally for a site x we can analytically compute an op-timum for the gauge transformation V ( x ). Therefore welook at the unitarity norm F [ U µ ( x )] = (C3) (cid:88) x,µ Tr (cid:2) U † µ ( x ) U µ ( x ) + ( U † µ ( x )) − U − µ ( x ) − I (cid:3) , which simplifies in our case U µ ( x ) = e iφ µ ( x ) to F [ φ µ ( x )] = (cid:88) x,µ (cid:104) e − φ Iµ ( x ) + e φ Iµ ( x ) − (cid:105) = (cid:88) x,µ (cid:2) | U µ ( x ) | − + | U µ ( x ) | − (cid:3) , (C4)depending only on the absolute value of the links. Howa gauge transformation changes (C3) depends thereforeonly on its absolute value. For a local V ( x ) we get the local minimum at | V ( x ) | = (cid:115) (cid:80) µ ( | U µ ( x ) | − + | U µ ( x − ˆ µ ) | ) (cid:80) µ ( | U µ ( x ) | + | U µ ( x − ˆ µ ) | − ) , (C5)which we apply between sweeps.Alternatively one doesn’t need to record the link valuesand we restrict ourselves to just recording the plaquettevalues, which are naturally bounded. We apply the up-date steps only in terms of changes in links affecting theirneighboring plaquettes.Reweighting to different orders of TMs can also be ap-plied here by injectively mapping the plaquette configu-rations to link configurations and aligning the TMs basesto the unit basis, corresponding to the individual links.Given a set of plaquette values θ ( x, t ) in the algebra, apossible mapping to link variables φ µ ( x, t ) would be1. For x = 0 and all t ∈ , . . . , N t − φ (0 , t ) = θ (0 , t ) .
2. For x = 1, all t , set φ (2 , t ) = − θ (1 , t ) .
3. Successively for x ∈ , . . . , N x − t , set φ ( x + 1 , t ) = φ ( x, t ) − θ ( x, t ) .
4. The last space-slice x = N x − φ ( x,
0) = 0and for t ∈ , . . . , N t − φ ( x, t ) = φ ( x, t − θ ( x, t − φ (0 , t − − φ ( x, t − . The last plaquette value θ ( N x − , N t −
1) is automaticallytaken into account up to a factor of 2 π by the periodicboundary conditions. Now the values obtained for thealgebra of the links are used to map onto the subleadingTMs. Appendix D: Proof of the matrix identity
We consider the system of equationsTr[
M T a ] = 0 ∀ a, (D1)where M is a general N × N complex matrix and T a aregenerators of the lie algebra su ( N ) or u ( N ). We showthat Eq. (D1) holds, iff M = c for an arbitrary c ∈ C for su ( N ) and M = 0 for u ( N ).First note, that the T a also form the basis of the com-plex lie algebras sl ( N, C ), gl ( N, C ). This is done by al-lowing complex coefficients effectively removing the anti-hermiticity of the elements. Consequently elements areonly defined by being traceless for sl ( N, C ) and we have7 gl ( N, C ) (cid:39) C N × N .The statement does not depend on the basis of the Liealgebra: Let T a and S a be two basis of the same Liealgebra. We show thatTr[ M T a ] = 0 ∀ a ⇔ Tr[
M S a ] = 0 ∀ a (D2)by expanding S a = (cid:80) b c ab T b and using the linearity ofthe trace Tr[ M S a ] = (cid:88) b c ab Tr[
M T b ] = 0 . (D3)Especially the statement is equivalent toTr[ M T ] = 0 ∀ T ∈ sl ( N, C ) or ∀ T ∈ gl ( N, C ) . (D4)Note, that the backward direction is trivial, sinceTr[ c T ] = c Tr[ T ] = 0 for T being traceless and Tr[0 T ] =Tr[0] = 0 for T ∈ gl ( N, C ).For the forward direction we use the more general state- ment (D4). For T ∈ sl ( N, C ), suppose first M ij (cid:54) = 0for i (cid:54) = j . Then we can choose T lk = δ lj δ ki , which isobviously traceless and haveTr[ M T ] = (cid:88) k,l M kl T lk = (cid:88) k,l M kl δ lj δ ki = M ij (cid:54) = 0 . (D5)For M ii (cid:54) = M jj , i (cid:54) = j we take T lk = δ li δ ki − δ lj δ kj , whichis also traceless and haveTr[ M T ] = (cid:88) k,l M kl T lk = (cid:88) k,l M kl ( δ li δ ki − δ lj δ kj )= M ii − M jj (cid:54) = 0 . (D6)Since i and j were arbitrary M = c for some c ∈ C .For T ∈ gl ( N, C ) sup sl ( N, C ), T is not traceless anymoreand can also be e.g. , which excludes the M = c caseand M has to be zero to fulfill the statement. [1] I. M. Barbour, S. E. Morrison, E. G. Klepfish, J. B.Kogut, and M.-P. Lombardo, Nucl. Phys. B Proc. Suppl. , 220 (1998), arXiv:hep-lat/9705042.[2] Z. Fodor and S. Katz, Phys. Lett. B , 87 (2002),arXiv:hep-lat/0104001.[3] C. Allton, S. Ejiri, S. Hands, O. Kaczmarek, F. Karsch,E. Laermann, C. Schmidt, and L. Scorzato, Phys. Rev.D , 074507 (2002), arXiv:hep-lat/0204010.[4] C. Allton, S. Ejiri, S. Hands, O. Kaczmarek, F. Karsch,E. Laermann, and C. Schmidt, Phys. Rev. D , 014507(2003), arXiv:hep-lat/0305007.[5] R. V. Gavai and S. Gupta, Phys. Rev. D , 034506(2003), arXiv:hep-lat/0303013.[6] P. de Forcrand and O. Philipsen, Nucl. Phys. B , 290(2002), arXiv:hep-lat/0205016.[7] M. D’Elia and M.-P. Lombardo, Phys. Rev. D , 014505(2003), arXiv:hep-lat/0209146.[8] S. Kratochvila and P. de Forcrand, Prog. Theor. Phys.Suppl. , 330 (2004), arXiv:hep-lat/0309146.[9] A. Alexandru, M. Faber, I. Horvath, and K.-F. Liu,Phys. Rev. D , 114513 (2005), arXiv:hep-lat/0507020.[10] F. Karsch and K. Mutter, Nucl. Phys. B , 541 (1989).[11] P. de Forcrand, J. Langelage, O. Philipsen, andW. Unger, Phys. Rev. Lett. , 152002 (2014),arXiv:1406.4397 [hep-lat].[12] C. Gattringer and K. Langfeld, Int. J. Mod. Phys. A ,1643007 (2016), arXiv:1603.09517 [hep-lat].[13] G. Gagliardi and W. Unger, Phys. Rev. D , 034509(2020), arXiv:1911.08389 [hep-lat].[14] J. Ambjorn, K. Anagnostopoulos, J. Nishimura, andJ. Verbaarschot, JHEP , 062 (2002), arXiv:hep-lat/0208025.[15] Z. Fodor, S. D. Katz, and C. Schmidt, JHEP , 121(2007), arXiv:hep-lat/0701022.[16] K. Langfeld, B. Lucini, R. Pellegrini, and A. Rago, Eur.Phys. J. C , 306 (2016), arXiv:1509.08391 [hep-lat].[17] F. Karsch and H. Wyld, Phys. Rev. Lett. , 2242 (1985). [18] G. Aarts, E. Seiler, and I.-O. Stamatescu, Phys. Rev. D , 054508 (2010), arXiv:0912.3360 [hep-lat].[19] G. Aarts, F. A. James, J. M. Pawlowski, E. Seiler,D. Sexty, and I.-O. Stamatescu, JHEP , 073 (2013),arXiv:1212.5231 [hep-lat].[20] D. Sexty, Phys. Lett. B , 108 (2014), arXiv:1307.7748[hep-lat].[21] F. Attanasio, B. J¨ager, and F. P. Ziegler, Eur. Phys. J.A , 251 (2020), arXiv:2006.00476 [hep-lat].[22] F. Pham, Proc. Symp. pure Math , 319 (1983).[23] E. Witten, (2010), arXiv:1009.6032 [hep-th].[24] M. Cristoforetti, F. Di Renzo, and L. Scorzato(AuroraScience), Phys. Rev. D86 , 074506 (2012),arXiv:1205.3996 [hep-lat].[25] A. Alexandru, G. Basar, P. F. Bedaque, and N. C. War-rington, (2020), arXiv:2007.05436 [hep-lat].[26] M. Cristoforetti, F. Di Renzo, G. Eruzzi, A. Mukherjee,C. Schmidt, L. Scorzato, and C. Torrero, Phys. Rev. D , 114505 (2014), arXiv:1403.5637 [hep-lat].[27] A. Alexandru, G. Basar, P. F. Bedaque, and G. W. Ridg-way, Phys. Rev. D , 114501 (2017), arXiv:1704.06404[hep-lat].[28] F. Bursa and M. Kroyter, JHEP , 054 (2018),arXiv:1805.04941 [hep-lat].[29] Y. Mori, K. Kashiwa, and A. Ohnishi, Phys. Rev. D96 ,111501 (2017), arXiv:1705.05605 [hep-lat].[30] Y. Mori, K. Kashiwa, and A. Ohnishi, PTEP ,023B04 (2018), arXiv:1709.03208 [hep-lat].[31] A. Alexandru, P. F. Bedaque, H. Lamm, andS. Lawrence, Phys. Rev. D , 094505 (2017),arXiv:1709.01971 [hep-lat].[32] J.-L. Wynen, E. Berkowitz, S. Krieg, T. Luu, and J. Ost-meyer, (2020), arXiv:2006.11221 [cond-mat.str-el].[33] A. Alexandru, G. Basar, P. F. Bedaque, andN. C. Warrington, Phys. Rev. D , 034513 (2017),arXiv:1703.02414 [hep-lat]. [34] M. Fukuma, N. Matsumoto, and N. Umeda, Phys. Rev.D , 114510 (2019), arXiv:1906.04243 [cond-mat.str-el].[35] F. Di Renzo and G. Eruzzi, Phys. Rev. D , 014503(2018), arXiv:1709.10468 [hep-lat].[36] S. Bluecher, J. M. Pawlowski, M. Scherzer, M. Schlosser,I.-O. Stamatescu, S. Syrkowski, and F. P. Ziegler, Sci-Post Phys. , 044 (2018), arXiv:1803.08418 [hep-lat].[37] F. Di Renzo, S. Singh, and K. Zambello, (2020),arXiv:2008.01622 [hep-lat].[38] A. Alexandru, G. Basar, and P. Bedaque, Phys. Rev. D93 , 014504 (2016), arXiv:1510.03258 [hep-lat].[39] K. G. Wilson, Phys. Rev.
D10 , 2445 (1974).[40] M. Cristoforetti, F. Di Renzo, A. Mukherjee,and L. Scorzato, Phys. Rev.