A simple approach towards the sign problem using path optimisation

We suggest an approach for simulating theories with a sign problem that relies on optimisation of complex integration contours that are not restricted to lie along Lefschetz thimbles. To that end we consider the toy model of a one-dimensional Bose gas with chemical potential. We identify the main contribution to the sign problem in this case as coming from a nearest neighbour interaction and approximately cancel it by an explicit deformation of the integration contour. We extend the obtained expressions to more general ones, depending on a small set of parameters. We find the optimal values of these parameters on a small lattice and study their range of validity. We also identify precursors for the onset of the sign problem. A fast method of evaluating the Jacobian related to the contour deformation is proposed and its numerical stability is examined. For a particular choice of lattice parameters, we find that our approach increases the lattice size at which the sign problem becomes serious from $L \approx 32$ to $L \approx 700$. The efficient evaluation of the Jacobian ($O(L)$ for a sweep) results in running times that are of the order of a few minutes on a standard laptop.


Introduction
The sign problem poses a challenge for a lattice simulation of many physical theories, ranging from widely-studied systems such as QCD and other theories with a finite chemical potential [1,2], through the simulation [3] of PT-symmetric theories [4,5], to the simulations [6] of more exotic systems such as string field theory [7]. It also appears in the simulations of various condensed matter theories as well as in other systems. In general, one faces a challenge whenever e −S is not a positive quantity. In such a case e −S cannot be identified as a (non-normalized) probability density and the Metropolis algorithm [8] cannot be applied without a modification. A simple modification is reweighting. However, for fast oscillations of the phase of the action the computational cost of evaluation using reweighting grows exponentially with the lattice volume. The sign problem is this exponential behaviour.
Several methods have been introduced in order to circumvent the sign problem: One could attempt to expand the problematic part of the integrand in a Taylor expansion such that all terms in the expansion involve only integrals with a constant phase [9]. While this method works for some systems, it introduces an error in the estimate of observables on top of the standard statistical error. Also, one could not expect it to be useful for large values of the parameters used in the expansion, e.g the chemical potential. Another approach is to analytically continue the problematic parameters to values that lead to no sign problem and then continue back [10,11]. While this method could also be used in some cases, numerical analytical continuation can be quite challenging and the approach suffers from problems similar to those of the Taylor expansion approach. A third approach is the complex Langevin method [12]. Here, stochastic dynamics is used for the calculation of observables and there is no reference to e −S as a probability density. Hence, the sign problem is avoided. Nonetheless, this method also has its limitations [13][14][15].
An important ingredient of the complex Langevin method is the complexification of the degrees of freedom. Complexification can be used in other ways as well: One could attempt to avoid the phase oscillation causing the sign problem by using (the multi-dimensional form of) Cauchy's integral theorem for deforming the original integration contour to one without (or at least with less) phase oscillations. When using Cauchy's theorem, one must take care not to pass any singularities of the integrand and to deform the asymptotic integration range only when it vanishes fast enough, e.g using Jordan's lemma. For most physical theories the analytical continuation of the action is regular. Hence, the first issue is of no concern for these theories. The asymptotic behaviour of the integrand, on the other hand, can become singular for specific contours. For potentials that behave asymptotically as V φ n with N variables (N equals the product of field components by the number of lattice points) there are generically (n − 1) N different homology classes of integration contours [16,17]. One should therefore be careful not to deform the contour away from the original homology class.
A possible prescription for the contour deformation is to use Lefschetz thimbles [17], which also give the steepest descent of the (real part of the) measure. The implementation of Lefschetz thimbles in lattice simulations was introduced in [18]. Despite the success of this approach, it is also not without faults: 1. While it is known that Lefschetz thimbles are manifolds, explicit expressions defining these manifolds are absent. This leads to complicated and expensive algorithms for verifying that the integration contour does not leave the thimble.
2. While in some cases only a single thimble contribution is relevant in the continuum limit [18], in other cases one needs many thimbles. Since the phase can differ among different thimbles, this could lead to reemergence of the sign problem as a "global sign problem", especially in light of the fact that the number of homology classes (the number of independent thimbles) goes to infinity in the continuum limit.
3. Lefschetz thimbles are constructed in a way that keeps the imaginary part of the classical action constant. However, when working with thimbles the integration measure changes. This leads to the "residual sign problem".
The above issues can be summarized by noting that Lefschetz thimbles improve, but do not completely solve the sign problem, while reducing the running cost from exponential to O(V 4 ). In some cases the computational cost of a thimble simulation can be further decreased [19]. Nonetheless, it is still advisable to look for improvements and alternatives to the Lefschetz thimble method. Indeed, generalizations of the thimble method have been proposed, in which the above mentioned issues are addressed by choosing a single integration contour as an approximation to the sum of all the relevant thimbles [20][21][22].
In this work, we propose a method for a different deformation of the integration contour. We do not attempt to approximate the thimbles. Thimbles are defined by a gradient flow, which defines a real N -dimensional submanifold of the complex N -dimensional space.
On the other hand, the requirement of a vanishing imaginary part puts only one real constraint and hence its solution is a (2N − 1)-dimensional real manifold. We intend to choose an N -dimensional subspace of this manifold, or an approximation to such a subspace, in order to obtain a proper integration cycle. This might appear to be too arbitrary as compared to a cycle defined by a gradient flow. However, there is much arbitrariness also in the definition of a gradient flow: In order to define the gradient operator, a metric should be introduced in the complex N -dimensional space. It is possible to choose this metric to be the flat metric, but this choice is not canonical [17]. The only relevant attribute of this metric is its compatibility with the complex structure, which guarantees that the imaginary part is constant along the flow, i.e. the metric must be Kähler [17]. While not all the manifolds with constant imaginary part and appropriate boundary conditions can be obtained by varying the metric (also, some metric variations lead to a reparametrisation rather than a change of the obtained manifold), some can be obtained this way. This illustrates that even "the thimble" is not a uniquely defined object. Moreover, an integration cycle which is not the thimble can be chosen to be a single contour that already takes into account the integration measure. Hence, with a proper choice of integration contour, both the residual sign problem and the global sign problem could be avoided.
For illustrating our method we consider the specific case of a Bose gas at a finite chemical potential. This model is often used for the purpose of examining methods of dealing with the sign problem, e.g. [18,23]. On the lattice (in lattice units) the action of the model is, Here, u k , v k are respectively the real and imaginary parts of a single complex field, the index k runs over all lattice sites and d is the dimensionality of space-time. While the signature is Euclidean, the "time" direction is special, since the interactions related to it differ from those of the other coordinates due to the presence of the chemical potential µ. Moreover, the terms related to the chemical potential are the ones responsible for the imaginary part of the action. The action enjoys a U (1) symmetry, even for a non-zero chemical potential.
Rescaling the fields by a factor of 2d+m 2 λ the action becomes, where we defined (1.4) For fixed d, this constant ranges from α = 0, which corresponds to the infinite mass limit to α = 1 2d , which corresponds to the zero mass limit. For small values of the chemical potential the vacuum of the theory is in an unbroken phase of the U (1) symmetry, while a spontaneously broken phase is obtained for cosh µ > m 2 2 + 1 . (1.5) In the continuum limit the above equation takes the simpler form µ > m.
For concreteness we concentrate in the rest of the paper on d = 1. In this case the first term in the second line of (1.3) is absent. Assuming periodic boundary conditions and N lattice points one can write As long as µ is small enough, taking the integration contour to be the product of the real u k and v k axes would not lead to a significant phase factor. As µ is increased a larger and larger phase appears along this integration contour and the sign problem appears. To illustrate this (well known) situation we attempted to simulate the system with the standard integration contour on lattices of size 16, 24, 32, ..., 72. For simplicity we set, as we do for all the simulations in this paper, λ = m = µ = 1 and estimated the number of configurations that should be sampled in order to obtain a standard error of S that would be about 2% of S . Since this is only for an illustrative purpose we were content with identifying a number of configurations that leads to a relative error in the range of 1.6% − 2.4%. We find that while for 16 and 24 lattice points we need roughly the same number of configurations, for larger lattices the required number of configurations grows fast. We performed a least squares fit for the logarithm of the number of configurations as function of lattice size for the range 24...72 and found that the number of configurations needed grows as e 0.21L where L is the lattice size. We plot the number of configurations needed together with the fit found in fig. 1. Assuming this exponential behaviour continues, we can estimate that the running time (on a standard laptop) with no contour modification needed for the stated precision for L = 96 would already be of the order of a couple of months. This is significantly longer than the running time with contour deformations studied below, for which the precision is also much better.
In light of the above we look for adequate contour deformations. We write the complexified fields as and look for a parametrisation of the contour of the form  that is, we intend to use the real parts of the variables as the space over which we integrate. This is not the most general form of a deformed contour: We do not allow contours that return to previous values of any of the x k , ξ k . However, as we already stressed, there is very large freedom in choosing the integration contour and our choice of parametrisation does not limit us. Finding an integration contour amounts to solving a single equation in 2N variables. One could get to the conclusion that this is a trivial task. However, the solution must be well defined and free from singularities in the whole space spanned by the x k and ξ k , since otherwise the "contour" would not be continuous and Cauchy's theorem would not hold. Moreover, as stated above, care must be taken to ensure that the asymptotic behaviour of the contour is the desired one. While it is still possible to allow the contour to approach other asymptotic regions before returning and heading to the correct one, this should only be allowed as long as the integral of e −S does not diverge along such paths. How should one choose an integration contour? We suggest to follow an old and important principle of physics: Choose a contour that is as simple as possible but not simpler than that. We should now decide upon our criteria for "simplicity" and upon the general considerations that should guide us in our search after good integration contours: 1. Explicit and functionally simple form: We would like to have expressions that can easily be implemented and whose computational cost is small.
2. Approximate expressions: The imaginary part of the action does not have to vanish exactly as in (1.9). A small variation of the phase that does not lead to a severe sign problem can well be tolerated. When possible, we can trade the accuracy of (1.9) for simplicity of the expressions.
3. Include the Jacobian if needed: In order to prevent the residual sign problem, the contribution to the phase coming from the Jacobian can be included in our expressions ab initio. If the residual phase is small we can follow the previous principle and ignore it, if it simplifies the expressions. 4. Fixed functional form: The functions used for defining the contour (1.8) should have the same (relative) functional dependence on the x k , ξ k , for all the y k , ζ k . Here, "relative" means that the dependence of, say, y 1 on, say, ξ 2 would be the same as that of y 3 on ξ 4 . The motivation behind this restriction in not so much the deep physical principle of homogeneity of space-time as the ease of implementation. In particular, since the lattice size can vary, we would like to have a prescription for the deformation that does not change as the lattice size changes.
5. Locality: We would look for contours for which y k , ζ k could depend on x k , ξ k and maybe on their nearest neighbours. Again, the main motivation is computational rather than the deep physical principle of locality: The evaluation of Jacobians can become significantly cheaper computationally for local expressions. Also, the parallelization of local expressions is much easier and much more efficient. Moreover, for contours that are defined in terms of some free parameters, locality allows us to fix their values on a small lattice and use the same values for larger lattices. Of course, the last issues is related to the fact that the original action is indeed local. Physical locality, which in the case of (1.1) manifests itself as the restriction that interactions are nearest-neighbour ones, also implies that a contour for which y k , ζ k only depend on nearest neighbours could grasp the essential contributions to the sign problem.
6. Symmetry: The theory we consider here enjoys a U (1) symmetry. It is not clear a priori that the best possible contours should respect this symmetry. After all, the theory has a spontaneously broken phase. However, when examining ansätze for the form of the contour, one could examine whether the symmetry is respected for ansätze with a small number of parameters. If this is the case, it could be more efficient to restrict the number of parameters at higher levels in a way that respects the symmetry.
We suggest some possible contours in section 2 and a more general ansatz for the contour in section 3. We then study the Jacobians of the various contours proposed in section 4. Simulation results are presented in section 5 and we end with some discussion in section 6.
Note: While this work was prepared the papers [24][25][26] appeared, in which some similar ideas were presented. Our approach shares with [26] the low computational complexity and the ability to use a relatively small number of parameters, while generalizing the ansatz for the contour in a way that enables taking into account the interaction between nearby lattice points.

Simple integration contours
In order to find an approximation for the contour we write the imaginary part of the action (1.6) in terms of the unknown functions y k ({x m , ξ m }) and ζ k ({x m , ξ m }). If we were to choose we would have obtained Instead, we expand the imaginary part of the action in powers of We assume that the functions we are after can also be expanded this way and set their order zero term to zero and write At first order we obtain the equation (at this level of approximation we set cosh(µ) = 1) Setting this expression to zero gives a single linear equation forỹ k andζ k . Such an equation has many solutions. However, we are interested only in solution that are continuous with proper asymptotic behaviour and for which the dependence of y k and ζ k on x k , x k+1 , ξ k , ξ k+1 is the same for all k. On the other hand, we do not need an exact solution. It could be enough, at least in some range of the parameter space, to minimize as much as possible the absolute value of this expression. We can then look for a solution that is perturbative also with respect to α 1 . We do not intend to actually find a full double expansion in terms of σ and α. Instead, we want to recognize the leading contribution and identify contours that cancel it. These contours could well be enough for various values of the parameters and of L. When they are not enough, they could still be used to guide us towards ansätze that generalize them. Such ansätze are introduced in section 3. The term in the first line of (2.5) is the inhomogeneous part of the equation. We can look for a solution in whichỹ k cancels the contributions of the first inhomogeneous term andζ k cancels the contributions of the second one. We further look for a systematic expression, for which the (lowest order with respect to α of the) expression vanishes even before summation, That is, we set,ỹ Similarly we obtain,ζ Note that this expression respects the U (1) symmetry of the theory: If we define we can write (2.7) as The U (1) symmetry rotates both φ k and ψ k in the standard manner. The expressions (2.7) are not bounded in the limits ξ k+1 → ∞ and x k+1 → ∞ respectively, with all other variables remaining fixed. This means that the asymptotic range covered by them might include also parts of the complexified space unrelated to the original contour. Nonetheless, since these expressions can be obtained by continuously deforming the original integration contour, they would lead to the correct results, as long as they do not pass asymptotic regions over which the integral diverges. One could examine this contour in a simulation. However, one could also modify it in order to obtain expressions that are everywhere bounded. A simple possibility is to use the following expressions, These expressions neither solve (2.6) nor do they respect the U (1) symmetry. However, they do it approximately, which might be enough. Plugging (2.10) into (2.5) we obtain, This approximation for the integration contour already reduces the imaginary part of the action significantly. Moreover, the reduction is largest when the real part of the action is smallest, i.e. near points with large contribution to the functional integral. The reduction in the case of the contour defined by (2.7) is even more dramatic. The imaginary part of the action in this case is Of course, in order to turn (2.7) into bounded expressions, one could use also variants of (2.10). In particular, one could extrapolate between (2.7) and (2.10) by introducing a constant κ and settingỹ Now, for κ = 0 we obtain (2.7), while (2.10) is obtained for κ = 1. We could try to look for the value of κ for which the sign problem is most significantly reduced. Eq. (2.13) is an ansatz for the contour depending on a single parameter. We examine all the contours introduced so far in section 5.

A general ansatz for the integration contour
In order to get a better behaved contour, it might be useful to generalize the expressions obtained so far, by introducing higher-order terms. We could also consider a somewhat less local functional dependence. We therefore consider a more general contour 2 , Here the sum in the numerator runs over all non-negative integers µ, ν, ρ, σ, τ, φ such that 0 ≤ µ + ν + ρ + σ + τ + φ ≤ p, and the sum in denominator runs over all even non-negative integers µ, ν, ρ, σ, τ, φ such that 2 ≤ µ + ν + ρ + σ + τ + φ ≤ q. We define (p, q) to be the order of the contour. The reason for restricting the sum to only even integers in the denominator is to ensure the denominator is always positive so that the contour always remains finite. This also requires that the b µνρστ φ are non-negative. However this may be overly restrictive; it may be sufficient that the contour remains finite for the fields actually reached during a simulation rather than for every possible field. The advantage of the rational function forms (3.1) is that they are easy to differentiate and so the Jacobian can be explicitly written. They also contain (2.7), (2.10), and (2.13) as special cases -note that we have chosen the relationship between the two equations so that the general contour has the same symmetry between u k and v k as in (2.7), (2.10), and (2.13). We can further restrict our ansatz by requiring that no terms will appear in the numerator with powers of any variable that does not appear in the denominator with at least the same power. This will result in contours that are everywhere bounded. A further possible restriction is to impose the U (1) symmetry, that is to consider only expressions that can be written in a way that generalizes (2.8). This will not necessarily lead to the best contours. However, if it turns out that the U (1) symmetry gives approximately the best contours at low contour order, it would make sense to simplify the optimisation procedure by imposing the U (1) as a restriction ab initio, reducing in this way the amount of free parameters needed to be optimised.
It is easy to include higher order terms in this contour simply by increasing the order. Similarly, longer-range terms could be included by adding to the ansatz expressions depending on x k±2 , ξ k±2 etc. However, this comes at a cost; we need to determine the parameters a µνρστ φ , b µνρστ φ . There may be a large number of these; already at order (2,2) there are 34 free parameters. Furthermore, the optimal contour will presumably depend on the chemical potential µ, the interaction strength λ, and on α. Another important restriction is to contours for which the Jacobian can be efficiently evaluated. We discuss this issue in the next section. For now we just mention that this restriction implies that one should consider in (3.1) either only terms with x k+1 and ξ k+1 or terms with x k−1 and ξ k−1 , but not both. We refer in what follows to the use of only x k+1 and ξ k+1 as "forward nearest neighbour".
Choosing the best contour sounds like a daunting task; on a large lattice where the sign problem is severe we would expect that for generic contour parameters a µνρστ φ , b µνρστ φ the mean phase will be so small that it will be statistically indistinguishable from zero. How can we then decide how to change the parameters to obtain a larger mean phase? The solution is to tune the contour on small lattices. Since the action is local, the imaginary part of the action is presumably mostly due to local contributions rather than long-range correlations. Our ansatz is also local. Hence, if we tune the parameters a µνρστ φ , b µνρστ φ to maximise the mean phase on a small volume we will obtain a contour that minimises these contributions. Therefore if we use the same contour on a larger volume the imaginary part of the action should still be small. If needed, one could start on a large lattice with the parameters fixed on a small one and update them, but we expect that this will generically not lead to a significant improvement.

The Jacobian
To do calculations with these contours we must parametrise u k and v k in some convenient way. The most natural thing to do is to use the real parts x k and ξ k . We must then include the Jacobian determinant in the evaluation, For the contour defined by (2.7) the Jacobian matrix takes the form, where we defined For the contour defined by (2.10) we have a slightly more complicated Jacobian matrix, where now we defined Similarly it is straightforward to derive the Jacobian for our general ansatz (3.1). We note, that both (4.2) and (4.4) are almost upper-block-bidiagonal with 2×2 blocks. The block structure comes from the fact that we have two fields in every lattice site and the upper-bidiagonal structure originates from the use of forward nearest neighbour in the definition of the contours. In both cases the Jacobian matrix has the form The only term that does not fit the upper-bidiagonal form is B N (marked in red). This term is there in light of the periodic boundary conditions. If it was not for this term, the Jacobian determinant would be the product of 2×2 determinants. There are several ways to deal with this term. One option would be to choose to impose other boundary conditions. An even simpler solution would be to decide not to deform the contours associated with u N , v N . This will generically enlarge the phase, but in an N -independent way that should not lead to a sign problem. If this resolution is chosen one should verify that observables at k = N do not behave in an anomalous way. While analytically this is obvious, the asymmetric treatment of a single point can in principle induce different numerical errors in this point as compared to the other ones. Nonetheless, it is also easy to evaluate the determinant explicitly without any modifications. If the blocks were just numbers we could have used elementary row operations in order to eliminate B N −1 then B N −2 and so on, obtaining a matrix whose determinant is given by In the current case where we deal with blocks instead of numbers, we can perform "elementary row operations" by left multiplication by "elementary block matrices", whose determinant is unity. This leads to a very similar expression, We reduced the evaluation of the determinant to the evaluation of a product of block matrices. Thus we resolved the complexity related to the evaluation of the determinant. One could still worry that since the determinant includes a product of N +1 terms of which the evaluation of the first involves a product of N matrices it would lead to a complexity of O(N 2 ). This is not the case. The locality of our prescription implies that when we change x k or ξ k we should only update A k−1 , A k , B k−1 , B k . Dividing the product by the previous value of det (A k−1 ) det (A k ) and multiplying by the new value of this factor is of O(1). As for the first term in the product, symmetry of the problem implies that it can be written in any cyclic order, This identity can also be proven explicitly, for example by expressing the determinant in terms of a trace and utilizing the cyclicity property of the trace. Thus, when we initiate the simulation we define We start by updating x 1 after which we change We then update ξ 1 after which we change 3 , We used the cyclicity property (4.10) in order to obtain an expression that is adequate for x 2 . We then continue with x 2 , ξ 2 and so on. All the manipulations described are of O(1) and so do not increase the complexity. The above algorithm might not always work: While the A k are always invertible the B k can be non-invertible for specific contours. Moreover, even when the B k are invertible, they might be close to being singular. This, together with the fact that these are complex matrices, could lead to stability issues for the algorithm. There can be several ways to handle this obstacle: 1. One could use a slow algorithm to evaluate the determinant without using the nearlybidiagonal structure, such as LU decomposition (O(N 3 ) per site so O(N 4 ) per sweep), or evaluate M using its definition (4.11) after every update, which leads to complexity of O(N 2 ). This solution is, of course, not the desired one, since even the O(N 2 ) would be much more restrictive than an O(N ) algorithm.
2. The problem described would emerge only for specific values of the parameters defining the general contour. One could then choose a different set of parameters, for which the B k are generically large and the problem does not occur.
3. It is possible to deform all the contours, but one. Suppose we decide not to deform u 1 , v 1 , this modifies the expression for the Jacobian to the simpler form, While such a choice certainly decreases the phase, we expect that this change would be N -independent and so would not change by much the range of validity of the given contour, i.e. the unmodified and modified contours should begin to suffer from the sign problem at not too different values of N .
4. Since the numerical instability would usually occur for small values of B k , one could expect that for large enough N , the matrix M would become extremely small and the approximation det(J) det (A 1 ) · ... · det (A N ) (4.15) would become exact within the simulation precision. This expression is almost identical to (4.14), but the constant decrease of phase is now absent. Moreover, since in this case (as well as in the previous one) there is no need to evaluate the matrix M , the computational cost is further reduced.
In the next section we examine all these possibilities. The discussion so far is inadequate for the general form for the contour (3.1), which leads to an almost tridiagonal-block-matrix. Similarly to the previous case, one could get a genuine tridiagonal-block-matrix by not deforming the integration contours of x 1 , ξ 1 , x N , ξ N . This would lead to a total computational cost of O(N 2 ). However, as in the discussion above, one could do better than that by deciding to take a slightly less general contour not involving the x k−1 and ξ k−1 terms. Moreover, our simulations suggest that the benefit for taming the sign problem from considering two-sided nearest neighbours is not very high. Therefore, such a restriction on the form of the contours could even enable us to get to a better resolution of the sign problem given a fixed number of free parameters.

Results
As mentioned in the introduction, we examine the method presented in this paper by simulating the system with λ = m = µ = 1. The simulation is performed using a Fortran code for L = 4, 8, 16, 24, 32, ..., 96. In all simulations we take 300,000 configurations. The jackknife method is used for the evaluation of statistical errors. We examine the observables S and u 2 k + v 2 k as a function of lattice size for the different contours suggested. Since we only illustrate the method we do not study other observables. Of course, one could examine in this way also the average density and other observables. The mean phase for the same contours is also studied, since it serves as an indicator to the severity of the sign problem.
We refer to the original contour without any deformations (2.1) as "contour 0". The simplest modified contour (2.7) is referred to as "contour 1", and "contour 2" is defined by (2.10). The contour that interpolates between the last two (2.13) is "contour 3". We optimise the value of κ in (2.13) on a lattice of size L = 8 for obtaining maximal value for the mean phase. This leads to κ = 0.25, which is the value we use in the simulations. The general ansatz (3.1) is referred to as "contour 5". In the previous section we suggested the possibility of using a version of contour 5 that only includes forward nearest neighbours. We refer to this variant as "contour 5f".
Since (3.1) defines a family of contours we should specify the value of the parameters used for defining contour 5. Working at L = 8 at order (1, 2), we find that the optimal (mean-phase maximising) contour lies near with all other contour parameters equal to zero. Notice that the parameters are very close to a U (1) symmetric contour, i.e. b 000002 b 002000 and b 000200 b 200000 . This could guide us when we go to higher (p, q) orders.
Our simulations show that while contours 1,2,3 do not suffer from the numerical instability described at the end of the previous section, contour 5f is affected by it 4 : We evaluated M using its original definition (4.11) at the end of each sweep and compared to the result obtained by successively using the update algorithm (4.13) 5 . In the case of contour 5f, the relative size of the components of M , evaluated in these two ways, differed from unity by more the one part per million (the threshold we used), already at L = 16. We did not find any problems of this sort with any of the other contours.
We present the main results for all the contours including a version of contour 5f, for which the numerical stability issue has already been resolved, in the next subsection, 5.1. We then return to examine the issue of the numerical instability and its possible resolutions in subsection 5.2.

The observables and the mean phase for the various contours
In fig. 2 we plot the action density, that is, the expectation value of the action divided by the lattice size, for the various contours described above. The different contours are horizontally separated to enable distinguishing them, as many of them overlap. Since it is hard to infer the standard errors from the plot they are also displayed separately. It is obvious that one cannot trust contour 0 for L 32. As for the other contours, they all give the same results, which are consistent with the expectation of obtaining a constant action density for large L. Given the deformation of the contours and the fact that the action is complex, we could have obtained general complex results. In fact, all the results are consistent with being real. The largest value for |Im S | L is 0.1 for contour 0 and about 0.004 for contour 3. It is smaller for the other contours. This gives further support to our method.
It seems that there is no significant difference between contours 1,2,3 as far as the standard error of the action is concerned. Initially it goes down, but then, around L = 40, it starts to go up again. The situation with contours 5 and 5f is less clear, since the tendencies are masked by noise. It is hard to infer whether they started going up already or not. In general, there can be several reasons for a dependence of the standard error on lattice size: Since we are not at a critical point, the evaluation at sites separated by more than a finite correlation length is roughly independent. We then expect the standard error of the action to scale like that of L/l corr independent measurements of the same object, that is, like √ L. Thus, if this was the only relevant effect, we would have expected the standard error of the action density to scale like L − 1 2 . However, there are two effects that increase the size of the standard error: the increase of the autocorrelation time and the sign problem. Thus, one could suspect that contours 1,2,3 are closer to developing a  sign problem than contours 5 and 5f, although it is hard to make this statement into a quantitative one based solely on the values of the standard error. In fig. 3 we plot the expectation value of u 2 k + v 2 k as a function of lattice size. The general tendencies identified for the action can be seen also here. The results of contour 0 are slightly better now, presumably since no higher order terms are involved.  Figure 3. Left: u 2 k + v 2 k as a function of lattice size. Different contours are separated horizontally for clarity. Right: The standard error of u 2 k + v 2 k . The various contours are colour marked.
In fig. 4 we plot the mean value of the phase as a function of lattice size. One can see that on this logarithmic scale the mean phase decreases linearly for all contours (the lines plotted are mean square fits of the data), until, for contour 0, it gets to the range, in which we already found out that the sign problem kicks in. This loss of linear behaviour (on a logarithmic scale) that signifies the sign problem happens around a mean phase of about 0.01. It is natural to anticipate that, at least for a given theory, similar values of the mean phase would be associated with the sign problem for different integration contours. Thus, by extrapolating the plotted lines to 0.01 one can attempt to predict the onset of the sign problem for the other contours. We therefore expect that the onset of the sign problem would occur around L = 200 for contour 1,2,3, around L = 700 for contour 5f, and around L = 1100 for contour 5.  In order to examine the suggested rule of thumb, we simulated contour 1,2,3 in the range L = 150, 200, ..., 500. The results are shown in fig. 5 together with the previously found fit. The exponential decrease of the mean phase continues until about L = 200, where it begins to fluctuate. In some cases the mean phase becomes negative and so it cannot be presented on the logarithmic scale used. Around the same value of L the action density begins to fluctuate and to develop large standard errors. Yet another (unplotted) consequence of the emergence of the sign problem is the appearance of large values of the imaginary part of the observables. All that illustrates our prediction for this case. Next, we examine the same predictions for contour 5f (we do not examine contour 5, for which we do not have a fast algorithm). The results are presented in fig. 6. The behaviour is essentially the same and we conclude that the rule of thumb suggested holds in this case as well. Such estimates can be very useful, since they can enable one to understand which (p, q) order to pick for obtaining reliable results on a given large lattice, by simulating other contours on much smaller lattices.

Resolving the numerical instability
As mentioned above, contour 5f suffers from numerical instability when simulated using (4.13). Four possible resolutions to the stability problem where offered at the end of the previous section. To these options we add "option 0", of completely ignoring the problem. We rewrite our options below and refer again to the end of section 4 for details: 0. Do nothing.
1. Use LU decomposition; this is exact and numerically robust but slow.
2. Change the parameters to avoid the problem.
3. Leave the contour of one specific point undeformed.
4. Use the approximate form for the Jacobian (4.15).
We examined the same observables as above for all these options. For option 2, we looked for other local maxima (in the space of parameters) of the mean phase on a small lattice. In light of the form obtained for the parameters defining the original contour 5f we considered only the 4-parameter family of U (1)-invariant forward nearest neighbour contours (in fact in this case there are 5 parameters at our disposal, but we found that there is no benefit from allowing a non-zero value for a 000000 ). One of the local maxima we found did not suffer from numerical instability. The parameters of this contour are We use these parameters for defining option 2 in what follows. For option 3, we verified that no anomalous behaviour is obtained at the special point. In fig. 7 the results for the action obtained by the different options are compared. We see that the "approximate" option 4, based on (4.15), is in fact exact beyond a given small critical value L c . Hence, in the previous subsection, we used the exact slow algorithm up to L c and the "approximate" one from this value on. One can still wonder why is the approximation so good and whether this would necessarily remain the case for all values of L > L c . It turns out that, within our working precision, this approximation is exact throughout this range. The reason for that can be traced to the values that can be obtained by the matrices S k (4.9), from which the matrix M is constructed. It is hard to evaluate these values analytically, but we found out that throughout the run, the norm of each row or column vector, within the S k , never gets values above 0.631. Then, using Cauchy-Schwarz inequality, the absolute value of all entries of S k S k+1 cannot exceed 0.631 2 = 0.398. Then, the entries of S k S k+1 S k+2 S k+3 cannot exceed 2 × 0.398 2 = 0.317. Continuing this way, one can deduce that as L grows the possible values of the entries of M become smaller and smaller, until for some critical value L c it gets below our precision threshold. Further increasing L can only improve the approximation. Note that the numbers written above can only be given for fixing an upper bound on L c . From the simulations we actually find that it is lower than this upper bound.
While using option 4 resolves the problem in the case at hand, it is still worth comparing the mean phase in all the cases, since in the general case the actual numbers might be different and it is important to know how the other options could be of help. In fig. 8 the mean phase as a function of lattice size is presented for the various options. We see that the phase of option 0 does not decay exponentially in the range in which this is expected to happen, i.e. the fit in this case is not particularly convincing. This is another The various options are colour marked. Ignoring the problem (red) leads to an erroneous increase of the action beyond L 40. All other options lead to similar and consistent results. As can be seen in the inset, in which the y coordinate is magnified, for L > L c ≡ 24 option 4 (using the approximate expression for the Jacobian) produces similar results to the ones obtained by option 1 (using the exact algorithm). In fact, the two options produce in this range identical results, within our working precision.
indication that something is wrong with this approach. We use again the rule of thumb for predicting when the sign problem would kick in for the various other cases. Extrapolating the numerical fits drawn in fig. 8 suggest that while options 1 and 4 (the "contour 5f" of the previous subsection) are expected to give reliable results up to about L = 700, option 2 should be reliable up to about L = 400, while option 3, which appears to be much lower than all the other ones in the plot, would give consistent results up to about L = 500. The reason for this relative success of this approach is that, as expected, the slope in this case is almost identical to the slope of the usual contour 5f (−0.00645 for contour 5f as compared to −0.00650 for option 3 on the logarithmic scale). The difference between the two cases is related to the large constant phase obtained from the one point whose contour was not deformed. We deduce that in cases in which one cannot claim, as we do here, that the approximation (4.15) is exact, treating one point differently from the others would probably be the best approach. Moreover, one can easily increase the initial low phase in this case by allowing a deformation of the contour in the special point that would be strictly local (that is, without a dependence on any nearest neighbour).

Discussion
In this work we illustrated, using a toy model, a viable method for addressing the sign problem. The method relies on several simple principles, sketched in the introduction (section 1). We demonstrated that with a relatively low number of configurations (300,000) and with a small number of free parameters, one can reliably simulate lattices of not too small size. We also suggested how to generalize our construction in a way that should give reliable lattice simulations also on quite large lattices. While the generalization to other physical systems would require repeating the steps used in the paper in order to construct an adequate ansatz, this construction is relatively simple in principle and presumably can be performed for various systems. The parameters describing the suggested contours depend on the values of the lattice parameters, but their small number together with the ability to find them on a small lattice, suggest that they can be easily found at the beginning of a simulation. Our method leads to a significant improvement in running times. While we estimated in the introduction that at L = 96 running time for contour 0 that suffers from a sign problem would be at least a couple of month on a standard laptop, our actual running time was much shorter and led to more accurate results: Without imposing the fast evaluation of the Jacobian, described in section 4, running time for L = 96 was about one week, while using these results it reduced to about two minutes! Results at L = 1000 were obtained within about 20 minutes. Longer running times on fast computers would enable obtaining results for very large lattices.
While the introduction of a fast algorithm for simulating forward nearest neighbours certainly helps, one could still worry about the fact that only one nearest neighbour could be included. This is not quite so. One could extend our expressions and allow the y k and ζ k to depend on the fields at k, k + 1, k + 2, ..., etc. This would lead to almost uppertriangular matrices. The evaluation of the determinant could then be performed either by setting some of the contours to the trivial ones, or by generalizing the ideas presented in section 4. In this way the approach can be extended to work also for terms induced by the secondary interaction of next to nearest neighbours, etc.
Another possible reservation would concern the generalization to higher dimensions: It seems that when lattice points are coupled to neighbours from different space-time directions one would not be able to obtain the almost upper-block-triangular matrix needed for the fast evaluation of the Jacobian. Again, there are several ways out: 1. At any number of dimensions the 0 direction is special, since this is the direction related to the chemical potential. Including only nearest neighbour terms in this direction would already take care of the most significant contribution to the sign problem.
2. One could generalize the prescription of options 3 of subsection 5.2, such that there would be no deformations of the contour for lattice points that lie along specific co-dimension 1 hyper-surfaces (one could allow some restricted deformations of the contour at these points). This would give an upper-block-triangular matrix 3. One could split the lattice points such that some points would be coupled to neighbours in the 0 direction and some others to neighbours in another directions, in a way that would still produce an almost upper-block-triangular matrix. The approximate continuity of the field on the lattice (enforced by the kinetic term) would induce an effective coupling of all fields in all directions.
4. One could generalize the method introduced here such that some non-locality would be allowed, but in a way that could still be evaluated using a fast algorithm.
We currently examine variants of all these possibilities and hope to describe the results in a future work.