Direct computational approach to lattice supersymmetric quantum mechanics

We propose a numerical method of estimating various physical quantities in lattice (supersymmetric) quantum mechanics. The method consists only of deterministic processes such as computing a product of transfer matrix, and has no statistical uncertainties. We use the numerical quadrature to define the transfer matrix as a finite dimensional matrix, and find that it effectively works by rescaling variable for sufficiently small lattice spacings. For a lattice supersymmetric quantum mechanics, the correlators can be estimated without statistical errors, and the effective masses coincide with the exact solution within very small errors less than 0.001%. The SUSY Ward identity is also precisely studied in compared with the Monte-Carlo method. Our method is not limited to a lattice SUSY quantum mechanics, but is also applicable to any other lattice models of quantum mechanics.


Introduction
Supersymmetric quantum mechanics (SUSY QM) has been extensively studied since it provides a good testing ground for supersymmetric field theories [1][2][3]. Supersymmetry breaking triggered by instantons has been deeply understood and the associated topological index (Witten index) are widely used in other theories. In early studies, a class of exactly solvable potentials, so-called shape invariant potentials (SIP), was discovered and has attracted much attentions in revealing various aspects of SUSY QM [4]. However, to learn SUSY models beyond those cases, numerical approaches including lattice theory could play a significant role.
Recently, a direct computational method has been proposed by Baumgartner and Wenger [36][37][38]. This approach could help us to solve the sign problem because a product of transfer matrix is directly estimated without any stochastic processes in principle. In this sense, it is similar to the tensor network renormalization (TNR) [39] which has developed over the last decade, including applications to field theory [40][41][42][43][44][45][46][47][48]. In TNR, the partition function and correlators are represented as a product of tensors, which are estimated by coarse graining of network. Those two methods share common idea in the sense that a product of local matrix or tensor is directly estimated by matrix operations. Keeping this point in mind, we may say that further developments in SUSY QM using transfer matrix are associated with TNR and applications to higher dimensional field theories.
In this paper, we propose an explicit and simple computational method in lattice (SUSY) QM, which is based on the transfer matrix representation of Euclidean path integrals. Although the transfer matrix is an infinite dimensional matrix, we express it as a finite dimensional ones by approximating the integrals of coordinate. The well-known Gaussian quadrature is used for the approximation. As we will see later, such simple quadrature effectively works with some scale transformation of variables. In this procedure, we can estimate the correlators in tremendous lattice volume with sufficiently small lattice spacings.
We test our method for two cases without supersymmetry breaking, φ 6 theory and a SIP potential (Scarf I). Although the energy spectra can be directly estimated from transfer matrices, we use the standard technique to evaluate the effective masses from correlators, which are used in the Monte-Carlo method, in order to compare our results with Monte-Carlo ones. The estimated masses reproduce the known results in high precision. For Scarf I potential, energy eigenvalues are exactly solvable, and the lowest one for A/α = 10, B = 0 in (2.7) is m exact /α 2 = 10.5. (1.1) As the effective mass, we have m/α 2 = 10.49998 (2).
Thus one can say that our approach effectively works compared with the previous Monte Carlo simulations. This paper is organized as follow. In section 2, SUSY QM is given in the path integral formulation with Euclidean time, and a lattice theory which partially keeps supersymmetry is also given. We propose a direct computational method using the transfer matrix in section 3. The numerical results for φ 6 -theory and Scarf I potential are given in section 4. We conclude this paper with some discussions in section 5.

Supersymmetric quantum mechanics
We start with reviewing N = 2 supersymmetric quantum mechanics which is defined in Euclidean path integral formulation. The lattice theory is then defined according to Refs. [23] in a form which retains one supersymmetry exactly on the lattice.

Continuum theory
The supersymmetric quantum mechanics is described by a coordinate φ(t) ∈ R and onecomponent Grassmann variables ψ(t) and ψ(t). We assume that they obey the periodic boundary condition with a period β for the Euclidean time t ∈ R. The action of SUSY QM is given by We refer to W (φ), which is an arbitrary function of φ, as a superpotential. 1 For any superpotential, (2.1) is invariant under supersymmetry transformation, where ǫ and ǫ are one-component Grassmann numbers.
In the path integral formulation, the partition function is given by 2 5) and the correlation functions are also defined in the usual manner. For periodic boundary condition, Z is just the Witten index Tr(−1) F which measures the difference of the number of bosonic and fermionic vacuum states. In the case of SUSY QM, we can perfectly learn from the index whether supersymmetry is broken or not: Z = 0 if and only if supersymmetry is broken. As well-known, this index is unchanged under any deformation of parameters which retains the asymptotic behavior of W (φ). In the case of |W (±∞)| = ∞, supersymmetry is unbroken when the signs of W (∞) and W (−∞) are different from each other, while supersymmetry is broken when they are the same [1,2]. In this paper, we focus on two cases which do not exhibit supersymmetry breaking: φ 6 theory and Scarf I potential. The φ 6 theory is given by the cubic superpotential, with the mass m and coupling g > 0, and Scarf I potential is given by where A > B ≥ 0 and α > 0. We have Z = 1 for these two cases. In section 4, we will test our computational method proposed in section 3 for (2.6) and (2.7). The former potential is used to compare our results with the previous Monte-Carlo results since it has been studied in many literatures. The later one provides a comparison between our results and exact solution because it is one of shape invariant potentials (SIPs) which are exactly solvable [3]. 1 For W = d dφ h, h is often referred to as a superpotential in several literatures. 2 The fermion measure is normalized such that Z = 1 at large m in the free theory.

Lattice theory
The lattice theory is defined by considering the bosonic and fermionic coordinates as variables φ t and ψ t living on the integer sites, t ∈ Z. We now take the lattice spacing a = 1 without loss of generality, and often express a-dependence of dimensionful quantities explicitly. We again assume that all variables φ t , ψ t , ψ t satisfy the periodic boundary condition of the period β = N a, where N ∈ N is the lattice size.
The lattice action is given by where the backward difference operator, is used to approximate the derivative ∂ t . More generally, we can use ∇ = ∇ + +∇ − 2 − r 2 ∇ + ∇ − with forward and backward difference operators ∇ + , ∇ − , and the Wilson parameter r. Here we have the backward one by setting r = 1 for simplicity.
We now consider lattice supersymmetry transformation which is defined by replacing ∂ t in (2.2)-(2.4) by ∇ ± as δφ = ǫψ + ǫψ, (2.10) For free theory, the lattice action (2.8) is invariant under this transformation. For interacting theories, it is, however, not invariant for any ǫ and ǫ; one-supersymmetry (ǫ = 0, ǫ = 0) remains as an exact symmetry, while the other (ǫ = 0, ǫ = 0) is broken at finite lattice spacing. The full N = 2 supersymmetry is shown to be restored in the continuum limit as seen in Ref [23,25,27,30,34] and also precisely shown in section 4. The partition function is given by (2.5), and after integrating fermionic variables, we have (2.14) The fermion determinant is well-defined since the size of matrix ∇ + W ′ (φ) is finite. In this case, we can explicitly write it down as 3 which is strictly positive for (2.6) and (2.7). The integrations of φ are given by the finite dimensional integrals on the lattice as Thus the Monte-Carlo techniques can be applied to this lattice model. As seen in the following sections, the transfer matrix representation of (2.13) is very useful to study this model. Using (2.15) and (2.16), we can easily show that where where Tr(A) = ∞ −∞ dφA φφ . We now note that numerical approaches are not applicable to this representation in a straightforward way since the traces in (2.20) can not be performed for the infinite dimensional matrices T φφ ′ and S φφ ′ with φ, φ ′ ∈ R.
The two transfer matrices S and T have specific meanings. In SUSY QM, there are two Hamiltonian H B and H F which act on the fermionic and bosonic states, respectively [3].
i and E F i are sorted in ascending order. In the continuum theory, one can show that The transfer matrices are interpreted as S = exp(−aH B ) and T = exp(−aH F ) on the lattice as discussed in Ref [36][37][38]. Then, we can estimate the energy spectra for bosonic and fermionic states from the eigenvalues of S and T .

Partition function
We first consider the partition function which is the simplest example to describe our computational method. In order to express the transfer matrices as finite dimensional matrices, we replace (2.16) by a numerical quadrature which is given as a weighted summation, 4 where g K (φ) is a weight function and S K is a set of K points which are used for the approximation. g K and S K is determined for each quadrature rule. For instance, the Gauss-Hermite quadrature is provided by and taking S K as a set of roots of K-th Hermite polynomial H K . We expect that in (3.1) the true value of LHS can be obtained within small systematic errors by evaluating RHS at sufficiently large K for well-behaved function F (φ) for which the quadrature effectively works.
The partition function (2.13) is now approximated by replacing each measure of φ by (3.1), and (2.17) turns out to be Since T and S are matrices of size K, (3.3) can be written as where "tr" is the trace of K by K matrices. Note that T and S are not uniquely determined because (3.6) is invariant under similarity transformations of them. We can estimate Z by computing the RHS of (3.6) at large K once an effective quadrature is found. One might think that conventional methods such as trapezoidal rule, Simpson's rule, or a kind of Gaussian quadrature are less effective for multiple integrations since the number of required sample points increases exponentially with multiplicity. Besides, (3.3) is likely ineffective since it is just a superposition of one-dimensional quadrature (3.1). However, numerical results shown in section 4 show that a simple Gaussian quadrature effectively works by using it with a scale transformation of variables, explained in section 3.3.

Correlation functions
The correlation functions are estimated in a similar manner as with the partition function.
As simple examples, we present two-point functions of φ t and ψ t . It is very easy to extend our method to general correlation functions which are given as a product of local operators O i (t). 5 The expectation value of an operator O is defined by For later use, let us define the fermionic part of (3.8) as In our approach, Z and O are individually estimated in contrast to the Monte-Carlo method which provides the expectation value O itself. We can easily show that φ t φ s is written as a product of matrices by expressing the operator insertion as impurity matrix D φφ ′ = φ · δ φφ ′ : where k = (s−t) mod N . 6 The translational invariance and φ t φ s = φ s φ t are manifestly follows from this representation. (3.11) is extended to a case of two local operators made are regarded as a matrix with the column φ and the row φ ′ . Similarly, ψ t ψ s is given in terms of a product of the transfer matrices. Since ψ t is connected to ψ s by the hopping term of (3.10), it is easy to show that, for k = 1, · · · , N , and integrating fermionic variables, (3.14) Instead of (2.15), we have (3.14) as a contribution of the fermion-part. Finally we find that Once the matrix representations like (3.11) and (3.15) are obtained, one can evaluate the expectation value directly by computing the matrix products and traces and dividing them by Z.

Some improvements
The efficiency of our method depends on whether the quadrature is compatible with the lattice action, e.g. W and the difference operator ∇. To improve the efficiency, let us consider a rescaling of φ in (2.13) as Approximating the integrals of Φ by a quadrature, the same formula (3.3) is derived for the modified transfer matrices, In the free theory, Gauss-Hermite quadrature could be effective for large masses ma ≫ 1 by taking s = ma/ √ 2 so that the damping factor is just a Gaussian with rate one, e −Φ 2 . Near the continuum limit ma ≪ 1, one might think that such rescaling is less effective because T and S have long tails in the field space, which corresponds to the zero kinetic term, Φ = Φ ′ . For interacting cases, the situation becomes more complicated since we have to find an effective quadrature for each interaction term.
The Gauss-Hermite quadrature with rescaling φ, however, effectively works by tuning s even for interacting cases, as seen in the numerical results in section 4. We can choose s in such a way that some exact relations are obtained with as high precision as possible. In section 4, s is determined such that |Z − 1| is minimized for each lattice spacing since the Witten index is exactly one thanks to exact supersymmetry. Then we find that the energy eigenvalues and correlators are obtained in high precision even for small lattice spacings.
The computational cost of our method basically behaves as O(K 3 logN ) for the partition function, while it is multiplied by N n−1 for n-point function. The lattice size dependence follows from computing a matrix product T 2n as (T n ) 2 for n = 1, 2, . . ., recursively. The diagonalization of transfer matrices often significantly reduces the cost. We show that (3.6) is expressed as where λ i and ρ i are eigenvalues of S and T , respectively. Solving all eigenvalues of square matrices with size K requires O(K 3 ) as the cost. We actually need a few eigenvalues for large N if the lowest eigenvalues are isolated. Furthermore, for Gaussian quadrature, T and S given as (3.4) and (3.5) are actually sparse matrices since g K (φ) ≃ 0 for |φ| ≫ 1.
In the case, the diagonalization is not a demanding task, and it is better to use (3.19) in estimating (3.6) for large N (small lattice spacings for fixed physical lattice size). The cost of computing correlation functions reduces in the similar manner. Although one-point functions are evaluated in the same way as correlation functions shown in section 3.2, we may use another way with the numerical derivative. For given local bosonic operator O(t), let us consider (3.20) The partition function Z ′ with the action S ′ is also expressed as the form (3.6) with modified T ′ and S ′ . We have 21) and the one-point local function can be estimated without using an impurity matrix. Also in the case, the cost reduces by diagonalizing S ′ and T ′ .

Numerical results
We have proposed a method of computing the partition function and correlators based on the transfer matrices in the previous section. In this section, we demonstrate our method for φ 6 -theory and Scarf I potential for which supersymmetry is unbroken.

φ 6 -theory
Let us consider the φ 6 -interactions given by the cubic superpotential, where m is the mass and λ = g/m 2 is the dimensionless coupling. Then, any dimensionful quantities are given in a unit of the physical scale m. For instance, ma is the physical lattice spacing and βm is the period (lattice size) in the physical unit. Figure 1 shows the partition function Z P against βm at fixed am = 0.01, K = 150. At first, using the Gauss-Hermite quadrature without rescaling φ explained in section 3.3, we use (3.6) with the transfer matrices (3.4) and (3.5) to compute Z P . We also plot Z AP , which is the partition function with anti-periodic boundary condition, evaluated from Z AP = tr(T N ) + tr(S N ). Since one exact supersymmetry leads to the correct Witten index even at finite lattice spacing, we find that Z P = 1 within the errors of O(10 −9 ) for all βm. This is because S only has the eigenvalue one and the other eigenvalues coincide with those of T . We also confirm that tr(S N ) and tr(T N ) rapidly converge as βm increases. So, in the following, we take βm = 30 which is large enough for the purpose of this section since the finite β effects are of the order of e −βm . The Gauss-Hermite quadrature with the scaling is employed in the computation. We have to choose the number of Gauss-Hermite points K and the scale parameter s for each ma. Although the approximation is expected to be improved as K increases, we find that K ≃ 150 − 200 is large enough to show the convergence of results, for instance, as seen in Figure 2. For fixed K, we optimize s such that the relative error δ = |Z P − 1| is smaller than O(10 −9 ). Figure 3 is an example of calculation for determining s. In this case, we have δ O(10 −9 ) for 0.21 s 0.24. Table 1 shows the parameters used to compute the energy eigenvalues and correlators of φ 6 theory, chosen by the procedures above. We concentrate on the case of This value of λ has been used in the previous Monte-Carlo studies. The continuum limit is achieved by taking ma → 0 for fixed βm. In the table, we basically choose N such that the lattice spacings are given in equal intervals as ma = (βm)/N . 7 Same as the free theory, the scale parameter s approaches zero as ma decreases. Although we could not find s that realizes Z P = 1 within errors of O(10 −9 ) for ma < 0.001 within K ≤ 200, the lattice spacings are smaller than those used in the Monte-Carlo studies. The energy eigenvalues can be read directly from the transfer matrices S and T . As already mentioned in section 2.2, S and T correspond to two Hamiltonians in SUSY QM, H B and H F . We can read the eigenvalues of bosonic states from S = exp(−aH B ) and ones of fermionic states from T = exp(−aH F ), and obtain E B i and E F i by diagonalizing the transfer matrices S and T . We find that the lattice results satisfy (2.21) within the relative errors of O(10 −9 ) in all parameters (finite lattice spacings), at least, for five smallest eigenvalues. This is because the one supersymmetry remains exactly on the lattice. Figure 4 shows the five smallest non-zero energy eigenvalues for several lattice spacings. To extrapolate the values to the continuum limit, we use a fit formula, Note that a 0 is the value of E/m at continuum limit. Table 2 shows the fit results. The main results are obtained from a fit for ten data points at small lattice spacings, while the errors are given by the largest difference between the main fit and other fits; one for   Table 2. Energy eigenvalues of φ 6 theory. The lowest five eigenvalues are obtained by using a fit, E/m = a 0 + a 1 (ma) + a 2 (ma) 2 , with a 0 the values of E/m at the continuum limit.
twenty points, and fits using cubic polynomial for ten and twenty points. As explained above, E i ≡ E B i = E F i and E B 0 = 0 within the relative errors of O(10 −9 ) for all lattice spacings. So the values given in this table are common for E B i and E F i because the errors of the extrapolation are larger than O(10 −9 ). The extrapolated values of E/m can be determined in very high precision, within the errors less than 0.001%.
In Figure 5, the correlation functions φ N φ t and ψ t ψ N are shown as the black points and the red ones, respectively, in the case of am = 0.01. Those points look like two curves since they are estimated for all t without statistical errors. The correlators clearly show the exponential damping, and the slope is expected to be E 1 /m. One can obtain any correlators in the same manner as φ N φ t and ψ t ψ N , namely, as very clear signals since there are no statistical uncertainties in our method. Although the energy eigenvalues can be read from the transfer matrices and were already shown above, let us apply the standard techniques to extract the effective masses from the correlation functions for the purpose of making a comparison our approach and Monte-Carlo method.  Figure 6 shows the effective masses defined by which are estimated for ma = 0.01 (N = 3000). The signals are also very clear without any statistical errors in comparison with the Monte-Carlo results. They are not constants since the high energy modes contribute near t = 0 and the finite size effects are seen around t = N/2. The degenerated plateaus are seen for 0 ≪ t ≪ N/2. Estimating the effective masses at t ≃ N/5 = 600, we have m fermi /m = 1.6678215773620(2) and m boson /m = 1.667821577 (1). The errors are estimated by varying t from 500 to 700. This degeneracy is realized in the accuracy of nine digits. In Figure 7, the finite size effect is shown by changing βm with the other parameters fixed. The masses are estimated at t = [N/5]. 8 We find that βm = 30 is large enough to obtain m eff having no finite size effects.
In Figure 8, an extrapolation to the continuum limit are shown for the effective masses estimated at t = [N/5]. We simply use the polynomials f (x) = a 0 + a 1 x + a 2 x 2 as a fit function. The systematic errors of the fit are estimated in the same manner as the extrapolation of energy eigenvalues. The resultant effective mass is m eff /m = 1.686497(3) which is completely same as E 1 /m shown in Table 2.
Finally we show the numerical result of the SUSY Ward identity in the present lattice model. The SUSY Ward identity is used to show whether the broken supersymmetry is restored as a symmetry in the continuum limit. Consider the lattice supersymmetry transformation of ψ N φ t − ψ t φ N . Then one find that where . We can test whether full N = 2 supersymmetry are restored in the continuum limit by examining R (1) (t) and R (2) (t) because δ(ψ N φ t − ψ t φ N ) vanishes in the continuum theory. Figure 12 shows numerical results of SUSY Ward identity for free theory (λ = 0) and interacting case (λ = 1). In the free theory, we observe that R (1) (t) = R (2) (t) = 0 within the error of O(10 −9 ) for all t since full N = 2 supersymmetry remain on the lattice. For interacting case, although R (1) (t) vanishes for all t, R (2) (t) does not. We found that it vanishes for 1 ≪ t/a ≪ N , and SUSY breaking effect does not survive at a scale lower than the cut-off. Thus we find that full N = 2 supersymmetry are restored in the continuum limit even for interacting case.

Scarf I potential
The Scarf I model is characterized by the superpotential, where α is a parameter with the mass dimension 1/2 and λ is the dimensionless coupling. This potential is one of SIPs, and the energy spectra are exactly solved as (4.10) We now take λ = 10, βα 2 = 5, (4.11) to test our method in strong coupling region λ ≫ 1. One find that the finite size effect is small enough by taking βα 2 = 5, since it is of the order of e −βE 1 and βE 1 ≃ 10βα 2 . Table 3 shows our parameters used in the computations of Scarf I model. We determine them using the same prescription as that of φ 6 theory. We also use the Gauss-Legendre quadrature which is more effective than Gauss-Hermite method for small lattice spacings in the sense that Z P = 1 for smaller K by choosing the scale parameter s.
In Figure 10, the energy spectra obtained by our method are shown. Same as φ 6theory, the eigenvalues of T and S satisfy (2.21) within the errors of O(10 −9 ).   shows the results of extrapolation to the continuum limit. We use a linear function to fit five data points at small lattice sizes, and the errors are estimated from the maximum of differences between the fit and ones for ten points and using quadratic polynomial. The obtained a 0 coincide with the exact solution within the small errors less than 0.001%. We also estimate the effective masses from correlators φ N φ t and ψ N ψ t by the standard techniques used in the Monte-Carlo studies. The detailed procedures were already explained in φ 6 -theory. Figure 11 shows the extrapolation of effective masses to the continuum limit. We employ f (x) = a 0 + a 1 x as a fit formula. The main fit and error estimation are performed in the same manner as energy eigenvalues in Table 4. Since the values at   Figure 10. Energy spectra against the lattice spacing for Scarf I model. The lowest five eigenvalues are estimated for each α 2 a at fixed λ = 10, βα 2 = 5. the continuum limit also reproduce the exact solutions with accuracy of 99.999%, we can say that the correlators are evaluated in high precision by our method. SUSY Ward identity is shown in Figure 12. R (1) (t) vanishes for all t thanks to the exact supersymmetry, while R (2) (t) vanishes for 1 ≪ t ≪ N as with the case of φ 6 -theory. Since two identities hold at a scale lower than the cut-off, full supersymmetry is restored in the continuum limit. In this way, one can confirm the supersymmetry restoration with high precision, compared to the Monte-Carlo method.

Conclusion and discussion
We have proposed a computational method based on the transfer matrix representation of lattice quantum mechanics. The numerical quadrature has been used to define the transfer matrices as finite dimensional ones, and Gaussian quadrature with rescaling variable went well even for interacting cases. We have tested our method for cubic and Scarf I potentials, and found that the energy eigenvalues and correlators can be evaluated in high precision within reasonable number of discretized points K. Compared to Monte-Carlo method, stochastic processes are not needed to estimate the expectation values, and the sign problem does not exist in this method. The systematic error from finite-K effect is well-controlled by taking larger K and tuning the scale parameter s. The computational cost grows up with the logarithm of lattice size. Thus one can perform the computations at very small lattice spacings by taking large lattice sizes. Our method is useful for improving lattice SUSY models since one can examine the SUSY Ward identity and determine the cut-off dependence of physical quantities precisely.
Although we demonstrate our method in a specific lattice model of SUSY QM, the main idea of this paper is not limited to that case. Also, for non-SUSY models, a finite dimensional transfer matrix is defined by the numerical quadrature, and one can improve the method by rescaling variables. The scale parameter s is chosen such that the Witten index reproduces the correct value for the present lattice SUSY model. However, a method of tuning s is unknown for general models and remains an open question.
The higher dimensional theories could be studied by extending our method since the transfer matrix approach is related with not only TNR but also worldlines and worldsheets representation of SU (N ) gauge theory with fermions [49][50][51]. These kinds of extensions are applicable to higher dimensional models, and the techniques established in this paper could be useful.