Restoration of supersymmetry on the lattice: Two-dimensional $\mathcal{N}=(2,2)$ supersymmetric Yang-Mills theory

By numerically investigating the conservation law of the supercurrent, we confirm the restoration of supersymmetry in Sugino's lattice formulation of the two-dimensional $\mathcal{N}=(2,2)$ supersymmetric SU(2) Yang-Mills theory with a scalar mass term. Subtlety in the case without the scalar mass term, that appears to ruin perturbative power counting, is also pointed out.


Introduction
Nonperturbative study of supersymmetric gauge theory is of great general interest but in the context of lattice formulation, no compelling evidence of supersymmetry has been observed so far. Spacetime lattice is generally irreconcilable with supersymmetry and one must fine-tune coefficients of relevant and marginal operators so that supersymmetric Ward-Takahashi (WT) identities are restored in the continuum limit. (An important exception is the four-dimensional N = 1 supersymmetric Yang-Mills theory (SYM) [1]; see Ref. [2] for a lot of effort went into numerical study of this system; see Ref. [3] for a recent attempt.) Recently, for two-and three-dimensional extended supersymmetric gauge theories, lattice formulations that require no (or a little) fine tuning have been proposed [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]. 1

It can
Email addresses: kanamori-i@riken.jp (Issaku Kanamori), hsuzuki@riken.jp (Hiroshi Suzuki). 1 For relations among these formulations, see Refs. [19,20,21,22,23,24,25]. be argued that exact fermionic symmetries in these lattice formulations, combined with other lattice symmetries, (almost) prohibit relevant and marginal supersymmetry breaking operators to appear [5,7]. Supersymmetry is then restored in the continuum limit without (or with a little) fine tuning. 2 The aim of the present study is to test this scenario of supersymmetry restoration, in Sugino's lattice formulation of the two-dimensional N = (2, 2) SU(2) SYM [7,8]. By a Monte Carlo simulation, we study a supersymmetric WT identity in the form of the supercurrent conservation law. For the reason elucidated in Sec. 4, we introduce a supersymmetry breaking scalar mass term to the lattice formulation. The expected WT identity hence takes the form of a partially conserved supercurrent ("PCSC") relation in which the breaking term is proportional to the square of the scalar mass. We numerically confirm the restoration of this PCSC relation in the continuum limit. Our result strongly indicates that the proposed scenario for the supersymmetry restoration is in fact valid and the target continuum theory (i.e., the two-dimensional N = (2, 2) SU(2) SYM with a soft supersymmetry breaking mass term) is realized in the continuum limit of the present lattice model.

Preliminaries in the continuum theory
The euclidean continuum action of the two-dimensional N = (2, 2) SU(N c ) SYM, that is obtained by dimensionally reducing the N = 1 SYM from four to two dimensions, is given by 3 2 See also Refs. [26,27] for alternative approaches that do not use exact fermionic symmetries. 3 In this paper, we adopt the representation and  , and set Ψ T = (ψ 0 , ψ 1 , χ, η/2) that corresponds to the conventional basis in the topological field theory [28,29]. See also Ref. [30].
where Roman indices M and N run over 0, 1, 2 and 3, while Greek indices µ and ν below run over only 0 and 1. Here, in all expressions, it is understood that ∂ 2 → 0 and ∂ 3 → 0 (dimensional reduction). We also define complex scalar fields, φ ≡ A 2 + iA 3 and φ ≡ A 2 − iA 3 . The auxiliary field H is related to the auxiliary field H in Refs. [7,8] by H ≡ H − iF 01 . We also introduce a soft supersymmetry breaking term This term is "soft" in the sense that it does not introduce new ultraviolet divergences compared with the supersymmetric theory S.
where the parameter ǫ has four spinor components, ǫ T ≡ −(ε (0) , ε (1) , ε, ε). We define components of the super transformation by Also, the total action S + S mass possesses following global symmetries: The and a global Z 2 symmetry, The last one is a remnant of the reflection symmetry with respect to the M = 2 direction of the four-dimensional theory before dimensional reduction.
A Noether current associated with the supersymmetry of S (the supercurrent) is given by These four spinor components correspond to fermionic transformations in Eq. (5), Q (0) , Q (1) , Q and Q, respectively. If the auxiliary field H, instead of H, is regarded as an independent dynamical variable as is the case below, F 01 = −F 10 in this expression must be understood as −iH.
We will consider a correlation function of the supercurrent and another fermionic operator. As a lowest-dimensional gauge invariant fermionic operator, we take This is gauge invariant because the scalar fields A 2 and A 3 , as well as the spinor Ψ, transform as the adjoint representation under two-dimensional gauge transformations. We in particular consider a "diagonal part" in the product of supercurrent (9) and operator (10). That is, we consider following four correlation functions It is useful to note that these four correlation functions are interchanged each other under transformations (7) with α = π/2 and (8). Since Eqs. (7) and (8) are symmetries of the present system, four correlation functions in Eq. (11) coincide in the continuum theory. On the other hand, the present lattice formulation is not invariant under either Eq. (7) or Eq. (8) (while Eq. (6) is manifestly realized) and a "degeneracy" of these four correlation functions can be used as a useful measure for how we are close to the continuum.
With the presence of the supersymmetry breaking term S mass (4), the supercurrent is not conserved and, defining we obtain the supersymmetric WT identity (the PCSC relation) in spinor components if the regularization preserves supersymmetry. We emphasize that these local relations hold irrespective of boundary conditions. One can derive these relations in the functional integral by employing a local change of variables that does not "touch" the boundary. See Appendix A. Therefore, these hold in particular with the antiperiodic temporal boundary condition for fermionic variables, although this boundary condition explicitly breaks supersymmetry (for example, the energy spectrum would not be supersymmetric with this boundary condition).
We now argue that, if the argument for the supersymmetry restoration in Ref. [7] is valid, relations (13)-(16) must be reproduced for x = 0 in the continuum limit of the lattice model.
According to the argument of Ref. [7], the lattice action in Ref. [7] provides a regularization, that becomes supersymmetric in the continuum limit, for all 1PI correlation function of elementary fields. In particular, possible ultraviolet divergent functions, tadpoles and the scalar self-energy at the one-loop level, take the form (in the continuum limit) expected in the continuum theory (i.e., tadpoles vanish and no self-energy correction at zero external momentum). In this sense, the present lattice regularization is supersymmetric. This almost implies that supersymmetric WT identities such as Eqs. (13)-(16) hold in the continuum limit.
Eqs. (13)- (16), however, contain composite operators and a definition of composite operators in a lattice formulation is to a large extent arbitrary. We thus must be sure that lattice artifacts, when combined with ultraviolet divergences arising from these operators, do not modify the WT identities. Fortunately, in the present two-dimensional super-renormalizable system, operators in Eqs. (9), (10) and (12) themselves are ultraviolet finite (i.e., no operator renormalization/mixing is required; this significantly simplifies consideration of the supersymmetry WT identity compared with four dimensions). 4 Only possible ultraviolet divergence in Eq. (16), for example, arises from a one-loop diagram obtained by connecting elementary fields in J µ and X ν . However, this ultraviolet divergence can readily be avoided by setting two points x and 0 to be separated. 5 In fact, by analysing this one-loop contribution, one has for Eq. (16) and similar corrections for Eqs. (13)-(15) (the constant c Q here can differ for each of relations, as c Q (0) , c Q (1) and c Q , depending on the regularization), assuming that the regularization for 1PI functions of elementary fields preserves supersymmetry. The constant c Q in Eq. (17) depends on the definition (or regularization) of composite operators J µ and X ν , but this de-pendence on regularization disappears for x = 0. 6 In summary, if the argument for the supersymmetry restoration in our lattice model is valid, relations (13)-(16) with x = 0 must be reproduced in the continuum limit. 7 One might wonder why we worry about the supersymmetry restoration in the present model; after all, supersymmetry is explicitly broken by the scalar mass term (and by the boundary condition that we will adopt below). It is very important, however, to distinguish three different possible sources for supersymmetry breaking in the present lattice model; the scalar mass term, the boundary condition and the lattice regularization itself. Our point is that the observation of WT identities (13)-(16) enables us to isolate the last source of supersymmetry breaking. What we are talking about is whether this supersymmetry breaking owing to the lattice regularization disappears in the continuum limit or not.

Results of Monte Carlo study
For details of the lattice formulation in Refs. [7,8], we refer the reader to the original references and Refs. [33,35]. The point is that the fermionic transformation Q in Eq. (5) is nilpotent (up to gauge transformations) and action (3) can be written as a Q-exact form, as topological field theory [28,29]. These nilpotent Q and Q-exact action, and thus the invariance of the action under Q, can be realized even in lattice gauge theory [7,8]. This sort of exact fermionic symmetry is realized also in lattice formulations of the present system in Refs. [4,5,10,12,13]. Other fermionic symmetries, Q (0) , Q (1) and Q, from which Eqs. (13)-(15) follow, are not preserved in the lattice formulation and the question is whether the invariance under these is restored in the continuum limit.
We consider a finite two-dimensional rectangular lattice where a is the lattice spacing and the physical size L is fixed to be √ 2/g. 8 6 Incidentally, the last term of Eq. (17) is not a genuine anomaly of supersymmetry. In fact, there exists a possible definition of J µ that leads to c Q = 1; see footnote 11. 7 Our argumentation here is not quite rigorous. For more satisfactory treatment, one first derives a WT identity associated with an appropriate (would-be) supersymmetry transformation in lattice theory and then shows all lattice artifacts to Eqs. (13)-(16) with x = 0 vanish in the continuum limit. 8 Note that in two dimensions, the gauge coupling g has the mass dimension 1. Fig. 7 below, all the results were obtained with the antiperiodic temporal (= x 0 ) boundary condition for fermionic variables (the x 1 -direction is always periodic). We found that this boundary condition leads to faster approach to the continuum and less noisy signals, compared with the periodic boundary condition. Recall that WT identities (13)- (16) in the continuum theory must hold irrespective of the boundary condition. 9

Except the result in
Corresponding to mass term (4), we add a scalar mass term to the lattice action of Ref. [8]. 10 For lattice transcription of operators in Eqs. (9), (10) and (12), we adopt a simple prescription that a field variable at a point x is replaced by the corresponding lattice field variable at x. Covariant derivatives for the scalar field are replaced by the forward covariant differences, are link variables. As noted in the preceding section, a precise way of lattice transcription of these operators should be irrelevant to our present analysis.
That is, all ambiguities associated with the operator definition are integrated in the constants c Q in Eq. (17) (and c Q (0) , c Q (1) and c Q ) in the continuum limit and they do not appear for x = 0. 11 One of us (I. K.) developed a Rational Hybrid Monte Carlo (RHMC) [36] simulation code (on the basis of the multi-shift [37] CG solver [38]) for the present lattice system with the gauge group SU(2), by utilizing available libraries and programs [39,40,41]. Some details of this simulation code have already been reported in Ref. [42]. (See also Fig. 3 of Ref. [43] that illustrates how the effect of dynamical fermions is properly included with this code.) We summarize on parameters in the present simulation: The parameter ǫ for the admissibility [8] was fixed to be 2.6. 12 The degree of the rational approximation in RHMC was typically 20-30. The multi time step acceleration [44] was used in the molecular dynamics (see Appendix B). The time interval of one trajectory was fixed to be 0.5 and time steps of leapfrog were chosen such that the acceptance was greater than 0.8. We stored configurations in every 10 trajectories. The autocorrelation was then estimated by jackknife analysis with binning. We discarded first 10000 configurations for thermalization and used subsequent configurations at every 50 configurations. In this way, we prepared uncorrelated configurations listed in Table 1.
In the present lattice formulation, the pfaffian of the Dirac operator resulting from the integration over fermionic variables, that is real positive in the continuum theory, is generally complex owing to lattice artifacts. For parameters in the present simulation with the antiperiodic boundary condition, we found that the absolute value of the complex phase of the Dirac determinant is 0.8 radians. We took into account the complex phase of pfaffian by reweighting the half of the phase of the determinant. See Refs. [33,35]. We observed that, however, this reweighting had practically no notable effect. Also, all quantities presented below, that are real in the continuum theory, are generally complex with finite lattice spacings. The imaginary part is however generally small (and consistent with zero) and we will present only the real part for simplicity.
For the correlation functions used in the following analyses, we took an average over the whole lattice points to increase the effective number of configurations. For example, J is computed by taking the average of J (0) 0 (y)X (0) 0 (z) over all y ∈ Λ and z ∈ Λ with y 0 − z 0 = x 0 mod β and y 1 − z 1 = x 1 mod L are kept fixed. 13 12 Although a precise choice of the parameter 0 < ǫ < 2 √ 2 should be irrelevant for results with small lattice spacings, we admit that we have not carried out systematic study on this point. As a simple consistency check, we performed a small experiment that measures the expectation value of the action density of the pure Yang-Mills part (Eq. (4.27) of Ref. [35]) for different values of ǫ, ǫ = 1.0, 2.0, and 2.6. The scalar mass squared is µ 2 /g 2 = 0.25, lattice size is 6 × 6, and the lattice spacing is ag = 0.2357. The number of uncorrelated configurations is 100 for each ǫ. We observed that the admissibility bound 1 − U (x, 0, 1) < ǫ (see Ref. [35]) was never exceeded in the molecular dynamics and the expectation values for each ǫ coincided within statistical errors. 13 With the understanding that the sign of the correlation function is flipped when y − z = x − β0, the correlation function is translationally invariant even with the As a typical example of correlation functions, we plot in Fig. 1 correlation functions (11) obtained by employing set IV (a) of Table 1 as a function of x 0 (x 1 = L/2). Four correlation functions are quite well-degenerated, except at points near x = 0 and its periodic images. (As noted, at x = 0 and its periodic images, the correlation functions suffer from lattice artifacts associated with the operator definition.) We thus have an indication that the lattice spacing ag = 0.2357 is already rather close to the continuum, at least for the correlation functions with µ 2 /g 2 = 1.0.
It is instructive to plot the left-hand side of WT identities (13)- (16). If supersymmetry is restored, the left-hand side must vanish except at x = 0 and its periodic images; recall Eq. (17). In Fig. 2, we showed the left-hand side of Eq. (13) obtained by set IV (c) of Table 1. Here, we used the symmetric difference, , as a lattice transcription of the total divergence, because we found that its use considerably reduces the statistical error. We see that the left-hand side is almost zero everywhere (the absolute value is less than 0.05 in the central portion) and it is sharply peaked at the origin x = 0 and its periodic images. This result strongly indicates that antiperiodic boundary condition.  Table 1: our reasoning in Sec. 2 is correct and, at the same time, that supersymmetry is certainly restored in the continuum limit. Recall that Eq. (13) is associated with the fermionic symmetry Q (0) that is not preserved in the present lattice model with finite lattice spacings.
For systematic quantification of the supersymmetry restoration, we consider the following four ratios: According to PCSC relation (13)- (16), these ratios must become µ 2 /g 2 in the continuum limit, except at x = 0 and its periodic images. We plotted  Table 1. Fig. 3. The statistical error in the ratio was estimated by jackknife analysis. This plot is for µ 2 /g 2 = 1.0 and in the continuum limit the points should lie on the dotted line. We see this expected tendency, while a deviation near x 0 = 0 ≡ β can be understood as the effect of approximate delta functions at x = 0 and periodic images elucidated above. Also, again, we see that the lattice spacing ag = 0.2357 is already rather close to the continuum for µ 2 /g 2 = 1.0, because the WT identity is fairly restored with this lattice spacing.

ratio (20) for different lattice spacings in
To quantify the restoration of supersymmetry, while eliminating the effect of lattice artifacts existing around x = 0 and its periodic images, we adopted the following procedure. First, we take a cylindrical region C in Λ, C ≡ {x ∈ Λ | β/2−∆ 1 < x 0 < β/2+∆ 2 , 0 ≤ x 1 < L}. We then apply a χ 2 -fit by a constant to ratios (20)-(23) obtained for x ∈ C. We change ∆ 1 and ∆ 2 such that the number of points in the region C becomes maximum while keeping the χ 2 per one degree of freedom (dof) reasonably small. The results of this procedure are summarized in Table 2 and in Fig. 4. In Table 2, columns indicated by "used" show the number of points contained in the region C determined as above. In Fig. 4, the χ 2 -fitted values of ratio (20) obtained by the above procedure are plotted. As shown in the figure, these values were then used for a linear χ 2 -extrapolation to the continuum.  23) obtained by the above procedure is plotted as a function of the parameter µ 2 /g 2 . The result is consistent with the straight line, that is a prediction of the supersymmetric WT identities (the PCSC relation). Some care should be paid for the interpretation of Fig. 5; the plot shows only statistical errors. There might be rather large systematic errors associated with the fitting procedure especially for the smallest mass µ 2 /g 2 = 0.04 case (set I). Figure 6 shows ratio (20) for set I and we see that a flat region is narrower compared with Fig. 3, that is for µ 2 /g 2 = 1.0. In fact, as Table 2 shows, the numbers of points we used in the fit are rather fewer for set I than those for other sets. From Fig. 6, however, the systematic errors for set I would be at most 0.05, thus the result is still consistent with the PCSC relation even if this systematic error is taken into account. For other values of the mass parameters µ 2 /g 2 , the behavior of ratios is more or less similar to that in Fig. 3 and the numbers of points we used in the fit appear sufficient; so we do not expect large systematic errors. Therefore, from the agreement with the theoretical expectation in Fig. 5, we infer that supersymmetry, that is broken in the present lattice regularization with finite lattice spacings, is certainly restored in the continuum limit. 14 Table 2 The quality of the χ 2 -fit of ratios (20) The antiperiodic boundary condition explicitly breaks supersymmetry while the periodic one preserves supersymmetry, so it is certainly of interest to carry out simulations (for the energy spectrum, for example) with the latter boundary condition. A typical example of correlation functions with the periodic boundary condition is depicted in Fig. 7. The degeneracy of four correlation functions is not quite realized. This indicates that the spacing ag = 0.2357 is not yet so close to the continuum for µ 2 /g 2 = 1.0 when we use the periodic boundary condition. In addition to this, we found that signals with the periodic boundary condition are generally rather noisy and, when translated to the ratio like Fig. 3, the statistical errors are too large for a reliable fit. For these reasons, we postpone a detailed study of cases with the periodic boundary condition to a future work.
firmed that in fact the results significantly differ from those obtained above and the expected supersymmetric WT identities are not restored. 4 Subtlety in the supersymmetric µ 2 = 0 case In the above analyses, we have introduced the scalar mass term (19) that generally suppresses the amplitude of scalar fields. One is of course interested in the original two-dimensional N = (2, 2) SYM that does not contain such a supersymmetry breaking term. We now explain why we had to introduce the scalar mass term in our simulation. Figure 8 shows correlation functions (11) without the scalar mass term. The lattice spacing is ag = 0.2357. One sees that the degeneracy among four correlation functions is enormously violated and thus it appears that we are quite far from the continuum theory. Even if we decrease the spacing to ag = 0.1768 (Fig. 9), the crude feature does not change and the degeneracy is not restored. This behavior prompts us to draw a conclusion that the degeneracy is not restored even in the continuum limit.
Moreover, from a closer look at Figs. 8 and 9, it appears that the nondegeneracy cannot be understood as a breaking of a certain symmetry in the continuum theory. In these figures, only J related each other by certain global symmetries. 15 For example, J 0 (x)X 0 (0) and J 0 (x) X 0 (0) are related by transformation (7) with α = π/2. The fact that these two functions are quite degenerated in the figures suggests that this global discrete symmetry is well restored with the present lattice spacings. However, J 0 (0) are not degenerated at all, although these two are also related by this symmetry in the continuum theory. Something quite strange seems happening.
Another observation we made is that if we modify a term in the lattice Dirac propagator that originates from the Yukawa interaction, by an O(a) amount (in the process of measurement), the effect of the modification is tremendous and four correlation functions become almost degenerated; the effect appears to be O(1).
We suspect that the above strange behavior in the absence of the scalar mass term is caused by very large expectation value of scalar fields along flat directions-continuous set of minima of the classical scalar potential. Here, by "very large", we mean the lattice cutoff scale, O(1/a). In fact, as Fig. 3 of Ref. [42] shows, the (gauge invariant) amplitude of scalar fields can be as large as ∼ 40/a, with the antiperiodic boundary condition. 16 Such a very large expectation value could ruin perturbative power counting, in which the operator aφ, for example, is regarded as O(a). That is, if the expectation value of the scalar field φ is as large as O(1/a), the combination aφ would have to be regarded as O(1). This phenomenon affects also the argument for the supersymmetry restoration. For example, the operator Q(a tr{φψ µ }) = a tr{ηψ µ } + a tr{φiD µ φ} is invariant under the gauge, U(1) A and Q transformations-manifest symmetries of the present lattice formulation-and thus could radiatively be induced (in the one-loop level). This is an irrelevant operator in the usual sense but this could be "marginal" 16 On the other hand, we numerically observed that the amplitude does not grow so much ( 1/a) with the periodic boundary condition. This seems strange at first glance because the periodic boundary condition does not break supersymmetry and one may expect that the flat directions are not lifted upon radiative corrections when supersymmetry is preserved. This phenomenon could be understood if we consider an effective potential for scalar zero modes that is obtained by integrating out all other modes. With the antiperiodic boundary condition, there is no fermion zero mode and the degrees of freedom that integrated out are balancing between bosons and fermions. On the other hand, with the periodic boundary condition, fermionic zero modes are left unbalanced and they can induce nontrivial corrections. See Refs. [45,46] for related observations.
if the scalar field φ is regarded as O(1/a). This O(1) operator, moreover, is not invariant under Q (0) , Q (1) and Q transformations, thus invalidates the argument for the supersymmetry restoration. 17 It is thus conceivable that the target theory, the two-dimensional N = (2, 2) SU(2) SYM, is not realized in the continuum limit of the present lattice formulation, unless we supplement the scalar mass term (or something that suppresses the amplitude of scalar fields). This point is further discussed in the next section. Also, in light of this observation, a study on the dynamical supersymmetry breaking in the two-dimensional N = (2, 2) SU(2) SYM in Refs. [35,47] must be reconsidered [48].

Conclusion and discussion
In this paper, for the first time to our knowledge, the restoration of supersymmetry in a lattice formulation of a supersymmetric gauge theory was clearly observed. The PCSC relation in Fig. 5, first and foremost, can be taken as a solid basis for the lattice model in Ref. [8] to be used for evaluation of various quantities in the two-dimensional N = (2, 2) SYM with a supersymmetry breaking scalar mass term. Also, our result for a two-dimensional model demonstrates validity of general reasoning for the supersymmetry restoration on the basis of lattice symmetries and power counting. It is quite conceivable that, therefore, various lattice formulations of supersymmetric gauge theory that are based on similar reasoning work as they are aimed at.
The above statement is for the lattice model in which there is no flat direction in the classical potential. The lattice model for the original N = (2, 2) SYM without the scalar mass term [8] possesses flat directions and, as we have argued in Sec 4, it is quite plausible that the target theory is not obtained by the continuum limit. Thus, we close this paper by summarizing the situation concerning the N = (2, 2) SYM without supersymmetry breaking terms:   (1) One may first introduce the scalar mass term and then take the limit µ 2 → 0. Physically, as discussed in Ref. [5], if the mass µ is sufficiently small compared with the inverse of the system size, 1/L, the effect of the breaking mass would be practically negligible because one cannot observe the wavelength longer than the size of the "universe". This provides a possible way of defining the N = (2, 2) SYM. This route of definition would be mandatory in lattice formulations in Refs. [5,12,24], in which a supersymmetry breaking scalar mass term must be introduced to stabilize the lattice spacing. Our result in Fig 5 suggests that this prescription works in the formulation of Ref. [8], because we have observed the restoration of supersymmetry even for µ 2 1/L 2 = 0.5g 2 .
(2) One may carry out Monte Carlo simulations without introducing the supersymmetry breaking mass term. This is possible for example in the formulations of Refs. [8,26]. In this approach, however, there will be a subtle problem we encountered in Sec. 4 that O(a) lattice artifacts seem to be amplified to O(1) by O(1/a) expectation value of scalar fields along flat directions. 18 Our observation suggests that it is generally quite difficult to realize supersym- 18 In the large N c limit, there is a possibility that the expectation value of scalar fields along flat directions is suppressed and one can evade the problem without breaking supersymmetry. See Ref. [49]. metric theories with flat directions as a continuum limit of a lattice model, without suppressing the amplitude of scalar fields, because lattice formulation generally cannot be free from O(a) discretization errors.
A profound question is, however, whether above prescription (1) really provides a "correct" definition of the target two-dimensional theory or not; the existence of the vacuum moduli is unlikely in two dimensions while the prescription would enforce scalar fields to localize around the origin. Unfortunately, we do not have the convincing answer to this question at present.
We would like to thank Masanori Hanada, Hikaru Kawai, Martin Lüscher, Hideo Matsufuru and Fumihiko Sugino for enlightening discussions. Our numerical results were obtained using the RIKEN Super Combined Cluster (RSCC). I. K. is supported by the Special Postdoctoral Researchers Program at RIKEN. The work of H. S. is supported in part by a Grant-in-Aid for Scientific Research, 18540305.
A Derivation of supersymmetric WT identities (13)- (16) In this appendix, supersymmetric WT identities (13)- (16) are derived with emphasize on their independence of boundary conditions. We start with the expectation value of operator (10): where Z is the partition function and dµ denotes a measure of the functional integral. We then consider a certain infinitesimal variation of integration variables δ ′ , Since the functional integral itself is independent of any relabeling of integration variables, we have the identity for any variation δ ′ , where combinations s µ and f are given by Eqs. (9) and (12), respectively. In this expression, the second line identically vanishes for any boundary condition because ǫ(x) = 0 at the boundary. Similarly, in the first line, we may perform integration by parts neglecting boundary terms again because ǫ(x) = 0 at the boundary. Finally, setting ǫ(x) to be the delta function (times a Grassmannodd constant spinor), we have Eqs. (13)-(16) as particular cases of Eq. (A.2).
The assumption that the measure dµ is invariant under δ ′ corresponds to, in the present formal treatment, the assumption that regularization preserves supersymmetry.

B Multi time step acceleration in our simulation
Let S B be the action consists only of gauge and scalar fields and S PF the action of pseudofermions. S PF is bi-linear in pseudofermions. In the molecular dynamics, if there exists a definite hierarchy between force originates from S B and that from S PF , and if the latter is smaller than the former, one may reduce computational cost for the latter (it is expensive requiring inversion of a Dirac operator) by using different time steps for each of force. This is the multi time step acceleration [44]. In this scheme, the time-evolution operator in one trajectory ∆τ = nǫ is symbolically written as where T (S, ǫ) denotes the time-evolution operator with respect to the action S in the time step ǫ. That is, one uses a k-times larger time step for force from S PF .
In our Monte Carlo simulation, as shown in Fig. B.1, we found that typically force from S PF is several times smaller than that from S B and the scheme is in fact very efficient. In the example in the figure, our choice was k = 3 and n = 6 in Eq. (B.1) (we set ∆τ = nǫ = 0.5). Since the variation of the force is large (typically the same order of magnitude as the average itself), we chose a smaller value of k than naively expected from the figure. As a general tendency, when a lattice spacing becomes smaller, force from S B becomes larger and, correspondingly, we could use larger k. In fact, we used k = 5, n = 6 for ag = 0.1768, and k = 8, n = 6 for ag = 0.1414.