Supersymmetric quantum mechanics on the lattice: III. Simulations and algorithms

In the fermion loop formulation the contributions to the partition function naturally separate into topological equivalence classes with a definite sign. This separation forms the basis for an efficient fermion simulation algorithm using a fluctuating open fermion string. It guarantees sufficient tunnelling between the topological sectors, and hence provides a solution to the fermion sign problem affecting systems with broken supersymmetry. Moreover, the algorithm shows no critical slowing down even in the massless limit and can hence handle the massless Goldstino mode emerging in the supersymmetry broken phase. In this paper -- the third in a series of three -- we present the details of the simulation algorithm and demonstrate its efficiency by means of a few examples.


Introduction
The reformulation of supersymmetric quantum mechanics on the lattice in terms of bosonic and fermionic bonds as derived in the first paper of our series [1] provides a perfect setup for Monte Carlo simulations. First of all, the reduction in complexity by going from continuous to discrete variables is enormous. More specifically though, expressing the Grassmann fields in terms of fermionic bonds avoids the expensive calculation of the fermion determinant and allows the use of special algorithms for which critical slowing down is essentially absent [2,3] and simulations are possible even in the massless limit [4]. This is of particular importance for systems with broken supersymmetry, since the physics of those is driven by the massless Goldstino mode. In the present paper -the last in a series of three -we describe in detail such an algorithm and demonstrate its efficiency. Since the model can be solved exactly at finite lattice spacing by means of transfer matrices, as discussed in the second paper of our series [5], there is in principle no need for numerical simulations. Hence, the present paper rather constitutes a feasibility study to test the practicability and efficiency of the proposed simulation algorithm for the quantum mechanical system in the bond formulation. In that sense it also serves as a preparation for the application of the algorithm, in particular the fermionic part, in more complex situations, such as in supersymmetric Yang-Mills quantum mechanics [6], in the N = 1 Wess-Zumino model [7][8][9] or in the supersymmetric nonlinear O(N ) sigma model [10]. The advantage of the application of the algorithm in the quantum mechanical model presented here is of course the fact that the correctness of the algorithm can be crosschecked with the exact results from the transfer matrix approach, and that the algorithm can hence be validated in detail.
There is another rather pedagogical reason which motivates to consider a new simulation algorithm for quantum mechanics in the bond formulation. Often, simple quantum mechanical systems such as the harmonic and anharmonic oscillator are used to introduce the path integral approach. Similarly, the systems also provide a pedagogical context in which various Monte Carlo simulation algorithms can be illustrated and discussed, see for example [11] for an early example. However, it turns out that the standard Metropolis algorithms and even more advanced algorithms such as the overrelaxation or heat bath algorithm become extremely inefficient towards the continuum limit. This has to do with the usual critical slowing down of the simulations towards that limit, and for the anharmonic oscillator also with the suppressed tunnelling at small lattice spacing. The algorithms presented here do not suffer from these deficiencies, because they eliminate critical slowing down. In addition, in the bond formulation the Z 2 -symmetry φ → −φ is exactly maintained for each bond configuration.
Last but not least, the numerical simulations presented here serve as a test of the practicability of the solution of the fermion sign problem proposed in [4] and discussed further in the first paper of our series [1]. The solution is based on two ingredients. Firstly, the lattice regulates the vanishing Witten index and therefore also the sign problem. Secondly, the fermion loop formulation provides a tool to handle the fluctuating sign, because it naturally separates the contributions to the partition function into topological equivalence classes, each possessing a definite sign. Nevertheless, it is a priori not clear whether the lattice artefacts and the statistical fluctuations can be kept under sufficient control in a practical simulation. The statistical fluctuations of the sign are essentially determined by the amount of tunnelling between the topological sectors, i.e., between the fermionic and bosonic vacuum. In order for the fermion update algorithm to be a true solution to the sign problem, it must guarantee a sufficiently efficient tunnelling rate. The results in this paper demonstrate that this is indeed the case. Not surprisingly, the open fermion string algorithm discussed here has proven to be extremely successful in the N = 1 Wess-Zumino model [9] in which the fermion sign problem is prevailing.
Of course, supersymmetric quantum mechanics has already been simulated on the lattice in various setups using standard algorithms, cf. for example [12][13][14][15][16][17][18][19][20]. However, the bond formulation together with the simulation algorithm presented here brings the numerical nonperturbative calculations to a new, unprecedented level of accuracy. In that sense, the results presented here and partly in [4] serve as a benchmark against which new formulations or simulation algorithms can be tested.
The present paper is organised as follows. In Section 2 we construct in detail an algorithm designed for updating the bosonic and fermionic bond configurations. The discussion includes the explicit update steps and the derivation of the corresponding acceptance ratios. Their evaluation requires the calculation of site weight ratios which turn out to become numerically unstable for large site occupation numbers. Therefore, in Section 3 we present a computational strategy which allows to evaluate the ratios for arbitrarily large occupation numbers. In Section 4, we then present the results obtained using the proposed algorithm. The simulations are for the same discretisation schemes and superpotentials we used in the previous two papers [1,5]. Since this section is merely meant as a validation of the algorithm, the discussion of the physics behind the results is kept short and we refer to the exact results in [5] for a more thorough discussion.

Simulation algorithm
We start our discussion from the partition function of supersymmetric quantum mechanics on the lattice written as a sum over all allowed, possibly constrained bond configurations where the fermion number F = 0, 1 is determined by the fermionic bond configuration {n f (x)} with n f (x) = F, ∀x, and the weight W F (C) of a configuration is given by Here, x denotes the sites of the lattice and i labels the various types of bosonic bonds b i with i ∈ {j → k | j, k ∈ N}. The corresponding bosonic bond weights are denoted by w i and n b i (x) ∈ N 0 is the occupation number of the bond b i connecting the sites x and x + 1. The site weight Q F depends on the site occupation number, i.e., the total number of bosonic bonds connected to site x, and is given by In Section 3 we will discuss in detail the computational strategy necessary to reliably evaluate ratios of these integrals for arbitrary and possibly large site occupation numbers. The type of bonds b i , the weights w i as well as the potential V (φ) and the monomer term M (φ) in eq.(4) depend on the specifics of the chosen discretisation and the superpotential P (φ). We refer to the appendix of our first paper [1] for a compilation of the discretisations and superpotentials considered in our series.
As mentioned above, the bond configurations C = {n b i (x), n f (x)} are possibly constrained. In particular we have the local fermionic constraints while the local bosonic constraints may or may not be present depending on the bosonic symmetries of the system.
The challenge of updating constrained bond configurations lies precisely in the difficulty to maintain the constraints while moving efficiently through the configuration space Z. In [21] Prokof'ev and Svistunov proposed to extend the constrained bosonic bond configuration space by introducing local sources which explicitly violate the constraints. The so-called worm algorithm then probes the extended configuration space by moving the local violations around the lattice, thereby sampling directly the bosonic correlation function corresponding to the sources introduced. The contact with the original configuration space Z is established when the violations annihilate each other, e.g. when moving to the same site on the lattice, such that the bond configuration fulfils again all constraints.
In [2] the idea has been extended to fermionic systems expressed in terms of fermionic bonds. The fermionic constraint in eq.(5) allows only either an empty or a completely filled fermion bond configuration. The difficulty for the direct application of the worm idea to the fermionic system lies in the fact that the introduction of the fermionic source term ψ x ψ x is incompatible with the presence of the fermion loop at site x. A simple solution is to allow the unphysical situation of the site x being occupied by a propagating fermion and two additional sources. Such a configuration violates the Pauli exclusion principle and does not contribute to any physical observable. In the Grassmann path integral such a configurations indeed vanishes trivially.
In order to be more explicit it is necessary to introduce the bond configuration spaces of bosonic and fermionic two-point correlation functions, G b F and G f , respectively, following the notation in our first paper [1]. Bond configurations in G b F contribute to the non-normalised bosonic two-point function according to while the configurations in G f contribute to the non-normalised fermionic two-point function as where F denotes the set of lattice sites belonging to the open fermion string associated with the fermionic correlation function. The key point of the bosonic and fermionic updating algorithm is that the bond configurations for g b F (0), g f (0) and Z F have identical bond elements. As a consequence, statistics for g b,f and Z can be accumulated in the same Monte Carlo process. If the bosonic constraints in eq.(6) are not present, e.g. for superpotentials with broken supersymmetry, the equivalence of bond configurations even extends to g b F (x), i.e., Z F = G b F . The movements from one configuration space to the other are induced by introducing or removing bosonic or fermionic sources according to the scheme given in figure 8 of our first paper [1]. For convenience we reproduce it here in figure 1.
In the following we will now discuss in detail the various updating steps which establish explicitly the connection between the bond configuration spaces G f , Z F , G b F and in addition move the system within G f and G b F . The moves are generated by a Monte Carlo process with probabilities given by the weights of the configurations in eqs.(2), (7) and (8). In particular, we derive the transition probabilities P X (C → C ) for the transition X from bond configuration C to C , which is then accepted by the usual Metropolis prescription P acc (C → C ) = min{1, P X (C → C )}.
(9) In order to simplify the discussion we select the update from Z F to G b F or G f with equal probability which is balanced by corresponding proposal probabilities to select between moving in G f and G b F or returning to Z F .

Updating the fermionic bond configuration
Here we discuss the various update steps which moves the system within G f and relate the bond configurations spaces Z 0 and Z 1 via G f .
Moves within G f are induced by shifting ψ by one lattice spacing from site x to site x + 1, and vice versa, while keeping the other source ψ fixed. Such an update step is graphically illustrated in figure 2 and is called 'shift' update step. A shift in forward direction, x → x + 1, Now let us consider the case when we propose to shift the source ψ forward to the site x where the source ψ is present, as depicted in the upper half of figure 3. x → x − 1 or x − 1 → x, respectively, and the 'put/remove' update step ∅ → x and ∅ → x. The sources are marked with a for ψ and a × for ψ. The bosonic background bond configuration is not drawn.
The forward shift update step x − 1 → x is balanced with the backward shift update step x → x − 1. This backward shift, however, is proposed with probability 1 instead of probability 1/2 since the shift of ψ from x → x + 1 would involve the creation of an open fermion string around the entire lattice. The asymmetry in the proposition probabilities is balanced by the choice of the probability p rm = 1/2 to remove the sources ψψ, such that we find the acceptance ratio for a shift in backward direction x → x − 1 to be the same as in eq.(10), namely The shift step is balanced with the corresponding one in forward direction with acceptance ratio as given in eq.(11), The step from G f to Z 0 and vice versa is induced by introducing or removing a pair of fermionic sources ψψ at site x, respectively. It is called 'put/remove' update step and is graphically illustrated in the lower half of figure 3. The removal of the fermionic sources is suggested with probability p rm = 1/2 and is balanced on one side by the probability to add bosonic sources, and on the other by the probability to shift one of the sources and hence move within C f . Because the 'put/remove' update step does not alter the fermionic bond configuration, we have Z 0 = G f (0). On the other hand it adds or removes a fermionic monomer term M (φ) at site x. The relative weight of the configurations with or without this term is given by Q 1 /Q 0 and the acceptance ratios on a lattice with L t sites become The factor L t compensates for the proposition probability to choose lattice site x out of L t possibilities when putting the sources, while the factor 2 compensates for the asymmetric shift proposal probability when moving ψ from x to x − 1, since the shift of ψ from x to x + 1 is not allowed.
Next we consider the shift update step for the case when the source ψ at site x + 1 is shifted backwards to site x which is already occupied by the sink ψ. The step is graphically illustrated in the upper half of figure 4. While the resulting fermion bond configuration is a valid one (it belongs to Z 1 ), the whole fermion configuration including the source and the sink represents an unphysical situation, and in fact does not contribute to any physical observable, as discussed before. Therefore, such a backward shift from G f to Z 1 , essentially closing the open fermion string, automatically induces the removal of the fermionic source and sink pair ψψ from site x as illustrated in the lower half of figure 4. Such a step is called a hybrid 'shift/remove' update step. Of course, the step is balanced with a hybrid 'put/shift' update step when the additional fermionic sink and source variables are put on a closed fermion loop at the site x. As usual, the acceptance ratios for the hybrid update steps can be read off from the weights of the configurations involved and yield The factor L t compensates for the proposition probability to choose the same lattice site x when putting the sources ψψ back on the lattice, whereas the factor 2 compensates for the proposition probability to shift in forward or backward direction when the fermion string is still open. Note that there are no ratios of Q-weights involved, since no monomer term is added or removed by the hybrid shift/remove update step.
To complete our discussion of the fermionic bond update, we note that the algorithm provides improved estimators for the fermionic two-point function g f (x) and the partition functions Z F . Because the algorithm samples directly the configuration space G f , every open fermion string configuration contributes unity to the stochastic Monte Carlo estimator for g f (x). To be precise we have where x 1 and x 2 are the end and starting point of the open fermion string, i.e., the positions of the sink ψ and the source ψ, respectively. Similarly, every bond configuration in Z F is generated with its proper weight and hence contributes unity to the stochastic estimator for Z F , i.e., the Monte Carlo estimator for Z F is simply Finally we note that the factors of L t appearing in the acceptance ratios above may become inconvenient in practice, especially towards the continuum limit when L t → ∞. The factors only occur when contact between Z F and G f is made, i.e., they are responsible for getting the relative normalisation between Z F and g f right. However, since we make use of translational invariance in eq.(18) the factors of L t are in fact cancelled and can hence be omitted.

Updating the bosonic bond configuration
In this section we now discuss the update steps which relate the bond configuration spaces Z F and G b F for a fixed fermionic bond configuration with fermion number F = 0, 1. We point out that for an arbitrary superpotential there are in general no restrictions on the bosonic bond configurations. This is for example the case for the superpotential P b which we consider in our series of papers, cf. eq.(53). In contrast, the superpotential P u in eq.(54) yields the local constraint N (x) = 0 mod 2 on the site occupation number, due to the parity symmetry φ → −φ. In the following discussion, we always present the generic case first, and then specify the modifications or simplifications due to the constraint. In analogy to the fermionic bond update, the 'put/remove' and the 'shift' updates are the main steps for updating the bosonic bond configurations. The 'put/remove' step introduces or removes one or two sources φ, while the 'shift' step shifts the sources by one lattice spacing. If there are no restrictions on the bond configuration, we are free to decide for each Monte Carlo step whether to proceed by a 'remove' update or a 'shift' update. With probability p rm , we propose to remove the sources from the lattice, while the proposition to continue the worm update with a 'shift' step is chosen with probability 1 − p rm .
The step from G b F to Z F (and vica versa) is induced by removing (or introducing) a bosonic source φ at sites x 1 and x 2 , with x 1 = x 2 not excluded. The step does not alter the bond configuration, but only the site occupation numbers at sites x 1 and x 2 . Thus, only the ratios of the site weights Q F are involved in the acceptance probability for e.g. the 'remove' step, The prefactor 1/(p rm L 2 t ) is motivated as follows. The factor 1/L 2 t balances the probability for the proposition of putting the bosonic sources at the sites x 1 and x 2 when re-entering the configuration space G b F , while the factor 1/p rm balances the proposition probability for the choice of proceeding by the shift update instead of the remove update, as discussed above. The acceptance ratios for re-entering the configuration space G b F from Z F are given by Two remarks are in order. Firstly, if there are no constraints on the bond configuration, one can in principle introduce just a single source φ which subsequently is shifted around.
In effect, the algorithm then samples the one-point function which in this situation is indeed nonvanishing. Secondly, we note that if the constraint N = 0 mod 2 is in place, the two sources can only be placed or removed when x 1 = x 2 . As a consequence, only the first of the two acceptance ratios in eq.(20) and eq.(21) are relevant, while the second ones are zero by definition.
Next, we discuss the bosonic 'shift' update. With this step we now change the bosonic bond configuration. Shifting the source from site x to a next neighbouring site y is always associated with an increase or a decrease of the bosonic bond occupation number between the sites x and y by one. Whether or not the occupation number is increased or decreased is decided with probability 1/2. Similarly, the source can move forward or backward, and we propose both directions with equal probability 1/2. In addition, when there are several types of bosonic bonds b i with i ∈ {j → k|j, k ∈ N}, we need to decide in each step which bond is updated. We do so by choosing the proposition probabilities p j→k with j,k p j→k = 1. However, because the proposals are completely symmetric, these probabilities do not affect the acceptance ratios. In the following, we will use the shorthand notation for the occupation number of the bosonic bonds b j→k between the sites x and y. The shifts x → y and n j→k xy → n j→k xy + 1 are balanced with shifts y → x and n j→k xy → n j→k xy − 1, which gives the acceptance ratios P sh (x → y, n j→k xy → n j→k xy + 1) P sh (x → y, n j→k xy → n j→k xy − 1) Of course, these generic ratios simplify considerably for the specific bonds b i , i ∈ {1 → 1, 1 → 2, 1 → 3} relevant for the superpotentials considered in our series of papers. For example, the acceptance ratios for updating the bond b 1→1 read Because the bond is symmetric, there is no need to distinguish whether y = x+1 or y = x−1.
To complete the discussion of the bosonic bond update, we point out that the algorithm again provides improved estimators for the bosonic two-point function g b F (x) and the partition functions Z F . As in the fermionic case, the algorithm samples directly the configuration space G b F with the correct weighting when the sources are present. Therefore, every configuration contributes unity to the stochastic Monte Carlo estimator for g b F (x), and we have where x 1 and x 2 are the positions of the two sources. Whenever the bosonic update decides to remove the sources, we have a configuration in Z F and hence a contribution of unity to the stochastic estimator for Z F , that is, we have In complete analogy to the fermionic update we note that the factors of L t appearing in the acceptance ratios of the 'put/remove' step can be compensated by adjusting the overall normalisation of the two-point function, e.g. by making use of translational invariance.

Calculation of the site weight ratios
In order to calculate the weight of a bond configuration, it is necessary to know the site weights where V (φ) and M (φ) depend on the superpotential and the discretisation employed, and F = 0, 1 is the fermion number, for arbitrary values of the site occupation number n. The values of n required in practice are usually limited to O(10 3 ). However, it turns out that even for moderate values of n of order O(100) the site weights Q F (n) can quickly grow larger than 10 100 or more. As a consequence, the calculation of the site weights quickly becomes numerically unstable for growing n. In fact, even for simple potentials when the weights can be calculated analytically in terms of confluent hypergeometric functions, the numerical evaluation of these functions is difficult for large n, and even specialised libraries such as the ones available in Wolfram's Mathematica [22] appear not to be accurate enough.
Fortunately, for the Monte Carlo simulations we only need ratios of the site weights, such as Q F (n + 1)/Q F (n), Q F (n + 2)/Q F (n) and Q 1 (n)/Q 0 (n), and these ratios usually do not become larger than O(10) even for large n. In addition, also the transfer matrix elements can be rewritten in terms of these ratios as discussed in the appendix of our second paper of the series [5]. Therefore, we now present a numerically stable computational strategy to calculate the site weight ratios reliably for arbitrary values of the site occupation numbers.
We start by defining an arbitrary polynomial superpotential and the corresponding bosonic self-interaction potential V (φ) as well as the monomer weight M (φ), Explicitly, the weights in each sector are then given by and For convenience we also define the ratios of the site weights Q F (n), which are used for the acceptance ratios in the Monte Carlo simulations. In principle, only the ratios R 1 (n) need to be calculated since all other ratios can be derived from those. For example, R 1 (n) can be expressed in terms of R 1 (n) as but since in some cases Q 1 (n odd) = 0 the introduction of R 1 (2n) is nevertheless necessary. R m (n) can be expressed via the ratios R 1 (n) and R 1 (n) and appropriate products thereof, R m (n) = m 0 +R 1 (n) (m 1 + R 1 (n + 2) (m 3 + . . .))+R 1 (n) (m 2 + R 1 (n + 2) (m 4 + . . .)) , (38) and the ratios R 0 (n) and R 0 (n) via R m (n), R 1 (n) and R 1 (n) by First, we now discuss how to gain numerical stability for the special case of an even superpotential P (φ). In a second step we will then adapt the idea to treat the somewhat more subtle case of an arbitrary superpotential.

Even superpotential
Unbroken supersymmetric quantum mechanics requires a superpotential P (φ) with deg(P (φ)) = 0 mod 2. In particular, in our series of papers we investigate the superpotential which is symmetric w.r.t. the parity transformation φ → −φ. As a consequence of the symmetry, Q F (n odd) = 0 for both F = 0, 1 and the ratios R F (n) need not be consideredinstead, it is sufficient to determine R 1 (2n) with n ∈ N 0 only.
For the the potential V (φ) we then have the form consistent with both the standard discretisation and the Q-exact one. To keep the integrals numerically under control, for fixed n we apply a variable transformation φ → φ/ φ to obtain rescaled weights Q 1 (2n) as Since we have Q 1 (2n) ≥ 0, we can choose the rescaling factor to be φ = Q 1 (2n) 1/(2n+1) and the rescaled weight becomes Q 1 (2n) = 1. Calculating the ratio of rescaled weights as where both integrals Q 1 (2n+2) and Q 1 (2n) are rescaled with the same factor φ = Q 1 (2n) 1/(2n+1) , we find that In addition, the rescaled weight Q 1 (2n + 2) is now of O(1) and can be evaluated reliably via numerical integration. So if we start by integrating directly the numerically stable site weights Q 1 (0) and Q 1 (2), we can recursively generate ratios R 1 (2n) with higher and higher n. Note that after each calculation of a ratio R 1 (2n), one needs to update the rescaling factor φ → φ . This can be achieved most easily via Our procedure guarantees that all involved quantities are of O(1). Once all ratios R 1 (2n) are known, one can calculate the ratios R m (2n), noting that for the specific superpotential we consider, eq.(38) simplifies to The calculation of the ratios R 0 (2n) as given in eq.(39) is then straightforward.

Arbitrary Superpotential
In the context of broken supersymmetric quantum mechanics, one encounters superpotentials with deg(P (φ)) = 1 mod 2. Therefore, we now adapt the procedure from above to superpotentials of this form. For simplicity, we restrict ourselves to the odd superpotential we consider as the example in our series of papers, If at least one of the coefficients p 1 and p 2 is nonzero, which is always the case for the superpotentials we use, V (φ) reads and at least one of the coefficients k 1 and k 3 is nonzero either. This has a two important consequences. Firstly, the moments defined in eq.(32) are nonzero for n odd, from which it follows that the ratios R F (n) defined in eq.(36) have to be calculated as well. Secondly, the weights Q 1 (n) are no longer necessarily positive. It turns out, however, that for all practical purposes it does not affect the simulations. We will discuss this further in Section 4.
For the evaluation of the integrals, we apply the same variable transformation φ → φ/ φ as before, such that we have rescaled weights Q 1 (n) given by We now choose φ = |Q 1 (n)| 1/(n+1) · sgn(Q 1 (n)). Then, the integral becomes Q 1 (n) = 1 again as before. Furthermore, defining the rescaled ratios R 1 (n) to be where both integrals Q 1 (n+1) and Q 1 (n) are rescaled with the same factor φ = |Q 1 (n)| 1/(n+1) · sgn(Q 1 (n)), we find R 1 (n) = φ R 1 (n). We proceed analogously to the case of the even superpotential by recursive iteration, with the only exception that we generate the ratios R 1 (n) instead of the ratios R 1 (n). The update for the rescaling factor φ → φ is done via Once all the ratios R 1 (n) are known, one can calculate the ratios R 1 (n) via eq.(37), the ratios R m (n) via eq.(38), and the ratios R 0 (n) and R 0 (n) via eq.(39) and (40), respectively.

Results of the Monte Carlo simulations
The results in this section are merely thought of as a proof of the feasibility of the algorithm and as a test of its efficiency. Comparing the Monte Carlo results with the exact solution of the system at finite lattice spacing provided in our second paper [5] of course also serves as a validation for the algorithm. We refer to that paper for a thorough discussion and physical interpretation of the results.
For the following Monte Carlo simulations, we consider the same superpotentials and discretisations as in the previous two papers. In particular, we simulate the system using the action with counterterm for both unbroken and broken supersymmetry as well as the Q-exact action for unbroken supersymmetry. Details for the various actions can be found in the first paper of our series. Here we only give the details of the superpotentials for broken and unbroken supersymmetry, respectively, and we recall that the continuum limit is taken by fixing the dimensionful parameters µ, g, λ and L while taking the lattice spacing a → 0. In practice, the dimensionless ratios f u = g/µ 2 and f b = λ/µ 3/2 fix the couplings and µL the extent of the system in units of µ, while aµ and a/L are subsequently sent to zero. In analogy to the number of sweeps for a standard Monte Carlo simulation, we count the number of times the algorithm is in either one of the two configuration spaces Z F , F = 0, 1. The statistics for a simulation is therefore given by First, we consider the standard discretisation with the superpotential P u such that supersymmetry is unbroken. As a first observable, we show the results for the bosonic and fermionic correlation functions for µL = 10, L/a = 60 and f u = 1 for Z a = 10 7 in figure 5. This is essentially the same plot as figure 10(b) in our second paper [5], but now with the additional data from the Monte Carlo simulation and plotted on a logarithmic scale. Note that we use the notation x = t in accordance with [5]. The simulation indeed reproduces the exact result within very small statistical errors which demonstrate the efficiency of the algorithm. The exponential error reduction is due to the use of the improved estimators for the two-point function which are available in the context of the worm algorithms. The improvement is particularly impressive for the fermionic correlator where the error reduction allows to follow the correlator over more than seven orders of magnitude without loss of statistical significance. In fact the relative error for the lowest value of the fermionic correlator is still only 4%.
As a second example, we show the mass gaps for different µL at a coupling f u = 1 with statistics of Z a = 10 6 in figure 6. The µL considered are in the region where thermal effects are negligible and essentially only Z 0 contributes to the total partition function, such that Z a Z 0 . We extract the masses from the asymptotic behaviour of the correlation function at large t, i.e., we extract the lowest energy gap. Because of the extremely good signal-to-noise ratio the asymptotic behaviour can be truly reached and, in doing so, systematic errors from contributions of excited states are essentially excluded. Of course, we know from our exact results that the overlap of the simple operators we use to construct the two-point function is close to maximal. This is clearly visible in figure 5 where we observe an almost purely exponential decay for all t/L. Because the energy gaps are independent of µL, they are expected to fall on top of each other for all values of µL at fixed lattice spacing aµ. This is indeed the case within our numerical accuracy, and the extracted masses, when expressed in units of µ, indeed extrapolate to the correct zero-temperature continuum limit. The inset of figure 6 shows a detailed comparison of the simulation results with our exact solution from [5] represented by the dashed line and we observe a beautiful agreement even very close to the continuum.
Next, we consider the action with counterterm and the superpotential P b for which the supersymmetry is broken. In this case we encounter an issue concerning the potential nonpositivity of the weights which we already mentioned in Section 3.2. This potentially dangerous sign problem is not of fermionic origin, but is instead related to the bond formulation of the bosonic degrees of freedom. As a matter of fact it occurs already in the purely bosonic system, independent of the dimensionality of the system. However, negative weights only occur in a region of parameter space which becomes irrelevant towards the continuum limit. In that sense, the sign problem is a lattice artefact and can be avoided straightforwardly. Nevertheless, in order to eliminate any systematic error we deal with this bosonic sign problem by incorporating the sign of the configuration into the observables, even though it has no practical consequences.
As a first observable in the broken case, we show the bosonic and fermionic two-point functions, φ t φ 0 and ψ t ψ 0 , for periodic and antiperiodic b.c. for µL = 10 at fixed coupling f b = 1 in figure 7 for a statistics of Z a = 10 8 . The exact results from [5] are shown as dashed lines. The simulation yields results which agree with the exact results within the very small statistical errors on the level of 1 . Note that the correlators for periodic and antiperiodic b.c. are constructed a posteriori from the simulation results in the bosonic and fermionic sectors Z 0 and Z 1 , respectively, and it is crucial to sample the relative weight between the two sectors correctly in order to get the final values right. The relative sampling is solely in the responsibility of the fermion simulation algorithm. Our results in figure 7 show that the open fermion string algorithm indeed transits sufficiently well between the two sectors.
This statement can be made more quantitative by looking at the ratio Z p /Z a which represents the Witten index in our field theoretic setup. From our exact results in [5] we expect a nonzero Witten index at finite lattice spacing which however extrapolates to zero in the continuum limit. So the behaviour of the algorithm towards the continuum limit is particularly interesting, because for vanishing lattice spacing the would-be Goldstino at finite lattice spacing turns into a true, massless Goldstino. In such a situation one usually encounters critical slowing down of the simulation algorithms, such that the errors on the results grow large and the results become unreliable. The massless Goldstino is directly related to the tunnelling between the bosonic and the fermionic sector, and the reproduction of a Witten index W = 0 in the continuum with small errors is hence a true demonstration of the efficiency of the open fermion string algorithm to transit between the bosonic and fermionic sector. In addition, we know from [5] that the lattice artefacts are exponentially enhanced towards zero temperature and it is interesting to see how the simulation algorithm handles this situation at coarse lattice spacing.
In figure 8 we show the ratio Z p /Z a as a function of the lattice spacing aµ for different values of µL at fixed coupling f b = 1. For this quantity, too, the simulation yields results which agree with the exact results within the small statistical errors. Moreover, the efficiency of the algorithm does not appear to deteriorate towards the continuum limit or for small values of µL where the Witten index is very close to zero. This can for example be seen from the fact that the errors obtained with fixed statistics essentially remain constant towards the continuum limit and are also independent of the system size. This nicely demonstrates the efficiency of the algorithm also for a system with broken supersymmetry.
The last system we investigate with the worm algorithm is unbroken supersymmetry formulated with the Q-exact action 1 . We first consider the ratio of partition functions Z p /Z a which in the limit of µL → ∞ yields the Witten index. From a simulational point of view, the ratio essentially calculates the fraction of configurations in sector Z 0 versus the ones in Z 1 . For unbroken supersymmetry the system is almost exclusively in the bosonic sector, and hence the ratio is very close to one except when the size of the system becomes very small, i.e., in the high temperature limit. Moreover, from our exact results in [5] we know that the lattice artefacts in this quantity are very small and the continuum limit is not very interesting. For these reasons, we consider in figure 9 the dependence of the ratio Z p /Z a on µL for different values of the lattice spacing a/L with a statistics of Z a = 10 8 .
Also for this quantity, we find that the results agree with the exact result within the very small statistical errors. Again, the open fermion string algorithm proves to be very efficient even close to µL 0 where the tunnelling from the bosonic to the fermionic sector and vice versa becomes important and dominates the behaviour of the system. Thus, even in this somewhat extreme situation of very high temperature, the algorithm does not show any signs of critical slowing down despite the fact that there is a quasi-zero mode in the system.
Note that the algorithm is capable of handling negative bare masses independent of the discretisation used and fig.9 is simply also an illustration of this fact.
The last quantity we calculate are the lowest bosonic and fermionic mass gaps for different µL at fixed coupling of f u = 1 from a statistics of Z a = 10 6 . The mass gaps are extracted from the two-point correlation functions exactly in the same way as before for the standard action, and in figure 10 we show the results of this analysis. As expected, the masses for the boson and the fermion are indeed indistinguishable within statistical errors. The degeneracy of the masses at finite lattice spacing due to the Q-exactness of the action emerges also for the results from Monte Carlo simulations. Note that the chosen values for µL lie well within the region where thermal effects are negligible and the masses extrapolate nicely to the correct zero-temperature continuum limit. The inset in figure 10 shows a detailed comparison with our exact results from [5] and we again observe beautiful agreement.

Conclusions
In this paper we present an algorithm for simulating N = 2 supersymmetric quantum mechanics on the lattice. The algorithm is based on the reformulation of the system in terms of bosonic and fermionic bonds, and in essence represents an efficient Monte Carlo scheme for updating fermionic and bosonic bond configurations. The updating of the fermionic degrees of freedom is of specific interest, because this is in general the most challenging part of a simulation. This is particularly true for systems with broken supersymmetry, where standard simulation algorithms suffer from critical slowing down due to the massless Goldstino mode. In addition, these systems inevitably also suffer from a sign problem related to the Goldstino and the vanishing Witten index.
In contrast, the fermion simulation algorithm proposed in [2] eliminates critical slowing down by directly sampling the fermionic two-point correlation function. It is based on introducing a fluctuating open fermion string which efficiently updates the bond configurations on all length scales up to the correlation length associated with the fermionic correlation function. As a consequence, the fermion string induces frequent tunnellings between the bosonic and fermionic vacuum when that correlation length becomes large. Since the two vacua contribute to the partition function with opposite signs, the frequent tunnelling guarantees sufficiently small statistical fluctuations for the average sign, and hence a solution to the fermion sign problem. In fact, the more severe the sign problem gets towards the continuum limit, the more efficiently the algorithm tunnels between the bosonic and fermionic sectors. This is of course due to the growing correlation length associated with the vanishing Goldstino mass. The bosonic degrees of freedom can be expressed in terms of bonds as well. Therefore, we also give the details of an updating algorithm for the bosonic bond configurations. Since we consider Q-exact discretisations in addition to the standard one, the algorithm involves updating generic types of bonds.
The simulation algorithm requires the calculation of the site weights Q F (N ). Their numerical evaluation, however, turns out to be numerically unstable for large site occupation numbers N . Hence, in Section 3, we devise a computational strategy which allows to reliably evaluate the ratios of weights for arbitrarily large occupation numbers. Since this is a generic problem occurring in the bond formulation of field theories with real scalar fields, such a computational scheme is useful also in other situations.
Finally, we present a selection of results obtained using the open fermion string algorithm. We concentrate on two specific realisations of supersymmetric quantum mechanics, one with broken and one with unbroken supersymmetry. In addition, we consider both the standard and the Q-exact discretisation. Since exact results are available at finite lattice spacing from our investigation in [5], we can benchmark our stochastic results and directly validate them. The calculation of the bosonic and fermionic correlation functions shows that they can be determined very accurately over several orders of magnitude. This allows for a very precise computation of the boson and fermion masses, the latter in many cases with a smaller error than the former. In general, a precision of 1 can be reached with a very modest computational effort. In systems with broken supersymmetry it is crucial that the simulation algorithm efficiently samples the relative weights between the bosonic and fermionic sectors. Our results for the partition function ratio Z p /Z a , i.e., the Witten index, show that this is indeed the case. For fixed statistics, the errors do not grow towards the continuum limit. In that limit the index gets very close to zero and the sign problem would therefore be most severe. Similarly, the error is essentially independent of the system size, which shows that the sign problem is truly solved.