A complex path around the sign problem

We review recent attempts at dealing with the sign problem in Monte Carlo calculations by deforming the region of integration in the path integral from real to complex fields. We discuss the theoretical foundations, the algorithmic issues and present some results for low dimensional field theories in both imaginary and real time.


Introduction
The Monte Carlo approach to field theories is based on the path integral representation of the expectation of operators given by the path integral: where φ stands for the fields in the theory and S is the action of the theory. After discretization the integral in eq.(1) becomes a finite dimensional integral. The numerical computation of an integral of this type is nearly impossible due to the highly oscillatory integrand. The standard procedure to deal with this problem is to use a Wick rotation of the time coordinate t → it which transforms the integral into where now S stands for the euclidean action. In many cases of interest, the euclidean action S is a real quantity. In this case the integral in eq.(2) has the form of a classical thermodynamic average and it can be evaluated by Monte Carlo methods using importance sampling. O is estimated as real and imaginary parts of the action as S = S R + iS I and considering the phase e −iS I as part of the observable with the help of the equation: This method ("reweigthing") works to the extend that the average phase e −iS I S R is not small and the expectation value does not result from a ratio of very small numbers. In the thermodynamic limit (and at low temperature), where the spacetime volume βV diverges, the fluctuations of S I also diverges and the average sign decreases exponentially as βV grows. This is the sign problem. Physical systems with sign problems include some of the most interesting physics problems. Systems at finite density typically suffer from the sign problem. That includes the study of cold, dense hadronic/quark matter, a central problem in nuclear physics and astrophysics. Many models used in condensed matter models, like the Hubbard model away from half-filling, also have a sign problem. In addition, the direct evaluation of time-dependent quantities, even in thermal equilibrium, based on real time integrals like in eq.(1), lead to integrals like the one in eq.(1) and have a particularly nasty sign problem 1 . Given the importance of all these physical systems it is no surprise that many ideas have been put forward to solve or ameliorate the sign problem, too many to even mention (some reviews can be found in [1][2][3][4]). There is even an argument showing that the sign problem is a NPcomplete problem in at least one specific model, implying that a general solution to the sign problem is extremely unlikely to be found [5].
For one-dimensional integrals, there is a well know way of dealing with the oscillatory integrals leading to the sign problem: the stationary phase of steepest descent method. In the stationary phase method one changes the contour of integration in the integral dz e −S from the real axis to a contour on the real plane going through one (or perhaps more) critical points (defined by the condition S (z) = 0) such that the phase of the integrand is constant. Since the change of contour, under some circumstances, does not alter the value of the integral, this change of contour provides a way of computing the original integral without having to deal with oscillatory phases. The stationary phase contour has one extra feature: it goes along the direction where the absolute value of the integrand decreases the fastest, hence the alternative name "steepest descent" for this contour. This makes the stationary phase/steepest descent contour also ideal for semi-classsical evaluations of the integral and it is in this context that it is better known among physicists. For the resolution of the sign problem however, it is the constant phase of the integrand on this contour that is relevant.
In this proceeding we will discuss a recent approach to the solution to the sign problem that can be viewed as an extension of the stationary phase/steepest descent contour method, but the method does not involve any semi-classical expansion. It is based on deforming the manifold of integration from real variables R N to another N-dimensional manifold embedded on the complexified field space C N . There are at this point in time a variety of choices for the precise deformation of the manifold of integration to be used. Even larger is the choice of algorithmic ideas suggest to accomplish the task in practice. We will discuss some mathematical generalities related to the change of manifolds of integration in section 2, one algorithm (and its variations) to compute the path integral in practice in section 3 and present two examples of the method in action in section 4. In the final section we briefly comment on newer ideas and future prospects.

Complexified fields, thimbles, anti-holomorphic flow
Consider N-dimensional integrals of the form encountered in lattice field theory: By a generalization of the Cauchy theorem familiar from complex analysis, the same integrals can be obtained by integrating over an N-dimensional submanifold M of C N parametrized by where φ i is now a complex variable, ζ j are the real variables parametrizing the manifold M, is the analytic continuation of the action and detJ = det ∂φ i ∂ζ j is the jacobian related to the change of variables from φ i to ζ j . Notice that φ i are complex and so is the jacobian J.
The change of domain of integration from R N to M does not necessarily lead to the same value of the integrals and care must be taken when performing this step. In the theory of complex functions of one variable the rule determining which deformations do not change the value of the integral are well known: the deformations allowed are the ones that can be made continuously from the initial to the final contour without crossing any singularity (pole, branch cut) of the integrand including "poles at infinity". When poles are crossed during the deformation the value of the integrand jumps discontinuously. For contours extending towards the infinity the integral is not even necessarily well defined. In fact, the integral exists only for contours that approach the point at infinity along certain directions. For the other directions, the integrand does not decay to zero fast enough and the integral diverges. This can be exemplified with the help of the integral ∞ −∞ dz e −z 4 . This integral only exists for contours beginning and ending along the shaded regions on fig.2. Since the integrand has no singularity at finite values of z we are free to move the contour around as long as the its beginning and end of the contour do not jump from one of the shaded regions to another. Doing that amounts to crossing the singularity at z = ∞. Thus, the real line is equivalent to the contour labelled 1 but differs from the one labelled 2. The integral is not even defined on the contour 3. Notice that in order to deform contour 1 into contour 2, the end of the contour has to necessarily pass through a region where the integrral diverges. Two contour are said to be equivalent if they lead to the same value of the integral. The set of all contours (in which the integral is well defined) can thus be divided up in equivalent classes (homology classes) dependent only on which of the shaded regions the initial and final points of the contour lie.
The integrands encountered in lattice field theory do not typically have any singularities. For the most part they are exponential of polynomials (the Boltzmann factor e −S ), perhaps multiplied by another polynomial (observables and/or fermion determinants) 2 . In theories where the field variables are compact (before complexification), the absence of singularities on the integrand is enough to guarantee that any continuous deformation of the domain of integration does not alter the value of the integral. In theories where the fields take values on a non-compact region, like R N there may be singularities at asymptotically large values of the fields. Just like one-dimensional integrals, the integrand is well defined only if the asymptotic values are approached from certain "good" directions. And like one dimensinonal integrals, the set of "good" directions are grouped in discrete sets (homology classes) and the integral takes on a different value in each of these groups [6][7][8][9][10]. As the domain of integration is changed in such a way to make the asymptotic regions go from one "good direction" to another, the value of the integral jumps discontinuously (a "pole at the infinity" was crossed). In between two "good regions" the integral diverges. Thus, a sufficient condition for a deformation not to change the value of the integral is that it can be continuously connected with the initial domain of integration and that, at all intermediate steps of the deformation, the integral is well defined. This conditions guarantee that the domains lie in the same homology class and can be used as a integration domain instead of the original R N . Now that we know that the domain of integration can be moved to a different one in complex space the question is whether this is advantageous. For instance, we may look for a N−dimensional manifold M embedded in C N on which the phase e −iS I is constant. In fact, the condition S I = constant is one single condition on a function of 2N real variables, defining a 2N − 1 submanifold of C N = R 2N ; there are certainly many N dimensional manifolds in that 2n − 1 dimensional subspace. But not all manifolds M where e −iS I is constant is suitable. For instance, if the Boltzmann factor e −S R defines a complicated multimodal landscape on M such a manifold will be unsuitable for Monte Carlo calculations despite the absence of a sign problem. We will first show, as a proof of principle, that manifolds with constant phase and well behaved probability landscapes exist. This construction can be thought out as a multidimensional generalization of the steepest descent/stationary phase contour of one-dimensional integrals.

Anti-holomorphic flow and thimbles
Consider the anti-holomorphic flow defined by the equations 3 where the bar denotes complex conjugation. The important property of the flow is that it increases the real part of the action while keeping the imaginary part constant. In fact, This property can also be seen by splitting the flow in its real and imaginary part we can see that it is the gradient flow for the real part of the action S R and the hamiltonian flow for the "hamiltonian" S I . The flow has some critical (fixed) points φ 0 where ∂S /∂φ i = 0, which we will assume to be isolated. We can understand the behavior of the flow around a critical point by looking at its linearized form: where ∆φ = φ − φ 0 is the displacement from the critical point and where the vectors e a and its "eigenvalues" satisfy: with real λ a (other solution exist with complex λ a ). It can be shown that • λ a , e a come in pairs: (λ a , e a ) and (−λ a , ie a ) . We will assume that critical point are non-degenerate, that is, all eigenvalues of H are non zero.
• e a form an orthogonal set: 1 2 (e a i e b i + e a i e b i ) = δ ab . • any complex vector in C N can be written as a real linear combination of the e a as v = 2N a=1 A a e a or as a complex linear combination of the e a with positive eigenvalues as v = N a=1 (A a + iB a )e a . From eq.(10) we see that the flow around a critical point is repulsive along N of the 2N (real) directions, the ones given by the eigenvectors e a with positive λ a . The flow is attractive along the other N directions. The points reached by the flow along the repulsive directions form the thimble of the critical point. It follows from the properties of the flow that the imaginary part of the action is constant on the thimble and that the real part of the action grows monotonically as the flow moves away from the critical point. As such, thimbles are the proper generalization of the steepest descent/stationary phase contour to the many dimensional case.
In general, there are many critical points and corresponding thimbles and the original integral over real fields equals to the integral over a certain combination of thimbles. There is a recipe to determine whether a particular thimble is part of this combination: a particular thimble is part of the combination of thimbles equivalent to R N is, by evolution of that thimble by the reverse flow intersects the original integration domain (R N ). An intuitive argument for this statement will be presented on the next subsection.
The idea of using thimbles as the integration region in order to solve the sign problem in path integrals was put forward in [11]. There are two main problems in transforming this idea into a practical method: 1. There is no local characterization of the points belonging to a thimble so there is no obvious, computationally cheap way of making Monte Carlo proposals that lie on the thimble.
2. Except for the simplest toy models it is impossible to find all critical points, their thimbles and to determine which thimbles contribute to the original integral.
Problem number 1 has been addressed with several algorithmic ideas. One is the solution of the Langevin equation on the curved thimble, with an appropriate way of projecting the noise term on the tangent plane to the thimble [12]. The projection is the computationally most expensive part of the algorithm. A second method is an adaptation of the hybrid Monte Carlo algorithm to sampling on a thimble [13]. Here again the difficult and expensive part is the projection that guarantees the Monte Carlo chain stays on the thimble. The third method, which we will describe in more detail on the next section, uses a map of the thimble (in fact, all the relevant thimbles) by points on R N . An effective action on the space of real fields is then used with some standard method, like the Metropolis-Hastings algorithm. This last approach has the advantage of being able, in principle, of computing the integral over all the relevant thimbles at once, even if the existence and location of these thimbles is not known a priori.
Regarding the problem 2 above, it has been argued [11] using universality arguments that the integral over a single thimble is likely to agree with the full result in the continuum limit. Also, the contribution of every thimble is suppressed by a Boltzmann factor e −S 0 R , where S 0 is the the action at the critical point, and one can imagine that at weak coupling or for large spacetime volumes the integral will be dominated by a single thimble. It is, however, very hard to put any of these arguments on a firm footing. The contribution of a thimble depends not only on e −S R but also on its "entropy". As the spacetime volume increases or the lattice spacing decreases the position and number of thimbles changes. Finally, there is a phase e −iS 0 I multiplying the contribution of each thimble; the contribution of several thimbles may cancel or, if many thimbles contribute significantly, a sign problem might reappear when combining their contributions. The plausibility of the one thimble dominance can also be gauged against some of its consequences. In QCD, for instance, an instanton configuration is a solution of the classical equations of motion and so there is a thimble associated with it. If the thimble associated with the trivial configuration dominates the partition function, the instanton contributions and its physical consequences have to disappear in the infinite volume/continuum limits. Fortunately, there are ways to proceed without making any assumptions about the dominance of one or more thimbles and, in fact, those arguments can be tested with numerical calculations. The few examples where the thimble structure was fully understood do not have a continuum or thermodynamics limit where these questions ca be addressed [14][15][16][17][18].

Flowing towards the right thimbles
Consider taking every point ζ of R N and evolving them by the flow in eq.(7) by a fixed time T . The end points φ T (ζ) of the trajectories form a N− dimensional manifold we call M T . Assuming the integrand has no singularities and the flow by time T is well defined for all initial conditions with finite values of the fields, the multidimensional Cauchy theorem applies and the equality of the integrals hinges on the singularities "at the infinity", that is, whether the two manifolds are in the same homology class. In many field theories of physical interest, like gauge theories and σ−models, the field variables are compact, the original domain of integration is not R N but some other N− dimensional compact manifold M 0 and there are not asymptotic directions where the integral might diverge. The integral over M T thus equals the one over M 0 . In the case of non-compact field variables we note that the flow monotonically increases the real part of the action and, consequently, decreases the weight e −S R . Since the integral was well defined over R N and the integrand only decreases as T increases, it is reasonable to assume that the integral remains finite and well defined for all intermediate values of T . Since the value of the integral can only change, as the manifold on integration changes, after going through a region where the integral diverges (at large values of the field), we conclude that the integral will not change in the non-compact case either.
We can now find out what happens to the manifold M T as T is increased. One point ζ 0 will approach the critical point φ 0 asymptotically (if there isn't such a ζ 0 that particular thimble is not part of the set contributing to the original integral). The surrounding points will flow towards points with larger and larger values of S R . However, no point can cross a thimble since the flow is tangent to them. They can only asymptotically approach the thimbles. As a result, the flow "automatically" maps the real fields into a manifold that is close to the right combination of thimbles equivalent to the integration over R N . As T is increased, M T becomes closer to the thimbles. This situation is pictured in fig.2.
There is a jacobian involved in the parametrization of M T in terms of the real fields (the initial conditions of the flow). If the result of flowing the point ζ i by a time T is the (complex) point φ i (ζ), the jacobian J = det(∂φ i (ζ)/∂ζ j ) can be computed by taking a basis (v (a) ) of R N and evolving each of its vectors by the flow dv (a) The computation of the jacobian is extremely expensive. Not only N sets of N computed differential equations need to be solved but the determinant of a N ×N dense complex matrix need to be computed, leading to a computational cost of order N 34 . Also, the jacobian is, in general, a complex number: J = |J|e iα . A Monte Carlo evaluation of an integral over M T requires then the reweighting of the "residual phase" α as well as e −iS I : At large values of T , the manifold M T is close to the (correct set of) thimbles and, since S I is a constant over each thimble we don't expect S I to fluctuate much unless many thimbles with different phases contribute almost equally to the integral. Experience shows also that for the variety of models and parameter space studied up to now the residual phase also fluctuates little and there is no practical problem in reweighting it. For weakly coupled theories (not necessarily perturbative) this is expected. In fact, free theories, which have a quadratic action, have flat thimbles with constant α. The experience accumulated investigating several models over a large swatch of parameters space suggests that the fluctuation of α and, for large enough flow time T , the fluctuations of S I also, are small and easily reweigthed. Of course, this always has to be verified in a particular calculation. An important theoretical question is whether this experience will extend to the continuum and infinite volume limits in 4 dimensional theories. where φ T (ζ) is obtained by integrating the flow equations eq.(7) with ζ as the initial conditions and solving eq.(9).
3. Accept or reject the new configuration ζ based on the real part of difference of effective actions between points ζ and ζ with the usual Metropolis probabilities.
4. Return to 2 until a suitably large number of (statistically independent) configurations is collected. 5. Compute the average of the phase e −iS I +iα and the observable Oe −iS I +α over the collected configurations using the imaginary parts of the action (ImS [φ]) and the jacobian (α = Im log detJ).
An important observation is that any value of the flowing time T leads to a manifold M T equivalent to R N . If T is too small, however, the phase e −iS I will fluctuate a lot and the uncertainty on the ratio Oe −iS I +iα / e −iS I +iα will make the method unpractical. In this case a larger value of T is required. There is, however, a problem that can arise at too large a value of T . At large T , a small region of R N is mapped on to a large region of the thimble. Suppose two different thimbles contribute significantly to the integral. Since every thimble is mapped by one small region of R N , so the region of R N contributing significantly to the integral will split up into two small, isolated regions. In other words, the probability distribution ∼ e −Re(S e f f ) is multimodal, which is a situation well known to be challenging for Monte Carlo evaluations. This poses a serious difficulty for the integration over M T mapped by the flow as done here. This will be exemplified later. It would be useful at this point to list some advantages and disadvantages of this algorithm over previously proposed ones. The main advantage is that no previous knowledge about the position or shape of thimbles, or whether they contribute to the integrals, is necessary. The main disadvantages are i) the high computational cost of the jacobian, ii) the multimodality of the distribution in parameter space ζ if T is too large and iii) the fact that isotropic proposals in ζ-space (R N ) result in highly anisotropic proposals in M T which do not match the probability distribution profile in M T and make the sampling inefficient. This last problem is much reduced in many, but not all models, by rescaling the size of the proposals along different directions according to the expectation given by a quadratic approximation of the action [20].
Some extensions of the basic algorithm were proposed to deal with the trapping of Monte Carlo chains due to the multimodality of the distributions in parameters space [22,23]. Some ways of dealing with the high cost of the jacobian are described in the next subsection.
We know that at large enough T , the manifold M T approaches the correct combination of thimbles where the sign problem is minimized (it doesn't disappear completely because more than one thimble, with different phases e −iS I , may contribute). It is reasonable to expect that at intermediate values of T the sign problem will be somewhat ameliorated. Before we delve into algorithmic issues, it would be useful to understand how in practice that occurs.
I may seem a little mysterious how the sign problem can be ameliorated by computing the path integral over M T instead of over R N . After all, the flow preserves the imaginary part of the action and therefore does not change the fluctuating phase e −iS I which is the origin of the sign problem. The apparent paradox is resolved when one realizes that only a small region of R N is mapped into points near a particular thimble. If ζ 0 is mapped to φ T (ζ) ≈ φ 0 , (φ 0 being a critical point of the action) a small neighborhood around ζ is mapped near the thimble attached to φ 0 . As T is increased the size of the neighborhood around ζ 0 mapping the same region near the thimble shrinks. The points of R N that are not mapped to point near thimbles flow towards infinity or other regions of C N where the action diverges (if any exists). The consequence of this picture is that there is little phase (e −iS I ) variation on points near the thimble and all the phase variation found in R N is pushed towards regions of M T with very large values of S R and, consequently, contribute little to the integral. In the limit T → ∞ only an infinitesimal region around ζ 0 is mapped on to the thimble. All remaining points in R N flow to regions where S R diverges (usually the infinity) and do not contribute to the integral. This situation is pictured in fig. 3.

Cheaper jacobians: estimators, Grady-style algorithm
One strategy to deal with the high computational cost of the jacobian is to have a rough estimator of it to be used to generate the Monte Carlo chain. The difference between the estimator and the correct jacobian can be included as a reweigthing factor, as it is done with the phase of the jacobian. We have found two estimators that seem to be useful in (some) of the models we have studied. The first estimator of detJ is W 1 defined by where e (a) is the (complex) basis formed by the solutions to eq.(11) with positive eigenvalues. The hessian ∂ 2 S ∂φ i ∂φ j is to be computed along the flow trajectory starting at ζ while the basis e (a) is computed from the hessian at one critical point of a "dominating" thimble. For quadratic actions (free theories), W 1 = J and, consequently, it is to be expected that W 1 is a good estimator in weakly coupled theories. Experience shows that it is an useful estimator even in theories which, by any measure, are strongly coupled. The second estimator W 2 is even simpler to compute: W 2 is the determinant of the matrix W 2 i j satisfying the equation 10 EPJ Web of Conferences 175, 01020 (2018) https://doi.org/10.1051/epjconf/201817501020 Lattice 2017 satisfied by J. This suggests that W 2 will be a good estimator if J is nearly real, which agrees with our experience. We emphasize, however, that the use of these estimators does not introduce an uncontrolled error in the calculation since, when making measurements, the difference between the estimator and the correct jacobian is reweigthed. A more elegant way of dealing with the cost of the jacobian, which we will call the Grady algorithm, is given by a modification of our algorithm in eq.(3.1) that includes a way of bypassing the computation of the jacobian by adapting the Grady algorithm [24,25] originally created to deal with fermion determinants. The main idea [26] is to make a modified proposal that is isotropic in the variable φ T (ζ), not in ζ, and then correct for that by modifying the accept/reject step. The proposal probability is given by where ∆ζ = ζ − ζ, J = J(ζ) is the jacobian matrix and η is a complex vector tangent to M T . The action of J on a real vector is the same as evolving the vector by the equations eq.(9). The result is, by construction, a vector tangent to M T : J∆ζ = η . However, flowing an imaginary vector, normal to R N , does not give a complex vector normal to M T . Thus, the proposal with probability above can be found by picking a random vector from the distribution ∼ e −η † η , solving iteratively the equation J∆ζ = η and then computing η = JRe(∆ζ). The action of J −1 on a complex vector is obtained by flowing the real and imaginary parts of the vector separately by the inverse of the flow in eq.(9). The backward flow has to be applied to the real and imaginary parts separately as the flow evolution is not a linear operator (it has an anti-linear part). The accept/reject step proceeds as follows. After a proposal ζ is chosen, a random real vector ξ is chosen with probability distribution p[ξ] = e −ξ † J † J ξJ † J √ detJ † J by iteratively solving J ξ = η, where η is chosen from a gaussian distribution ∼ e −η † η . The probability of acceptance of ζ is given by min(1, e −S e f f [ζ ]+S e f f [ζ] ). It is easily shown that this procedure satisfies detailed balance [26].
4 It all comes together in one example:

The 1 + 1D Thirring model
After quite a bit of theoretical discussion it will be good to see how these ideas apply -and how well they fare -in an actual numerical calculations. We choose the 1 + 1 dimensional Thirring model at finite density (and temperature) as a suitable testing ground. It has the advantage of having less degrees of freedom than a 4-dimensional theory would, while sharing some features with QCD like asymptotic freedom, a formulation where the sign problem appears (only) at finite density and a complex fermion determinant (bosonic sign problems were studied using similar approaches in [12,[27][28][29] ). The model is defined in the continuum by the (euclidean) action where the flavor indices take values α, β = 1, . . . , N F , µ is the chemical potential and the Dirac spinors ψ, ψ have two components. Both Wilson and staggered discretizations of the fermions were used in [30]; we will describe here the results with Wilson fermions. In that case the discretized action chosen is with Notice we chose to use auxiliary bosonic variables that are periodic so the original manifold of integration of the path integral is (S 1 ) N . This choice makes the model close to gauge theories and avoids some of the discussion we had in section 2 regarding the asymptotic behavior of M T . The integration over the fermion fields leads to This action describes N F Dirac fermions in the continuum. For µ 0 the determinant det D(A) is not real so this model cannot be simulated by standard Monte Carlo techniques. From here on we choose N F = 2. In the continuum this model is solvable but, since we will work not too close to the continuum. The thimble structure of this model is partially understood, specially in the subspace A 0 (x) = Reφ + iImφ = constant, A 1 (x) = 0 where it is similar to the one of the 0 + 1 dimensional Thirring model [14,31,32]. There are a number of critical points in the A 0 (x) = φ, A 1 (x) = 0 complex plane. There are also points where the action diverges. They arise because the fermion determinant vanishes for some values of the A µ fields. They do no imply any singularity in the integral, they actually correspond to zeros of the integrand and the contribution of the fields around them is suppressed. Since the action approaches S R → ∞ there, they are attractors of the anti-holomorphic flow. The picture that emerges is that these singularities of the action arise at the border of the thimbles and different thimbles approach these singularities from different directions. The manifold defined by detD(A) = 0 has 2N − 2 (real) dimensions. Two manifolds, each with dimensions N, generically meet only at isolated points with dimensions 0. But thimbles are tangent to the anti-holomorphic flow that "seeks" points with divergent action so they end up meeting at a larger dimensional space. This structure is generic to models with fermions and it is depicted on fig.4.
Let us consider some specific values of the parameters for illustrative purposes. With m = −0.25, g = 1.0 we find, by performing standard Monte Carlo at zero density on a 10 × 10 lattice, a fermion mass am f = 0.30(1) and a bosonic mass (for the state generated in the continuum byψ a γ 5 τ 3 ab ψ b ) of am b = 0.44(1). Since m b << 2m f these parameters correspond to a strongly coupled model. An attempt at simulating this model with usual techniques, namely, an integration over (S 1 ) N reveals that there is a terrible sign problem for all chemical potentials larger m f . In fig.5 we show the average sign (right panel) and the fermions density (left panel). It is evident that, even on this small lattice, there is a severe sign problem.
The critical point with the smallest action is the one nearest to (S 1 ) N at φ = φ c = purely imaginary and, at least at weak coupling, is expected to dominate the path integral. The tangent plane to the thimble at φ c is purely real. If we take the domain of integration to be this tangent plane we expect the sign problem to improve as compared to the one on (S 1 ) N since the action, and therefore its imaginary part, varies only quadratically with the deviation from the critical point. In other words, the tangent plane is a good approximation to the thimble for points near the critical point. Since the phase e −iS I is fixed on the thimble, we expect it to vary little on the tangent plane. This expectation is supported by actual calculations shown in fig. 5. Notice that calculations on the tangent plane are just as computationally cheap as the usual integration over real variables. Still, for values of the chemical potential high enough (µ > 3m f ) the sign problem returns, even in these modest lattice sizes. We now use the algorithm described in section 3 except that we start the flow not from real values of the field A µ (x) but from the A 0 (x) = φ c , A 1 (x) = 0 plane in order to take advantage of the improvement obtained on the tangent plane. A cross section of the manifold M T obtained that way with several a N t a N β flow times is shown in fig.4. As discussed above, the higher the value of T the closer M T is to the (correct combination of) thimbles. The integrations over each one of these manifolds M T give the same result but at different computational cost as the sign problem is much milder at larger T . This can be seen on fig.5. The average sign becomes much larger already at T = 0.4, almost completely solving the sign problem and leading to the results on the left panel of 5. There is a price to pay, however, if T is too large. The regions of M T contributing significantly to the integral are mapped by regions of (S 1 ) N that are well separated by barriers of low probability ∼ e −S R . This situation is hard to be dealt with by any Monte Carlo calculation (although "tempered" algorithms have been used to deal with it [22,23]). If not extra care is taken (and the computational resources spent), the Monte Carlo chain becomes trapped in one of the modes of the distribution and only the region of M T close to one of the thimbles is properly sampled. In cases where the contribution of the other thimbles is sizable (given the statistical uncertainty of the measurements), the final result is erroneous. On the other hand the trapping due to excessive flow can be used to separate the contribution of different thimbles to a given observable. In [21], the 0 + 1 dimensional Thirring model was studied and it was shown that the trapped calculation, sampling only one thimble, gives results that differ by small but statistically significant, amounts from the exact analytic result (the single thimble calculation discussed in [33] also agrees with our trapped calculation.). The approach to the thermodynamics and continuum limits in this model does not seem to bring any extra difficulties in this model [21]. In particular, the plateau structure on observables due to shell effects in a finite volume are very visiby and correctly described by the calculations. As pointed out in [15], this structure is washed out if the contribution of only one thimble is taken into account.

Real Time Dynamics
As is well known, a number of observables like particle masses and matrix elements, can be extracted from the imaginary time correlators. Other observables like transport coefficients can be found only from real time correlators. While in thermal equilibrium the real time correlators can be, in principle, obtained from the knowledge of the imaginary time ones, the process requires an analytic continuation from the discrete, imaginary values of energy available in the Matsubara formalism to continuous real energies. Although the analytic continuation is possible (and unique), the analytic continuation is numerically unstable and the statistical uncertainty of a Monte Carlo calculations makes it nearly impossible in practice. An alternative approach would be a direct calculation in real time. A framework for calculations of thermal averages in real time exists, the so-called Schwinger-Keldysh formalism [34,35]. It has been almost exclusively used in conjunction with weak coupling expansions (the exceptions we are aware of are the complex Langevin studies in [36][37][38][39] ).
Consider a system described initially by a density matrix of the form ρ = e −βH /Z 5 and whose evolution in time is governed by the hamiltonian H. We will be interested in observables of the form O 1 (t 1 )O 2 (t 2 ) = tr (ρO 1 (t 1 )O 2 (t 2 )). In the case H = H this is an equilibrium correlation function and it depends on the difference t 1 − t 2 . Fully non-equilibrium situations can be dealt with the same formalism but, for now, we will restrict ourselves to the equilibrium case.
Expectation value of observables are given by: where U(t, t ) = e −iH(t−t ) and T is an arbitrary time T > t. The case of 2-point functions or more is dealt with similarly. This expression describes the propagation from 0 to t, then to T , down to T − iβ/2, then backwards to −iβ/2 and down the along the imaginary direction to −iβ. It has a path integral representation in terms of an action defined on a time contour in the complex t plane depicted in fig.6 6 : The contour action S c is purely imaginary along the part of the contour lying on the real axis. The discretization we used was: where t and n indexes the lattice along the time and spatial directions, a is the spatial lattice spacing and a t is the time lattice spacing: and N t , N β are the number of lattice points on the real and imaginary axis, respectively. The integral in eq.(25) has "perfect" sign problem: the average sign is not small but strictly zero, as can be seen by considering the variation of the contour action when the value of the field on a real value of t is considered. The real part of the action is not changed by this variation and there is no damping of the Boltzmann factor as φ(t) is varied.
The application of the generalized thimble method to real time problems faces some extra difficulties. Consider, for instance, a 1 + 1 dimensional φ 4 theory [26]. The first problem is that the tangent plane to the critical point φ c = 0 is not real and, in some directions in field space, points towards the imaginary direction. Therefore, the tangent space is likely not in the same homology class as R N and  is not a suitable integration domain. The computationally cheap improvement of the sign problem obtained by simply integrating over the tangent space to the thimble is not available in this problem. We instead use the usual flow and perform the integration on the manifold M T defined by a substantial amount of flow. The leads to the problems listed in section 3.1 being exacerbated. We use then the Grady method of section 3.2 and validated by performing calculation at small λ = 0.1 coupling and comparing it with the results from perturbation theory. For large values of the coupling the theory is no longer perturbative although the higher momentum modes still are. We show in fig.7 some of the results for the parameter set N x = 8, N t = 8, N β = 2, m = 1, a = 0.2 and λ = 1. Our calculations were limited by the need to use very large flow times when t is large. This leads to trapping of the Monte Carlo chain. Using a smaller flow time requires prohibitively large statistics. This is a problem that must be solved before the computation of transport coefficients becomes feasible.

Future Prospects
Considering the pace of development in the last couple of years, there is still plenty of room for algorithmic development in the thimble approach to the sign problem. Even the ideas discussed in this paper have not been explored fully; they should be applied to more models, implemented more efficiently, perhaps in more powerful computers (all calculations shown here were performed in runof-the-mill laptops/desktops). Certain features and trade-offs are only understood when larger, more realistic calculations are attempted. Even the models that were already studied were studied on a cursory way with the goal of developing techniques, not understanding the physics (no effort was put into approaching the continuum or large volume limits). In addition, a few novel ideas were recently put forward and are still largely undeveloped.
In order to bypass the computational cost of "flowing" points from R N to M T , it would be convenient to have a analytic parametrization of the manifold M T and the associated jacobian. It is far from obvious how to find such a form. One possibility is to compute φ T (ζ) with a suitably large T for a set of representative values of ζ, interpolate them and then use the interpolating function as the definition of M T . This multi-dimensional interpolation is not an easy feat. In [40] this problem was attacked through the use of machine learning techniques. More precisely, a feed-forward network was setup that takes as an input the real values Reφ T of the flowed field in M T and returns (an approximation