Reweighting Lefschetz Thimbles

One of the main challenges in simulations on Lefschetz thimbles is the computation of the relative weights of contributing thimbles. In this paper we propose a solution to that problem by means of computing those weights using a reweighting procedure. Besides we present recipes for finding parametrizations of thimbles and anti-thimbles for a given theory. Moreover, we study some approaches to combine the Lefschetz thimble method with the Complex Langevin evolution. Our numerical investigations are carried out by using toy models among which we consider a one-site z^4 model as well as a U(1) one-link model.


I. INTRODUCTION
QCD at vanishing and finite temperature is one of the best tested theories in high energy physics. At the present moment theoretical predictions from first principle lattice simulations match remarkably well with experimental data for instance from heavy ion collisions, see e.g. [1][2][3]. However, at finite chemical potential lattice simulations suffer from the sign problem, which a priori prohibits simulations based on importance sampling. This severely limits the access to the largest part of the QCD phase diagram. By now, there are many approaches towards a solution of the sign problem. Overviews addressing developments in finite density QCD over the last years can e.g. be found in [4][5][6][7][8].
Amongst those are Taylor expansions [9], simulations at imaginary chemical potential [10,11], reweighting [12], the density of states method [13], dual formulations [7], the Complex Langevin method [14,15] and the Lefschetz thimble method [16,17]. So far none of those methods have been able to give reliable results for µ/T 1. In this work we focus on the Lefschetz thimble approach. In Euclidean space-time the Lefschetz thimble approach has been applied to bosonic theories as well as to (lowdimensional) QCD in [17][18][19][20][21][22]. Recent applications to fermionic theories (such as the Thirring model) can be found in [23,24]. Moreover, field theories in Minkowski space-time formulated on the Schwinger-Keldysh contour have been studied using the thimble formalism in [25]. Algorithmic improvements to the holomorphic gradient flow method were proposed in [26]. Recent contributions in the field more generally involve complex manifolds close to Lefschetz thimbles that are optimized such that they ameliorate the sign problem, see e.g. [27][28][29][30].
As we see later, the CLE and the Lefschfetz thimbles are closely related. Studies investigating the interplay and the connection between the two approaches can be found in [31,32].
The Lefschetz thimble method relies on a deformation of the integration path. By construction the imaginary part of the action is constant on the transformed paths, which are called Lefschetz thimbles. There are two ba-sic algorithmic frameworks [17,23] providing recipes for Monte Carlo simulations on the Lefschetz thimbles. The first employs Monte Carlo simulations directly on the thimbles. The latter continuously deforms the original integration path close to the actual thimbles to lessen the sign problem.
In this work we address a few key challenges to the the Lefschetz thimble method and propose algorithmic improvements. One of the main problems with Monte Carlo simulations on Lefschetz thimbles is to determine the weights of the thimbles relative to each other. This difficulty arises as the original path integral is decomposed into a sum of integrals over multiple thimbles. We show that this difficulty can be overcome by a standard Monte Carlo determination of the ratios of the real partition functions on the thimbles. This is facilitated by a novel reweighting procedure which is generally applicable to field theories. In this work we assume prior knowledge on a parametrization of the contributing thimbles. To find this parametrization we propose two algorithms which can be generalized to higher dimensional theories. However, the reweighting procedure does not rely on knowing a parametrization. Our ideas are put to work in simple models, i.e. ordinary integrals. Among those we consider a one-site quartic model with a λ 4 z 4 term as well as a U (1) one-link model.
The paper is organized as follows. We start by briefly revisiting the idea behind Lefschetz thimbles, see Sec. II. In Sec. III we propose two algorithms to find thimbles and their parametrizations necessary for Monte Carlo integration. In Sec. IV we present our idea of sampling on multiple thimbles taking into account the relative weights of different thimbles. Sec. V introduces the toy models we use for numerical investigations together with our results. We conclude this paper in Sec. VI. During the research for this paper we have also developed many ideas to combine the Complex Langevin evolution and the Lefschetz Thimble method. While none of those approaches lead to generally applicable algorithms, they still provide some useful insight into the structure of the models, hence we give some of those ideas and corresponding results in App. C.

II. THE LEFSCHETZ THIMBLE METHOD
In the following we briefly revisit the Lefschetz thimble method. The idea behind the approach is to rewrite the path integral measure over a real manifold [16] to circumvent the sign problem by allowing Monte Carlo sampling on this manifold [17]. Here we explain the Lefschetz thimble method by using the example of simple one-dimensional integrals. Consider a complex action of a real variable S(x). Next, extend the real axis to the complex plane R → C, i.e. x → z = x + iy. Given the stationary points z σ of S(z) one can define a real path in the complex plane D σ ⊂ C as the solution of the steepest descent equation ending at z σ this path is called a Lefschetz thimble. The action has constant imaginary part along the thimble. The integral can be decomposed into integrals over all D σ where D ⊂ R is the original real domain and n σ is the intersection number of the steepest ascent path (unstable thimble) with the original domain D. One can now formulate Monte Carlo algorithms based on (3), see e.g. [17,23,32]. Observables are then computed in the usual way with the only difference, that it has to be computed on every thimble. There are two practical problems with this approach: 1. Finding all contributing thimbles can be a challenging task.
2. For the case of multiple contributing thimbles there is no simple way so far to access the relative weights, i.e. the ratio Z σi /Z σj .
Both are challenging problems without a general solution so far. An approach to the first problem is the holomorphic gradient flow, which approximates the thimble structure, by simulating on a deformation of the original domain close to thimbles [23]. In [21] the second problem regarding relative weights has been approached by using known results in some parameter regions. In the following we propose general solutions to both of those problems which do not rely on approximations. In this work we demonstrate our solutions by means of simple models.

III. FINDING THIMBLES
In this section we propose two algorithms which can be used to systematically find contributing thimbles. This is put to work in simple one-dimensional integrals. Generalizations to higher dimensions might be expensive. The first algorithm scans the real axis in search of intersecting anti-thimbles, while the second algorithm projects points in the complex plane onto thimbles, in order to determine a numerical parametrization of the thimbles. Both algorithms also apply in higher dimensions, however the numerical costs may rise exponentially with the number of lattice points. This is currently investigated for gauge theory in [33].

A. Axis scan
Since the only contributing thimbles are those with non-zero intersection number of the anti-thimble with the original manifold, one can find all contributing fixed points by scanning the manifold for such intersections. This can be a challenging problem in higher dimensional theories, however importance sampling by Monte Carlo methods in parameter regions without a sign problem or in the phase quenched theory might give good starting points for such searches. In the following we describe the searching algorithm for the case of simple integrals, i.e. the original manifold is an interval [a, b] ∈ R. The algorithm is the following [34], 1. Choose a starting point on the real axis.

Solve the steepest ascent equation
using the starting point as an initial condition.
3. If the derivative of the action becomes small the flow is close to a fixed point of the action.  (5) is solved using a starting point on the real axis close to the anti-thimble, once it is close to the fixed point (i.e. the derivative of the action is smaller than some value δ, visualized by the blue circle), the flow is switched according to equation (8), and will end in the fixed point (solid black arrows). The steepest ascent without switching close to the fixed point will asymptotically approach the thimble (black dashed arrow).

4.
Depending on the structure of the fixed point, one can now reach it by looking at the Langevin flow and changing the sign according to the following prescriptioṅ All those cases have to be tested, and one of them will end in the fixed point.
This algorithm is visualized in Fig. 1.
Once the fixed points are known, the numerical parametrization of the thimbles can be computed. In the case of one dimensional integrals, this boils down to solving one dimensional differential equations. We do so by solving the normalized steepest descent equation with opposite sign starting close to the fixed point [31]. The reason for the normalization with the absolute value will become clear later, however it also helps with numerical stability when solving the steepest descent equation.
Note that this normalization is simply a rescaling of the flow parameter τ .

B. Thimble cooling
In this section we propose a straight-forward algorithm to find parametrizations of all thimbles and anti-thimbles for a given action where the only required information is the knowledge about the fixed points. We show results for the z 4 model with parameters σ = 1, λ = 1, h = 1 + i. For a definition of this model see Sec. V A. The idea is to minimize the distance of any set of points in the complex plane to curves on which the imaginary part of the action is constant. This leads to the following definition.
where σ labels the stationary points. We call M σ (x, y) the cooling function. The numerical minimization procedure is facilitated by the following gradient equationṡ By construction the gradient equations orthogonally project a given point on a curve with constant imaginary part of the action. As initial conditions for (11) we use a grid of points in the complex plane. This is shown in the left plot in Fig. 2. The choice of a random grid is arbitrary. We could have also chosen a regular grid. In the right plot of Fig. 2 we show the result of the solution to the gradient equations. The procedure works well for a large box of initial conditions (red) as shown for the right fixed point (black). The resulting set of points lies on the thimble, the anti-thimble and an additional curve on the left without physical relevance. Alternatively, we can start with a small rectangle around the fixed point, see the blue points in the left plot. From the resulting points flowed to thimble and anti-thimble we can choose the next set of initial conditions along e.g. the thimble and repeat the procedure iteratively. Note that also from looking at the flow lines of (11) we can determine suitable areas for initial conditions. This method provides a useful tool to find parametrizations to the thimbles (and anti-thimbles) by interpolating the flowed set of points. Moreover, by knowing each anti-thimble we can map out if it intersects with the original integration manifold thus enabling us to determine whether and how much the corresponding thimble contributes. In particular, the method could be applied in higher dimensional theories where the minimization procedure is combined with importance sampling around the fixed point. From there thimbles and anti-thimbles can be successively parametrized as illustrated in blue in Fig. 2.
Thimble cooling has the potential advantage over the procedure in Sec. III A that one does not have to solve the holomorphic gradient flow in many directions but one directly flows to the (anti-)thimbles. The cooling method put forward here can also be used to reduce numerical discretizations artifacts (see Sec. V) in the thimble parametrization.
We remark, that a generalization of thimble cooling to higher dimensional integrals may in general prove difficult due to the dimensionality of the hyper-surface parametrized by Im[S(z)] = const. In App. C 3 we propose a combination of Lefschetz thimble and complex Langevin, which samples around all thimbles. This can also be used as a starting point for thimble cooling.

IV. MONTE CARLO SIMULATIONS ON LEFSCHETZ THIMBLES
The previous section dealt with finding parametrizations for the thimbles. In this section we propose an algorithm for simulating on thimbles provided its parametrization is known. We also show how to compute the ratio of partition functions from within Monte Carlo simulations.

A. Reweighting on thimbles
Once the parametrization of the thimble is known, we can simply rewrite the partition function Z σ on the thimble aŝ Dσ dz e −Re[S(z)] =ˆb a dτ e −Re[Sσ(τ )] J σ (z(τ )) , (12) where S σ (τ ) = S(z(τ )) is the action evaluated on the thimble D σ . We have also rewritten the integral to run over the flow parameter in (9) and introduced integral boundaries, which are defined by the domain of τ . This introduces the complex Jacobian J σ = ∂z/∂τ on D σ . For the right hand side of (12) we can apply a real Langevin simulation or Monte Carlo sampling along the thimble. The Jacobian is dealt with via reweighting, Monte Carlo sampling now produces samples according to the distribution So far we have dealt with a single thimble. The reweighting equation (13) for multiple thimbles becomes where we have defined Note that in (15) the thimbles are weighted with their partition functions which have to be determined within the simulation.

B. Computing the partition function weights
Now that we have a simple algorithm for computing observables on the thimbles, we can proceed to the problem of how to compute the weights. With the above definition of S σ and considering only two thimbles for simplicity, we look at the ratio of their partition functions, i.e. we choose one thimble as a "master" thimble and divide the numerator and denominator of (15) by its partition function. The following identity states the ratio of partition functions provided (i) the integrals over the thimbles have the same boundaries and (ii) the flow parameters τ on both thimbles can be identified -if the latter does not hold, an additional Jacobian must be taken into account. For a derivation of (18) see App. A. (i) can be enforced by using suitable variable transformations, see App. B, while (ii) is guaranteed by normalizing the steepest descent equations, as we did in (5). Hence, it is possible to compute the ratio during the Monte Carlo simulation, which is necessary in higher dimensional integrals, e.g. in field theories.

V. APPLICATIONS
We investigate different models with varying complexity to test our algorithms. First we look at a model with only one contributing thimble, which is a good test case for setting up the Monte Carlo simulation. Next we address a model with two contributing thimbles, whose flow parameters run over the same interval τ ∈ [−∞, ∞]. Finally we investigate the U(1) one link model, which is a model for a simple gauge theory with fermions and has thimbles that end in poles. This model is quite general in the sense that it contains all features that are to be expected in more complicated cases such as field theories.

A. One-site z 4 model
The one-site z 4 model generically consists of three thimbles, which end in different asymptotic regions at infinity. This structure makes the model a rather simple test case for the algorithms we propose. The model is given by the action for more details see e.g. [31]. We can choose the models parameters such that there are one or two contributing thimbles, i.e. with n σ = 0: 1. For σ = 1, λ = 1/3 and h = 1 + i there is only one contributing thimble.
Both cases are shown in Fig. 3. The distributions exp (−Re (S)) on both thimbles are shown in Fig. 4, where they have been mapped onto the interval [0, 1], see App. B. Note that this is not necessary for the simulation in this model, since the original domains already overlap. There we see that both distributions fall off exponentially which guarantees numerically stable simulations. Numerical results for the observable z 2 as well as the ratio of partition functions  from (18) are given in Table I, those simulations have been performed without employing the variable transformation for simplicity. For the simulations we have collected (10 10 ) data points. Errors have been estimated via a standard Jackknife analysis.

B. U(1) one link model
After having demonstrated the viability of our methods in a simple model, we look at a more difficult model, namely the U(1) one link model with a finite chemical potential µ. In addition to having multiple contributing thimbles, this model has thimbles that end in poles at finite values of the flow time making it necessary to transform the integrals such that the the distributions on the thimbles overlap. This model is a suitable testbed for our methods, since it contains general features that will also be present in more realistic theories. The model's action reads where κ = 2, β = 1 and µ = 2. Its thimble structure is shown in Fig. 5. This model was studied by means of complex Langevin in [14]. Its thimble structure has been studied in [26,31]. It has three different contributing thimbles, of which two are connected by symmetry, see Note that in this model, only periodic observables which are analytic on U(1) make sense, hence we study the analogue of the Polyakov loop, its inverse, the plaquette and the density, which are given analytically in [14] U = e ix , Simulation results are given in Table II. For the simulations we have collected (10 9 ) measurements for the U(1) model. Again, errors have been estimated via a standard Jackknife analysis.   Our simulation results agree with the exact results from [31]. In Table II we only provide statistical errors. Systematic errors arise from numerical discretization artifacts along the thimble. In the simple cases at hand the latter can be quantified by comparing the exact solution with the result from integrating along the discretization obtained from the gradient flow (11). In the case of the observables given in (21) the deviation is of order 10 −6 . Therefore the systematic error is comparable to the statistical error. The ratio of partition functions seems to be particularly sensitive to this effect. However, by taking into account the statistical and the expected systematic error all quantities agree with the exact result within the error.

VI. CONCLUSIONS AND OUTLOOK
We propose a new method based on Lefschetz thimbles for solving theories with a sign problem. This method works with two steps: First we find all contributing fixed points by scanning the original manifold for intersecting anti-thimbles. We then project a grid of points in the complex plane onto the thimbles -which requires the knowledge of the fixed points -to obtain a numerical parametrization. Second we simulate on these parametrizations and determine the relative weights within the simulation by means of a reweighting procedure. This reweighting can be tuned such that there is no overlap problem. We remark that within this method finding numerical parametrizations of the thimbles may be costly in higher dimensions. The reweighting procedure on the other hand straightforwardly generalizes to field theories and can be combined with other simulation algorithms for thimbles.
The above method has emerged from discussions and investigations of the idea of simulating a complex Langevin evolution, that is directed to the Lefschetz thimbles. This procedure of a Lefschetz-cooled Langevin update has been partially, but not fully, successful. In our opinion such a combined approach still has its potential, more details can be found in App. C.
Our proposal for simulating on Lefschetz thimbles has been illustrated on a one-site z 4 model in Sec. V A. It is successfully tested in a U(1) one-link model at finite density. The results are discussed in detail in Sec. V B.
Interesting future applications are field theories, e.g. the Schwinger model and higher-dimensional gauge theories, see [33]. As there are many complicated steps for simulations on thimbles, it is desirable to find simpler alternatives, which at best can be applied blindly. One natural idea [32,[34][35][36][37] is to combine the complex Langevin evolution with the Lefschetz thimbles. A combination of both equations is only consistent after an appropriate coordinate transformation. The latter can adaptively be generated during the combined Langevin and gradient flows. Due to its similarity to standard cooling algorithms as well as the gauge cooling we call this process Lefschetz cooling.
Despite its full success described below in particular for simple Gaußian models it is only partially successful in more complicated models, notably already the z 4 model. While the method is not fully successful yet, in our opinion it is still a very interesting one to pursue. Its potential power is the self-adaptive local nature of the simulation steps. However, this also poses the biggest conceptual question: How does such a local procedure capture the global nature of the intersection numbers n σ in (3) correctly? Note that besides its formal importance this question could be practically less important as it seems: in most models under investigation so far we have n σ = 1.
Below, we list some ideas putting Lefschetz cooling to work and discuss their viability and applications. While none of those proposed ideas so far have managed to give quantitative correct results for observables, they provide useful insight into possible realizations of the approach.

Variable transformations
We aim to make the complex Langevin evolution compatible with constraints characterizing the Lefschetz thimbles by means of variable transformations. The latter have been investigated in combination with the complex Langevin evolution in a different context in [38]. There it was shown that the complex Langevin evolution including a transformation can give correct results while failing in the original formulation of the problem. Here we pursue the idea of having flow time-dependent variable transformations to transform the complex Langevin evolution towards the thimbles of the theory. This approach is natural in the sense that complex Langevin should already be able to sample the relevant fixed points [31], note that this is in a similar spirit as for the contraction algorithm [23]. By forcing the complex Langevin evolution close to thimbles the sign problem should be weakened and parameter regions that have been inaccessible so far may be reached. In the following discussion we consider again one-dimensional integrals. One rather general ansatz for such variable transformation is the Möbius transformation This rather general ansatz has four τ -dependent parameters that have to be determined during the simulation. This turns out to be a rather challenging task. We find that the transformation (C1) seems to introduce repulsive structures destabilizing the evolution. Hence, we focus on a special Möbius transformation, namely a rotation where u takes the role of the (complex) field variable and θ is a τ -dependent parameter. Consider a point in the complex plane sufficiently close to the thimble. Then, a rotation suffices to map this point even closer to or onto the thimble. This indicates that the transformation (C2) is both necessary and sufficient for fulfilling the constraints mentioned above. For the remaining part of this paper we always refer to the rotation (C2) when discussing variable transformations.

Lefschetz cooling
Thimbles are curves passing through the fixed points and along them the imaginary part of the action is constant. This gives rise to various constraints which we impose onto the complex Langevin evolution by including the variable transformation (C2).
Let u ∈ C and θ ∈ R. The transformed action reads The procedure here is to be understood as a passive transformation, see the appendix in [38]. Hence, the complex Langevin equation in the transformed variables becomes where η ∈ R.
In the following we investigate how the τ -dependent transformation parameter θ evolves under the dynamics induced by different constraints.
a. Lefschetz cooling the transformed thimbles First, we formulate the additional constraint completely in the transformed theory [35], i.e. we demand This constrains the evolution close to the thimbles in the transformed theory. By taking the total τ derivative of this equation we obtain This rotates the thimble in the original theory precisely onto the real axis of the transformed theory. Next, we investigate how θ changes with the flow induced by (C12). Hence, we take the total τ derivative of (C12), yielding This leads to Im ∂ 2 S u ∂u 2u + Im from which we finḋ upon inserting the drift term foru. In this form, the dynamics are unstable. For instance in the Gaußian model, where the (anti-)thimble is a straight line this manifests itself by a diverging evolution along the anti-thimble to infinity. (C15) rotates the thimble onto the real axis and transforms the anti-thimble into the imaginary axis. Inserting the action of the Gaußian model into (C18) we findθ = 1 σ r sin(θ σ + 2θ) .  Figure 9. Scatter plot of the combined complex Langevin evolution and thimble constraint dynamics in θ with the reversed sign mapping, see (C20) and the discussion in the text. The proposed method is being applied to the z 4 model and the evolution is represented in the original variable z. The scatter plot is color coded, where black corresponds to low density and yellow corresponds to high density. The mapping enforces the sampling of all relevant thimbles, as well as allowing for transitions between them and guarantees stability.
Examining the numerical solution to the previous equation we easily see that the thimble is repulsive whereas the anti-thimble is attractive, see Fig. 8 for an illustration represented in the transformed theory.
To render the thimbles attractive we consider again (C17). Replacingu by the Langevin drift with a reversed sign +∂S u /∂u reverses the sign in (C18) yieldinġ This also enables and enforces hopping between the thimbles and inverts the stability properties of the fixed points. Note that the evolution in u is still governed by the complex Langevin equation C4). Therefore, inserting the sign-reversed drift for u in (C18) can be understood in the sense of a mapping with the following properties: it guarantees that the evolution stays close to the thimbles, as well as allowing for transitions between different contributing thimbles. This approach yields the correct result for the Gaußian model. However, once the thimble structure becomes slightly more complicated, the values of observables are not computed correctly. This has been tested for different actions in [34,36]. While the procedure samples all thimbles, it does not correctly take into account their relative weights. Fig. 9 shows a scatter plot of the τ -evolution applied to the z 4 model. Clearly, this algorithm samples all thimbles. But the contributing ones are being sampled with a higher weight, see the yellow regions in the scatter plot.