Simulating Yang-Mills theories with a complex coupling

Jan M. Pawlowski, 2 Manuel Scherzer, Christian Schmidt, Felix P.G. Ziegler, and Felix Ziesché Institut für Theoretische Physik, Universität Heidelberg, Philosophenweg 16, D-69120 Heidelberg, Germany ExtreMe Matter Institute EMMI, GSI, Planckstraße 1, D-64291 Darmstadt, Germany Fakultät für Physik, Universität 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)


I. INTRODUCTION
The notorious sign problem hampers numerical simulations of many interesting physical systems, ranging from high energy to condensed matter systems. A sign problem is faced in numerical calculations of statistical models whenever the action becomes genuinely complex. Hence, standard Monte Carlo methods and in particular importance sampling drastically lose their efficiency with increasing lattice volume.
Examples of theories with a sign problem include real time calculations in lattice-regularized quantum field theories, i.e. lattice QCD in Minkowski space-time, lattice QCD with a nonzero vacuum angle θ, and lattice QCD with a nonzero baryon chemical potential µ B . For the latter two cases, many methods have been developed that potentially circumvent or solve this problem in the continuum limit. These methods include reweighting [1,2], Taylor expansions [3][4][5], analytic continuation from purely imaginary chemical potentials [6,7], canonical partition functions [8,9], strong coupling/dual methods [10][11][12][13], the density of states method [14][15][16], and complex Langevin dynamics [17][18][19][20], see [21] for a review of recent developments. However, all these methods so far face severe limitations that restrict their applicability in the continuum limit.
In the past decade deformations of the original integration manifold into the complex domain have been introduced, based on complex saddle points of the action, the Lefschetz thimbles [22,23]. If these deformations are chosen well, all physical expectation values obtained from an oscillatory integral remain unchanged but the sign problem is drastically alleviated. Thus, a numerical evaluation of the theory is accessible. By definition, the imaginary part of the action (phase of the probability density) is stationary on the thimble. A first Lefschetz thimble algorithm in the context of the QCD finite density sign problem was introduced in [24], for a recent review see [25]. Despite its great potential for beating the sign problem, simulations with Lefschetz thimbles have to overcome the following intricacies: (i) A parametrization of the thimble is a priori unknown and has to be obtained as the numerical solution of a flow-equation. The parametrization and the necessary Jacobian of the variable transformation are numerically demanding.
(ii) The curvature of the thimble manifold introduces a so-called residual sign problem through the Jacobian.
(iii) In most cases, one has to include relevant contributions from a large number of thimbles. Their relative weight may give rise to a further, residual, sign problem. In any case, the probability density becomes multi-modal.
These shortcomings have been addressed by formulating particularly efficient methods for the calculation of the Jacobian [26,27], by optimizing the deformation of the integration domain either based on a model ansatz [28,29] or by means of machine learning [30][31][32] and finally by using tempering [33,34] and other advanced strategies to foster transitions between thimbles and take into account contributions from multiple thimbles [35][36][37].
In the present work we analyze the structure of the Lefschetz thimble decomposition of pure gauge theories with gauge groups U(N ) and SU(N ) and complex coupling β. We propose to sample on the tangential manifold attached to the thimble of the main critical point, i.e. the critical point with the smallest action value. As expected, we find a full hierarchy of critical points (saddle points) that have to be considered depending on the coupling parameter β. We include successive subleading saddle point contributions via reweighting. The reweighting procedure is set up by using linear mappings from the main tangential manifold to the tangential manifold attached to the thimble of the target critical point.
Our approach has the following advantages over the flow-based generalized Lefschetz thimble approach [38]: (A1) During the sampling procedure, we neither have to flow our configuration nor have to calculate a Jacobian since it is constant on the tangential manifold.
(A2) As contributions from subleading thimbles are taken into account by reweighting [36], there is no ergodicity problem due to potential barriers between thimbles. The reweighting procedure does also not introduce any overlap problem since critical points are mapped onto critical points.
The disadvantages are: (D1) Sampling on the tangential manifold rather than on the thimble itself does not completely resolve the sign problem, however, there is no residual sign problem from the Jacobian [26].
(D2) The critical points need to be known.
The advantages are clearly demonstrated for the benchmark case of an U(1) gauge theory, as summarized in Fig. 1. There we compare our results for the real part of the plaquette with that produced by the standard reweighting procedure. For the same numerical costs the error bars of our results are reduced by orders of magnitude.
This will be discussed in more detail in this work, which is organized as follows: In Sec. II we derive necessary formulas for our setup. Emphasis is devoted to discussing the critical manifolds of the Yang-Mills action as well as the Lefschetz thimbles. Moreover we explain our update and reweighting procedures. In Sec. III we apply our algorithm to the case of a 2-dimensional U(1) gauge theory. In Sec. IV we discuss prospect of our approach. Finally we conclude in Sec. V.

A. Overview
We consider the standard discretization of the Yang-Mills action [39] The summation is done such that each plaquette is counted with only one orientation. The link variables U µ (x) are elements of the gauge group, which we consider to be a Lie group and in particular U(N ) or SU(N ).
A sign problem is introduced by choosing general complex couplings β.
We aim to update a given gauge field configuration such that the imaginary part of the action Im(S) varies slowly and hence the remaining sign problem is mild.
In order to achieve this goal the link variables are complexified, i.e. they are allowed to take values within the larger group GL(N, C) or SL(N, C), respectively. By generalized Picard-Lefschetz theory, there is a smooth middle dimensional manifold connected to each complex critical 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 already stated, updates that stay on these thimbles are computational very demanding and are invariably global updates [24,38,40,41]. To reduce the computational demand, it has been argued that one does not have to stay exactly on the thimble in order to reduce the sign problem significantly [29,42]. Any deformation of the original integration domain which has the correct asymptotic behavior will be sufficient. Here we construct an integral deformation, which is the union of all tangential manifolds to the critical points. As for the pure gauge theory all critical manifolds are known, this manifold is straightforward to parametrize. After discussing the critical manifolds in Sec. II B and II C we discuss properties of Lefschetz thimbles and tangential manifolds in Sec. II D and II 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 union of flow lines generated by the steepest descent equation which end in a non-degenerate critical point of the action. The degeneracy of critical points due to gauge symmetry necessitates the application of Witten's concept of generalized Lefschetz thimbles [43, 3.3].
The critical manifolds can be described in terms of plaquette variables. This is seen by examining the gradient of the action Here t a are the generator matrices of the related Lie-Algebras su(N ), u(N ). In our notation those are Hermitian, and for N > 1 they satisfy the normalization condition Tr[t a t b ] = 1 2 δ ab . The derivative with respect to the gauge link U in the direction of t a is defined as ∂ a f (U ) := ∂ ∂ω f (e iωta U )| ω=0 . Negative signs in the subscript of the plaquette variables refer to reversed directions in their orientation, e.g.
A necessary condition for a critical configuration is a vanishing gradient of the action. In the following we derive relations that constrain possible plaquette values from a critical configuration. Eq. (3) vanishes ∀ a, if the matrix in round brackets is proportional to 1 for plaquette values in SL(N, C). For plaquette values in GL(N, C), the matrix has to be zero. For a proof, see the App. D. This criticality condition yields relations for adjacent plaquettes sharing one link. Note that in d dimension one link is shared by 2(d − 1) plaquettes.
We exemplify the d = 2 case: For plaquette values in GL(N, C), we can directly read off the relations and P 1,2 (x) = P −1 −1,2 (x) or P 1,2 (x) = −P −1,2 (x) . (4) If we assume that our critical configuration consists of commuting (Abelian) link variables, i.e. if link variables are diagonal, the above relation simplifies to with ν ∈ {1, 2}. It follows that each critical configuration exhibits at most two distinct plaquettes values. If plaquettes take values in SL(N, C) we get, e.g., for a link in1-direction for an arbitrary α ∈ C. For the Abelian case, this equation reformulates again to a relation between adjacent plaquettes. Restricting this further to the original group SU(N ), i.e. the maximal torus of SU(N ), we find For d > 2, we obtain the same constraint on the imaginary parts of diagonal entries of adjacent plaquettes, but naturally there are more adjacent plaquettes. From the periodic boundary condition we can derive further constraints. Under the assumption, that the relevant critical configurations are Abelian, the product of all plaquettes in an arbitrary 2-dimensional hyperplane Λ µν must be one, i.e. x∈Λµν P µν (x) = 1 . (8)

C. Locality and critical manifolds
We are left with a fairly large number of critical configurations, which we will boil down to a set of basis configurations, getting the others by transpositions and symmetry relations. The ultimate aim is to obtain a homotopic covering of the original integration domain [U(N )] dV or [SU(N )] dV .
As a guiding inspiration, we start with a discussion of the one-plaquette model, i.e. we just have only one plaquette degree of freedom, as Eq. (1) is a sum of local plaquette terms. The action is defined as where an irrelevant constant has been omitted. The Lie derivatives with respect to P are given by Observations are: (i) Eq. (10) vanishes for self-inverse plaquettes P in U(N ). Therefore, the critical points consists only of matrices whose eigenvalues are +1 and −1.
(ii) For P ∈ SU(N ) Eq. (7) implies that the imaginary part of all eigenvalues must be identical. For vanishing imaginary part we obtain the self-inverse elements of SU(N ). The constraint det[P ] = 1 implies an even number of (−1)-eigenvalues, denoted by N (−) . For non-vanishing imaginary part, the center elements of SU(N ) are solutions. For N ≥ 6, we have roots of unity apart from ±1 with the same imaginary parts. So a mixing of these yielding a unit determinant is a valid solution.
Different critical points indicate different action values. Their importance with respect to the weight factor w ≡ e −Re (S) may be exponentially suppressed. For the oneplaquette model we find the following hierarchy of critical points: For P ∈ U(N ), we have The importance of a critical point thus shrinks with the number of (−1)-eigenvalues. For SU(3), we find 6 critical points with three different action values Critical points from the one-plaquette model can be used to construct certain critical configurations for the full lattice theory. We pick critical plaquette values from the one-plaquette model and distribute them in accordance with Eq. (4) and (6) over the lattice. Each critical configuration obtained in this way is one representative of a critical manifold, which consists of its gauge orbit and additional zero modes of the action. For this simple procedure we obtain an additional constraint from Eq. (8): The product of eigenvalues at every position over every two-dimensional hyperplane must be one. We are therefore limited to configurations, where the number of (−1)-eigenvalues at a given diagonal entry is even (the U(N ) case) in every hyperplane or in principal their overall product is one including additional roots of unity.
Note however, that not all critical configurations can be found by the above prescription. Recall that Eq. (4) restricts the plaquette values in a hyperplane at any given position to only two possible values. We can construct a further critical configuration by setting one (or more) plaquette values in that hyperplane to e i , while choosing e i(π− ) for all remaining plaquettes in the hyperplane. Possible values are constrained by Eq. (8), and hence where k is the number of plaquettes with the specific diagonal entry e i(π− ) . The 2πl factor stems from the 2π periodicity. For d = 2, N = 1, l is the actual topological charge. It is constant on the thimble, since the anti-holomorphic gradient flow can be seen as an analytic continuation of the classical gradient flow, which leaves the topological charge invariant [44]. Note especially for k = V /2, that has to be a real number, so the critical configurations are all in the original group space. Picard-Lefschetz theory tells us now, that if the tangent space of the thimble is not normal at this point to the original group manifold, then the intersection number is non-zero (in our case one). We will see, that for Im β = 0, this is the case and they all contribute. For SU(N ), we have to add the constrain that the determinant equals one, effectively reducing the number of critical configurations.

D. The Takagi decomposition and generalized Lefschetz thimbles
Next, we construct the tangent spaces at each critical manifold described Sec. II C. To that end we solve the Takagi equation Here H denotes the Hessian of the action evaluated on the critical manifold. Modes ξ associated with positive λ, point in the direction of the thimble. In the following we refer to such a mode as Takagi vector. Correspondingly, a mode ξ associated with a negative λ points in direction of the anti-thimble. It is called anti-Takagi vector (see e.g. [38]) The λ = 0 vectors do not change the action and are therefore referred to as zero-modes. They result for instance from the gauge degrees of freedom. The Hessian of the action can be written as where the plaquettes appear in different orientations depending on the position of the referred link. Since by construction we have considered only those representatives of our critical manifolds which have diagonal links.
They are Abelian and we can permute our variables to bring all generators to the right.
For critical configurations, where the plaquettes are elements of the original group, the Hessian splits into a real matrix with a complex prefactor H = βM , whose eigenvectors v and eigenvalues α can be computed. For α = 0, the solutions to Eq. (14) take the form where the ξ (+) (ξ (−) ) indicate the thimble (anti-thimble) directions. The α = 0 eigenvectors correspond to the λ = 0 solutions of Eq. (14) and can have an arbitrary complex prefactor.
We observe that for most of our critical configurations. We will come to the exceptions in Sec. III B. The Hessian does not change under field transformations in the direction of these zero modes. Consequently the critical manifold {U crit,0 µ (x)} is independent of those. Therefore, we can deduce that the projection of the subspace spanned by its zero modes in the Lie algebra is the critical manifold itself.
The k index enumerates the different zero-eigenvectors v of matrix M . Forb k ∈ R, this is a compact manifold. In complexified space withb k ∈ C, we have non-compact imaginary directions. Nevertheless, the manifold is still critical. To obtain the generalized Lefschetz thimble [43], we start with a compact submanifold (a cycle) of real dimension #(α = 0) = n (0) . Therefrom its tangent space is spanned. The compact submanifold is commonly called gauge orbit. At every point of this cycle, we use the α = 0 Takagi vectors to span the rest of the tangent space. Since the latter is invariant under zero-modes a point on this tangent space can be directly written in terms of Here the summation index k runs over all eigenvectors.
The real dimension n of the thimble thus splits into zero and non-zero mode directions as n = n (0) +n (+) = dV N g , where N g is the dimension of the Lie algebra. This construction can be generalized to more complicated Hessians, see App. A. If we do not conform to that construction, e.g., tilt the real vectors by a complex factor, we loose homotopy to the generalized Lefschetz thimble. The resulting manifold is non-compact, see Fig. 2. However, the noncompact directions refer to zero-modes, i.e. they leave the action invariant. This corresponds to multiplying the partition sum with an infinite volume factor, which drops out for all physical observables, when calculating expectation values. The prerequisite for these observables is, that they are independent from the zero-modes, i.e. that they are gauge-invariant observables.
The global minimum among our critical manifolds is Since this is a minimum 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 also for the real zero-modes. Since the eigenvectors of M span the whole R n space and since we chose the same complex prefactor for all directions, we can make a basis transformation to the unit basis {e 1 , . . . , e n }. Moving in one of these directions e i corresponds to a change of a single (local) color degree of freedom. This basis can thus be used to construct a local update algorithm on the generalized main Lefschetz thimble, see Sec. II H. As we are interested in a single homotopic covering of our original compact integration domain, we need to limit the tangent spaces. It turns out that the construction of a continous manifold that is homotopic to the original integration domain is the main conceptual difficulty in this approach. For the U(1) one-plaquette model, we have only two critical points and tangent spaces. They give, when glued together at their intersections, a manifold homotopic to U(1), see Fig. 3. We introduce boundaries on the main tangential manifold by identifying intersections with the tangential spaces of the other subleading critical manifolds. We will call these limited tangent spaces, tangential manifolds (TM) from now on. For the full lattice theory we discuss several possible choices of boundaries throughout this work.

E. Hierarchy of critical manifolds
The choice of critical manifolds introduces a natural hierarchy which is reflected by the values of the action. Their importance decreases according to their weight factor e −S with increasing action. Given our choice of critical configurations with diagonal plaquettes from the original integration domain (U(N ) or SU(N )), we can easily express the action in terms of the plaquette eigenvalues. We find where the φ (k) x is the angle of the k-th eigenvalue of the plaquette P µν (x). These critical action values are minima of the attached thimbles, since the real part of the action naturally increases if one moves away from the critical manifolds for Re (β) > 0 in Takagi direction. This is still true, if one considers the Thimble tangent space in a region around the critical manifold, limiting it to a TM. The main critical point (φ (k) x = 0 ∀x, k) defines the global minimum of the action in this many TMs scenario. Consequently with increasing Re (β), certain TMs become exponentially suppressed and we obtain a pronounced hierarchy with the main TM as leading order. For purely imaginary values of β, every thimble contributes equally.

F. Update algorithm on a TM with a Takagi basis
Next, we discuss possible sampling algorithms that are restricted to a single TM. Following our strategy to span the tangent space by the Takagi vectors and real zeromodes, we readily have a parametrization at hand. According to Eq. (18) we can express each configuration on the tangent space by real coordinates b k , specifying a vector in the Lie-algebra. Note that for the zero-modes (α = 0) we have no complex prefactor. Using these coordinates, one can think of applying various different update procedures on the tangent space, starting from a single random walk Monte Carlo (crude Monte Carlo), over a Hybrid Monte Carlo to a trained flow-based neural network [45,46]. The important steps are: 1. Choose a critical configuration, to specify the tangent space on which to carry out the updates. This configuration may also serve as starting configuration.
2. Calculate the real Hessian M and determine its eigenpairs (α i , v i ). This has to be done only once for a given critical configuration, independent of the coupling β.
3. Propose a new configuration by drawing a set of real coordinates b k . The proposed configuration is than specified according to Eq. (18). Perform an accept/reject step based on the real part of the action difference between the old and new configuration. If the proposed configuration is outside the boundaries of the TM, we assign to it a zero probability and reject the proposed configuration, recording the old configuration like in the Metropolis algorithm.
4. Finally take the remaining sign problem into account by reweighting with the imaginary part of the action.

G. Leading and subleading thimbles
So far, we have restricted the sampling to a single TM attached to a specific critical configuration. Ultimately, we are interested in sampling a compact manifold that is homotopic to the original integration domain. In the standard thimble decomposition the combination of various Lefschetz thimble leads, by definition, to a multi-modal probability distribution. Sampling such distributions by a Monte Carlo procedure is difficult. As a solution, a tempered sampling procedure was proposed [33,34]. Another possibility are independent Monte Carlo processes on each thimble. However, in this case the relative weights between thimbles need to be known. One possibility to infer these values is by using prior knowledge of a physical observable for normalization [35].
Here we construct a homotopic manifold by piecewise definition, where we use the TMs as building blocks. As our construction deviates from the thimble decomposition especially close to the boundaries, we do not have infinite action barriers between the patches. Hence, a sampling procedure that proposes configurations across boundaries would be in principle possible. However, we found that it is most convenient to sample the main tangent space with a single Monte Carlo chain and take all remaining patches into account via reweighting. We exemplify this procedure for a system of two TMs, τ 0 , τ 1 .
Calculating expectation values over this extended region requires the relative weight Z 1 /Z 0 . Here Z 1 denotes the partition function corresponding to the subleading tangential manifold and Z 0 refers to the corresponding quantity on the main tangential manifold. It holds Following the method proposed in [36], we introduce a mapping It maps configurations from one of the two patches to the other. With this mapping we can express the ratio Z 1 /Z 0 as It remains to find a suitable f . Since we consider only tangent spaces, f is linear and can be constructed as a basis transformation from the Takagi basis of τ 0 to τ 1 . The Jacobian det[df ] is therefore a constant factor. As the eigenbases of the real Hessian M are orthogonal and can be chosen to have determinant one, the Jacobian depends consequently only on the the different sets of complex prefactors c k of the Takagi vectors. These depend in turn only on the sign of the eigenvalues α k . Therefore we find for the Jacobian where n (+) and n (−) are the number of positive/negative eigenvalues of the real Hessian M (see Sec. II D) at the respective patch. Here we assumed that τ 0 is the patch attached to the main critical manifold having only positive eigenvalues being the global minimum. If the patches τ 0 and τ 1 are of different size, one may introduce an independent scale parameter to the mapping f , which is than also reflected as factor in the Jacobian. As the main tangential manifold is usually the largest such a factor is not necessarily needed in practice. It may however be used for optimization purposes.
H. Alternative updates on the main tangent space As outlined in Sec. II D we can sample the main tangent space by just setting each complex factor c k to β * /|β|. In practice this tilts every link, since in this case all eigenvalues are positive or zero. Having identical complex prefactors, we can perform a basis transformation of our coordinates to the unit basis. Therein each coordinate corresponds to a single gauge degree of freedom. This allows for applying a local heat bath or Metropolis algorithm. The situation discussed here is shown at the bottom of Fig. 2. The unbounded critical manifold in imaginary direction has to be dealt with.
We limit the plaquette values in accordance with the intersection points of the tangent spaces in the one plaquette-model, creating TMs (see Fig. 3). Similarly as before, we define the region outside the boundary to have zero probability.
But since link variables can still diverge in imaginary direction we have to apply a second limit or just record variables, which are unaffected by the zero modes. In our theory, these are the plaquettes variables. An alternative is to adopt gauge cooling [47] to make this problem milder. A severe limitation of this method stems from the fact that only observables invariant under all zero-modes can be measured. Taking pure Yang-Mills theory it is not possible to measure the Polyakov loop. It is invariant under gauge transformations but not under a global zero mode. The latter is represented by changing all links in the same direction and amount leaving the plaquettes and therefore the action invariant.

III. APPLICATION TO A 2-DIMENSIONAL U(1)-GAUGE THEORY
In this section we apply the above outlined method to two-dimensional pure U(1) gauge theory with periodic boundary conditions. Recently, complementary studies on this theory with a sign problem have appeared in the literature. In [48] the complex Langevin method is employed. The complex action here is caused by a non-zero vacuum angle. On the other hand in [49] the theory with a complex gauge coupling is investigated by means of the path-optimization method.

A. The effective degrees of freedom
We formulate the theory in terms of its effective degrees of freedom. Eq. (8) can be verified easily by noting that every link appears twice in all plaquettes. Hence the links cancel each other when being multiplied altogether. Consequently, one plaquette can be expressed in terms of all others. The action of the two-dimendional theory is rewritten as follows neglecting constant terms. The last term is called toron term. One can reformulate the full theory in terms of these plaquettes variables having a reduced partition sum, which gives the same expectation values for observables that are invariant under zero modes

the integral factorizes in plaquettes variables yielding
The plaquette expectation value is then which is exactly the the same as for the one-plaquette model. There is no volume dependence. Since the difference between periodic and open boundary conditions vanishes in the infinite volume limit, this is the expected value.
We construct the approximation scheme by successively including TMs as integration domains. Thereto we split the integral according to the critical points P = ±1 of the one plaquette model and get Z = τ0 dP e β/2(P +P −1 ) + τ1 dP e β/2(P +P −1 ) We can map this to the lattice by k denoting the number of plaquettes being −1 at the critical configuration. V k is the number of such combinations on the lattice.
Eq. (28) allows to calculate approximate values for the comparison with numerical simulations (see Fig. 4).
Complementary, there exists a formal solution for the for the full lattice theory with periodic boundary conditions involving no approximations. We can write the partition sum as being a series in modified Bessel functions I n (β), where V is the number of plaquettes [50,51]. The leading order of this series corresponds to our approximation. The following orders take finite volume effects into account and yield the exact result provided that the series converges for the given value of β.

B. Critical manifolds and their hierarchy
For an even number k in Eq. (28) the critical configurations in the approximation and in the original lattice theory coincide. In contrast, for odd k, this does not hold due to the periodic boundary conditions. However, it is possible to construct critical manifolds being (arbitrarily) close to the corresponding configurations in the approximation. For odd k < V /2 − 1 we set k plaquettes to e i(π− ) and the rest to e i such that Eq. (13) is satisfied. This leads to the choice which we refer to together with the critical configurations for even k as basis configurations from now on. Furthermore as discussed in Sec. II C, Eq. (13) gives rise to a multiplet of configurations being related to the basis configurations by a symmetry. To see this, note that these configuration have topological charge k 2 . Now, we can change the topological sector by adding/ subtracting 2π V −2k from . We observe that if | | < π/2, the Takagi vectors do not change. Consequently we can apply this transformation directly to measured plaquettes values by multiplying them with e ±2π/(V −2k) , where the sign is cho-sen whether we have a e i(π− ) or e i plaquettes. This transformation does not only apply to odd k but also to even k. By this procedure, we reach 2(V /4 − k 2 ) critical manifolds in different topological sectors for odd k and 2(V /4 − 1 − k 2 ) + 1 for even k, respectively. Second, we can go from the k-th TM to the (V−k)-th TM by changing the plaquettes from the critical configuration according to The corresponding real Hessian M → (−M ) changes sign and has the same eigenvectors and zero modes. Only the non-zero eigenvalues α → (−α) change sign. Consequently and using the fact that eigenvectors are unique up to a non-zero scalar multiplication we get the Takagi vectors for the (V − k)-th TM by multiplying the Takagi vectors from the k-th TM by √ −1 = i. For k = V /2, V /2 ± 1, all zero modes do not change the plaquettes. Writing a plaquette configuration on the k-th tangent space as P crit 1,2 (x)e i∆ϕ(x) , we can write the mapping to the configuration on the opposite (V− k)-th tangent space as having again a transformation for directly measuring the plaquette values. Having these, the same procedure to reach the different topological sectors can be applied to these plaquette values.
An exception is the case k = V /2, which has an additional zero mode, which can be parametrized by P 12 (x) = e iϕ and P 12 (x) = e i(π−ϕ) , for V /2 plaquettes on either side. This transformation necessarily leaves the action invariant while changing actual plaquettes values. The configurations, where all plaquettes are ±i, which would be assigned to k = V /2 ± 1 in our scheme are included in this critical manifold and are therefore being left out. This transformation reduces the combinatorial factor to 1 The critical manifolds form a hierarchy depending on their associated value of the action for the basis configurations. Depending on Re (β) the critical manifolds differ in importance for the partition sum since on thimbles and suitably bounded TMs the real part of the action is minimal at the critical manifold. Consequently, if β is purely imaginary, every thimble or TM contributes equally. Otherwise one can obtain an approximate result by taking only a few thimbles or TMs into account as the others are exponentially suppressed. considering subleading TMs up to order k = 4 restricted to the basis configurations. The simulation setup used for the shown data is described in detail in Sec. III D and III F.

C. Ensuring homotopy
To ensure global homotopy, we need to make sure, that the TMs form a patchwork covering [U(1)] 2V . Therefore, we have to create boundaries, which match each other exactly. Strictly speaking, this is not generally possible in higher dimensions, since there is no theorem that tells us, that these tangent spaces have to intersect in this manner as it is the case for thimbles. But we can at least get close to something alike minimizing the systematic error introduced by homotopy violations as much as possible. Our approach comes with thinking in effective degrees of freedom being plaquettes variables with a toron term as shown in Sec. III A. Looking at the different (subleading) tangent spaces, we have applied several schemes, each based on different criteria: 1. Real plaquette boundaries based on the one plaquette model: These are shown in Fig. 3. The plaquette variables take values in a region bounded by the intersection points of the two tangents in the oneplaquette model. In the full lattice theory, we still find these intersection points for the transitions in configurations where k is even. We therefore limit the real parts of the plaquettes by these intersections. For = 0 configurations the transition looks different and we take it into account by a shift of the boundaries. Having transition points does not mean that we get full homotopy. We would need to identify intersections with real dimension 2V − 1 which do not necessarily exist.
2. Imaginary plaquettes bounds: The limit can also be applied to the imaginary part of the plaquettes preventing them from drifting to far off in imaginary direction. This allows a larger space to be explored than for real plaquettes bounds. Moreover these boundaries are far easier to implement since all real critical manifolds are part of the original group and we only need one value to specify these boundaries. However, this doesn't care so much for homotopy like the real boundaries, since the plaquettes in the full lattice theory only approximately lie on the same tangent spaces found in the one plaquette model. Consequently overlapping TMs are not excluded in the case of imaginary plaquettes boundaries, while the real plaquette boundaries still guarantee that there is no Overlapping.
3. Action boundaries: On a Lefschetz thimble, there exists a coordinate system with the critical point at the center, where S R is a simple rising quadratic function [52,Lemma 2.2]. Since the TM is close to it in the vicinity of the critical manifold, we observe the similar rising behaviour up to some distance illustrated for the local action in Fig. 6. We limit the plaquettes variables by the local maxima of S R . Since configurations are exponentially suppressed by it and the fact that this distance goes beyond the intersection points mentioned before allows a larger part of configuration space to be explored, while the systematic error stays small.

Spherical boundaries:
This is the most conservative choice being independent from the coordinate system. We observe that wandering along the Takagi vector with the highest eigenvalue α = 8 (this is true for all even lattice sizes including and greater than 2 × 2) for the main TM turns it into the ultimate subleading TM (whose critical configuration has all plaquettes equal to −1) intersecting with its counterpart (α = −8) from there. This can be understood in terms of the effective d.o.f. as a diagonal connection in a hypercube. The intersection marks a corner in the cube containing the main TM. We now choose the radius of the inner sphere of the cube to make sure, we do not intersect with other TMs, where we can assign a radius in the same way. We get as the radius of the effective sphere (the distance from the critical manifold is calculated in terms of the non-zero eigenmodes). A disadvantage is that here the curse of dimensionality hits directly into the calculations, since with rising dimension, the volume of the inner sphere becomes negligible in comparison to the cube and we explore only small portions of configuration space. One can counter that by scaling the radii paying the price of overlapping TMs introducing another systematic error.

Im(S) boundaries:
The idea here is that the TMs go into the other TMs by their intersections. So the imaginary part of the action of one TM has to change continuously to the one of the other (On thimbles Im S is constant and changes abruptly at their singular intersection, which in our case is at infinity.). Having a pronounced hierarchy allows us to set the boundaries using this feature by e.g. allowing one tangent space only to vary in an interval of Im S and the other one on a subsequent interval. Problematic is the fact that we have intersections of one TM with multiple TMs in different orders. So here we possibly limit the explored space too much.
All in all, we have found a combination of real plaquettes boudaries with shifts and action boundaries most promising. To that end we calculate the former and correct them, if they extend futher than the action boundaries, which is especially important for TMs with high . Otherwise we can fall into the regions with negative action, 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 to generate configurations on the main tangential surface.
(1) Diagonalize the Hessian at the main critical point where all plaquettes are +1.
(2) Construct the parametrization of the main TM surface. The eigenvectors from (1) with non-zero eigenvalues correspond to thimble directions. Tilt those as described in Eq. (16) above to obtain the Takagi basis {ξ i i = 1, . . . , 2V }. (ii) A proposed configuration outside of the specified boundaries is being rejected.
The measured configurations are stored to be used for the reweighting on the subleading TMs. This part of the algorithm is described in detail in Sec. III F.
Like already mentioned in Sec. II H, there is also a local update algorithm (see App. C).

E. Numerical results on the main TM
Looking at the approximation in Fig. 4, we expect the main tangent results to roughly follow the zero order approximation. For high β R , the other tangent spaces are exponentially suppressed allowing for convergence to the full result. Noticing that the full result for this range is slightly above the one plaquette model due to finite volume effects, we hope to see the same from the simulation, which is the case, see Fig. 7. This proves that we do not accidentially just simulate the one-plaquette model by our procedure.
Another important point, is that the sign problem should be reduced in comparison with standard reweighting, which is also the case (see second plot in Fig. 7).
So, we can expect for high enough β R (e.g. here larger than 3) correct results for simulations at complex beta with a lesser sign problem. To extend this range, we need to take the subleading TMs into account, see Sec. III F.
For standard phase reweighting, the average sign should exponentially decrease with increasing space-time volume. Since our simulation works similarly just on a tilted space, we see the same happening here. The difference is that our average sign is higher than the one for standard reweighting and the slope is less steep, see Fig. 8. Therefore, higher lattice volumes are more easily accessible in our approach. Indeed, we needed to increase statistics for the reweighting simulations, while the number of samples for the different volumes in the TM simulations remained the same.

F. Reweighting of subleading TMs
For the reweighting, the observable from Eq. (20) becomes where n denotes the number of subleading tangent space orders, one wants to incorporate, m k are the multiplicities and f k are linear transformations projecting the main TM τ 0 onto the subleading tangent space τ k . This is done by aligning the Takagi basis using the real vectors v j to calculate a transformation matrix which we apply to the subleading Takagi basis where the c The z (k) i are now our new aligned basis vectors for τ k . Therefore we can directly project from the leading to subleading tangent spaces. As already mentioned in Sec. II G, depending on the other boundaries of the subleading tangent space, one can apply scaling factors to the variables. We observed, that for β R , β I ≥ 0, the main TM is always larger (or equally large) as the other subleading TMs. Instead of rescaling, we use an indicator function χ τ k (U ) for the boundaries: If the projected  Takagi ord0  Takagi ord5  Takagi ord10  Takagi ord15  Takagi ord20  Takagi ord25  Takagi ord30  Takagi ord31  Takagi Takagi ord0  Takagi ord5  Takagi ord10  Takagi ord15  Takagi ord20  Takagi ord25  Takagi ord30  Takagi ord31  Takagi with g(U ) an arbitrary function on the subspace. So, we simply set the integrand to zero in the reweighting process, if the projected configuration is out of bounds. We use several symmetries to incorporate the different topological sectors for each order as well as a mapping from the k-th to the (V − k)-th TM already discussed in Sec. III B. In the end, we have to diagonalize only V /2 Hessians to get every contribution. For practical reasons, we subsumed the V /2 TM under order V /2 − 1 and V /2 using a turn by ϕ = π of the additional zero mode (33). During the reweighting process, we need to check for every contribution, if it is in its predefined boundaries, since they also differ over topological sectors (see Sec. III C).
We applied the procedure on an 8×8 lattice at constant Im (β) = 1 and compared different orders of reweighting with the BesselNLO result (see end of Sec. III A). The results are shown in Fig. 10 as well as the average sign in Fig. 9 depending on how many orders of TMs one takes into account. The boundaries chosen are real plaquette boundaries based on the one plaquette model, since they prevent overlapping of the TMs, while allowing a large space to be explored.

IV. DISCUSSION
As we have stated already in the introduction, the particular choice of our deformation has two important advantages: (A1) Due to the flatness of the patches a parametrization in terms of real coordinates and basis vectors is easily constructed. It is thus straight forward to realize a sampling procedure on a particular patch.
In contrast to the generalize Lefschetz thimble approach, there is no need to solve a flow equation to propose a new configuration nor to evaluate a Jacobian at each sampling point.
(A2) Since our coordinate systems originate from a critical configuration on each patch, which is the configuration that receives the largest weight on that patch, a reweighting from one patch to another is possible without any severe overlap problem.
In the case of a pronounced hierarchy among critical configurations, we have demonstrated in the case of the 2d U(1) theory, that a well controlled approximation scheme emerges if successive contributions from suppressed patches are taken into account. On the other hand there are two -as we believe less severe -disadvantages: (D1) Sampling on the flat tangential manifold rather than on the curved thimble, does not erase the sign problem completely. It would be interesting to analyze whether the optimal deformation of the original manifold, which minimizes the combined sign problem of the action on the integration domain and the one introduced by the Jacobian, is closer to the thimble decomposition or our decomposition from flat patches. We leave this for future investigations. We have demonstrated that the resulting sign problem in the case of the 2d U(1) gauge theory is, although sill exponential in volume, very mild compared to the standard reweighting procedure.
(D2) In order to construct our deformation, we need to know all relevant critical configurations. Depending on the space-time dimension, volume and gauge group, this can be a very large amount of critical configurations.
It is one of the main results of this work that we have identified a large number of critical configurations, which are characterized by specific distributions of plaquettes across on the lattice. In particular, we have demonstrated that we can chose these critical configuration from the maximal torus of the original gauge group. We have shown further that there exists a large number of degenerate critical configurations and patches due to lattice symmetries. For the reweighting procedure we thus need to take only one of these patches into account if appropriate combinatorial multiplicity factors are used. The results look promising and systematic errors by nonhomotopy vanish with growing lattice size. We want to emphasize here that the choice of a complex coupling β is not at all a pathological choice. In the limit of a purely imaginary coupling, the kernel of the discussed partition sum can be seen as the real time propagator of the theory. For the evaluation of real time correlation functions in a thermal bath, the Schwinger-Keldysh formalism is usually applied. Our sampling strategy might also be applied to the Schwinger-Keldysh contour, even though in this case an additional sign problem arises from the edges of the contour.
Moreover, the critical configurations we have identified here are not only critical in the case of the Yang Mills action with complex coupling β. The same configurations remain critical when we introduce fermionic matter fields with a chemical potential µ. In this case the effective action might be written as S eff (µ; U ) = βS G (U ) − Tr ln D(µ; U ) .
Hence, the action gradient is the sum of a contribution from the gauge (S G (U )) and the fermionic part (S F (µ; U )) of the action. That the action gradient of the gauge part vanishes at our critical configurations has been discussed in detail above. The fermionic contribution to the gradient is given as ∂ x,ν,a S F = i Tr D −1 e µδν,0 U ν (x) − e −µδν,0 U † ν (x +ν) t a . (40) For those critical points that are not only chosen from the maximal torus of the gauge group but are also constructed from center elements of the original gauge group, the fermionic contribution vanishes as well. As all link variables are proportional to the unit matrix, the (none sparse) inverse of the fermion matrix D −1 contains N ×N diagonal blocks which are also proportional to the unit matrix. We conclude that the matrix multiplying the generator t a is proportional to the unit matrix and as such the whole expression vanishes. We thus hope that our strategy might also prove useful in this case, i.e. the field theoretical description of dense matter, including QCD at net-baryon number density.

V. CONCLUSION
We have put forward here a novel nonperturbative lattice approach for Yang Mills theories with a complex gauge coupling β. The approach is based on a deformation of the original integration domain of the theory into complex space. Guided by the Lefschetz thimble decomposition of the partition sum we have chosen a new integration domain. This is constructed piecewise from patches of tangential manifolds to the relevant Lefschetz thimbles.
For the numerical implementation we have chosen to set up a Monte Carlo procedure on the main tangent space only. All further contributions from other patches are taken into account by reweighting. We have tested our approach by applying it to the case of 2d U(1) gauge theory. Here, it far out-performs our benchmark simulation with standard reweighting, see Fig. 1.
Based on this observation we plan to apply our approach to general U(N ) and SU(N ) gauge groups in 4d. While our approach has shown to feature an exponentially better performance than standard reweighting, it remains to be shown that the sign problem stays numerically manageable also in these cases. We hope to address these questions in a forthcoming publication. Moreover, we envisage simulations of fermionic matter fields at finite chemical potential as well as real-time lattice theories with expectation values along the Schwinger-Keldysh contour. Since i and j were arbitrary M = c1 for some c ∈ C.
For T ∈ gl(N, C) sup sl(N, C), T is not traceless anymore and can also be e.g. 1, which excludes the M = c1 case and M has to be zero to fulfill the statement.