Tensor network formulation for two-dimensional lattice $\mathcal{N}=1$ Wess-Zumino model

Supersymmetric models with spontaneous supersymmetry breaking suffer from the notorious sign problem in stochastic approaches. By contrast, the tensor network approaches do not have such a problem since they are based on deterministic procedures. In this work, we present a tensor network formulation of the two-dimensional lattice $\mathcal{N}=1$ Wess-Zumino model while showing that numerical results agree with the exact solutions for the free case.


Introduction
Supersymmetric field theories have attracted great attention because they provide a deep insight about the non-perturbative physics [1][2][3] and have a close relation with the gravitational theory [4]. The lattice simulations are promising approaches to obtain a further understanding of them. However, it is generally difficult to use the standard Monte Carlo techniques for the lattice supersymmetric theories on account of the sign problem, and the theories with the supersymmetry breaking may be the most difficult cases as suggested from the vanishing Witten index [5]. In this paper, we apply the tensor network approach, which is free of the sign problem, to the two-dimensional lattice N = 1 Wess-Zumino model in order to make a breakthrough on the issue.
The two-dimensional N = 1 Wess-Zumino model is a supersymmetric theory in which a real scalar interacts with a Majorana fermion via the Yukawa term originate from the superpotential [6]. The supersymmetry is spontaneously broken for the supersymmetric φ 4 theory in a finite volume [5], and the Witten index becomes zero because the fermion Pfaffian has both the positive and negative signs. For the infinite volume case, the absence of the non-renormalization theorem suggests that the breaking may occur even at the perturbative level [7,8], and the theory has a rich phase structure which should be clarified by numerical methods free from the sign problem. tensor network representation for the fermion part and the boson part are individually constructed. By combining those two results, the tensor network representation for the total partition function is also given. Section 4 shows the numerical results for the free case, and we compare them with the exact ones. A summary and a future outlook are given in section 5.

Continuum theory
Two-dimensional N = 1 Wess-Zumino model is a supersymmetric theory that consists of a real scalar field φ (x) and a Majorana fermion field ψ (x). In the Euclidean space-time, the corresponding action is given by where γ µ is the gamma matrix which satisfies The Lorentz index µ takes two values 1 or 2, and the Einstein summation convention is used throughout this paper. Showing the indices in the spinor space explicitly, γ µ and ψ (x) are written as (γ µ ) αβ and ψ α (x) for α, β = 1, 2. The spinor index α and the space-time coordinate x are often suppressed without notice. W (φ) is an arbitrary real function of φ, which is referred to as the superpotential in the superfield formalism, and gives the Yukawaand φ n -type interactions with common coupling constants. W ′ (φ) is the first differential of W (φ) with respect to φ, that is, W ′ (φ) ≡ (d/dφ) W (φ). The Majorana fermion ψ satisfiesψ where C is the charge conjugation matrix which obeys For any W (φ), the action in eq. (2.1) is invariant under the supersymmetry transformation δφ (x) =ǭψ (x) , (2.5) where ǫ is a global Grassmann parameter with two components andǭ satisfies eq. (2.3).

Lattice theory
Let us consider a two-dimensional square lattice with the lattice spacing a and the volume V = aN 1 × aN 2 , where N 1 , N 2 ∈ N. In this paper, a is set to unity, and the lattice sites are simply expressed by integers: Γ = { (n 1 , n 2 ) | n µ = 1, 2, . . . , N µ for µ = 1, 2 } . (2.7) All of the fields live on the lattice sites n ∈ Γ and satisfy the periodic boundary conditions in both directions. The forward and the backward difference operators, ∂ µ and ∂ * µ , are given by whereμ is the unit vector along the µ-direction, and the symmetric difference operator is given by ∂ S µ = ∂ µ + ∂ * µ /2. We define the lattice Wess-Zumino model according to ref. [18]: where the lattice Dirac operator D which acts as Dψ m = D mn ψ n is given by with the nonzero real Wilson parameter r. In the following, S B denotes the pure boson part of the action: (2.12) Note that the kinetic term of φ is given by the symmetric difference operator instead of the forward one in the naive boson action and an extra Wilson term is included in the boson sector. In this paper we refer to the bosons with the Wilson term as the Wilson bosons in the same sense as the Wilson fermions. While S B,naive has only the nearest-neighbor interactions, S B has the next-nearest-neighbor ones that cause difficulties in constructing the tensor network representation of the partition function. This point will be discussed later.
In the free theory with the action in eq. (2.10) is invariant under a lattice version of the supersymmetry transformation even at a finite lattice spacing because S B has the similar structure with the Wilson-Dirac operator D in eq. (2.11) in contrast to the naive one. For the interacting cases, however, the invariance is explicitly broken owing to the lack of the Leibniz rule for the lattice difference operators. The broken supersymmetry is shown to be restored in the continuum limit, at least, at all orders of the perturbation [18]. The associated partition function is defined in the usual manner: with the path integral measures Here dψ n,α is a measure of the Grassmann integral defined in the following. The Grassmann variable ξ i and its measure dξ i (i = 1, . . . , I) satisfy The Grassmann integral is then defined by which suggests that dξ i is equivalent to ∂/∂ξ i . In the free theory, the boson and the fermion are decoupled from each other, and the respective partition functions are given by where p µ = 2πn/N µ (n = 0, 1, 2, . . . , N µ − 1) and the product in eq. (2.22) is taken for all possible momenta [35]. Note that Z B = ∞ (Z F = 0) for m = 0, −2r, −4r when N µ is an even integer because the first term and the second term in the square root in eq. (2.22) simultaneously vanish for certain combinations of p 1 and p 2 . Thus we find that the Witten index, which is defined as the partition function with periodic boundary conditions in a finite volume, reproduces the continuum one, sign {m}, for |m| ≪ 1.
After integrating the fermion field, the partition function can also be written as where the Pfaffian of a 2I × 2I anti-symmetric matrix A is defined by for Grassmann variables {ξ i } and corresponding measures {dξ i }. The fermion Pfaffian Pf (C * D) flips its sign depending on the scalar field in the interacting cases. To overcome this sign problem we employ the TRG method, whose first step is to represent eq. (2.17) as a network of uniform tensors, which is explained in the next section.
3 Tensor network representation of partition function

Fermion Pfaffian
We construct a tensor network representation for the fermion part of eq. (2.17) which yields the Pfaffian after integrating the fermion field as found in eq. (2.25). The basic idea follows from refs. [32,34] which deal with the Dirac fermions. We describe the procedure for the Majorana fermions with any value of the Wilson parameter r. Now we use the following representations for γ µ and C that satisfy eqs. (2.2) and (2.4): where σ i is the standard Pauli matrix. The method presented in this section is applicable to any possible choice of γ µ and C, and they just lead to different tensors. Then the Majorana spinor takes the form ψ n = ψ n,1 ψ n,2 ,ψ n = ψ n,2 , −ψ n,1 , and we obtain − 1 2 n∈Γψ n Dψ n = n∈Γ 1 + r 2 ψ n+1,2ψ n,1 + ψ n+2,2 ψ n,1 which are local transformations of the field variable ψ n .ψ n,α is introduced only to write eq. (3.4) as simple as possible. Note that the second term in eq. (3.7) We will see that u n , v n , p n , q n , which take 0 or 1 because of the nilpotency of ψ n,α (and ψ n,α ), are regarded as the indices of tensors. The four types of hopping factors have the same structure as Ψ n+μ Φ n , where Ψ n+μ and Φ n are single-component Grassmann numbers.
It is straightforward to show where new independent Grassmann numbers θ n ,θ n+μ and the corresponding measures dθ n , dθ n+μ satisfy eqs. (2.20) and (2.21) with the periodic boundary conditions. By applying this identity to each hopping factor in eq. (3.7) individually, one can make a tensor network representation.
Then the fermion part of the partition function is represented as a product of tensors where ξ n ,ξ n , χ n ,χ n , η n ,η n , ζ n ,ζ n , and those with bars are single-component Grassmann numbers introduced in the manner of eq.
for all possible indices with single-component Grassmann numbers Ψ, Φ, By integrating Ψ and Φ by hand, we can obtain the tensor elements. In the case of r = 1, note that the indices v n and q n vanish and that the Grassmann fields χ n and ζ n are decoupled because the second term in the RHS of eq. (3.4) is absent. In that case, the tensor network representation becomes much simpler: where It is rather straightforward to show that eq. (3.7) is reproduced from eq. (3.9) with eqs. (3.10) and (3.11) and from the identity in eq. (3.8). We now note that the eight Grassmann measures in the RHS of eq. (3.10) should be in this order and that the set of measures at the site n commutes with ones at different lattice sites because they are Grassmann-even as a set for non-zero elements of the tensor given in eq. (3.11).
The indices x n ≡ (u n , v n ) and the Grassmann fields ξ n , χ n carry the information of the hopping factors with µ = 1 as indicated by the last factors in eq. (3.9) while t n ≡ (p n , q n ) and η n , ζ n are related to the hopping with µ = 2. In this sense, x n , t n , x n−1 , t n−2 , which are the indices of the tensor in eq. (3.9), can be interpreted as being defined on the four links which stem from the site n. Since each index is shared by two tensors which are placed on the nearest-neighbor lattice sites (see eq. (3.9)), we can find that the partition function Z F is expressed as a network of the tensor T Fx ntnx n−1 t n−2 on the two-dimensional square lattice Γ.
If one uses another representation of γ µ and C, then the same partition function is given by a different tensor. This means that the tensor network representation is not uniquely determined.

Boson partition function
The tensor network representation is also constructed for the pure boson part of eq. (2.17) with S B in eq. (2.12). It is, however, not straightforward to construct a simple representation because S B has the next-nearest-neighbor interactions and φ is a non-compact field. A popular way to avoid the former issue is to rewrite S B in a nearest-neighbor form with the aid of auxiliary fields. For the latter, we employ a new method using a discretization for the integrals of φ. 2 After these procedures, we find that a discretized version of eq. (3.14) can be expressed as a tensor network for arbitrary discretization schemes.
Since the formulation is actually irrelevant to the details of the scalar theory, we will derive a tensor network for a general theory: where Dϕ = ∞ −∞ n∈Γ dϕ n,1 dϕ n,2 · · · dϕ n,N . We assume thatS B (ϕ) is invariant under the PT-transformation on a two-dimensional square lattice and has the interactions up to the nearest-neighbor, and that ϕ n is a non-compact real field with N components. As seen in section 3.2.5, it is very easy to extend it to the non PT-symmetric case.
We will show that eq. (3.14) can be expressed in the form of eq. (3.15) with N = 3 in section 3.2.1. After decomposing the hopping terms ofS B in section 3.2.2 and introducing a formal discretization for the integrals of ϕ in section 3.2.3, we give the tensor network representation for a discretized version of eq. (3.15) in section 3.2.4.

Introduction of auxiliary fields
The boson action S B in eq. (2.12) is transformed into a nearest-neighbor form using two real auxiliary fields G and H: with S B,naive given in eq. (2.13), α = (1 − 2r 2 )/2, and β = 1/ √ 2. Note that α is real for |r| ≤ 1/ √ 2 but becomes a pure imaginary for |r| > 1/ √ 2. The integral measures for G n and H n are defined in exactly the same way as φ n in eq. (2.18). Although, in general, two auxiliary fields are necessary for the next-nearest-neighbor interactions in two directions, it is somewhat surprising to find that G is decoupled from the other fields for particular values r = ±1/ √ 2, and the required auxiliary field turns out to be only H. It is clear thatS B has only the on-site and the nearest-neighbor interactions which are invariant under the PT-transformation Defining a three-component field variable we find that eq. (3.16) is just eq. (3.15) with N = 3.

Symmetric property of local Boltzmann weight
In the previous section 3.2.1, we found that eq. (3.14) is a special case of eq. (3.15). Hereafter we will try to derive a tensor network representation of general one (3.15). Before that, let us see the hopping structure of the local Boltzmann weight, which is an important building block of the tensor as shown in following sections 3.2.3 and 3.2.4. It can be easily shown thatS B is expressed as where L µ is symmetric in the sense that L µ (ϕ, ϕ ′ ) = L µ (ϕ ′ , ϕ) which is a consequence of the PT-invariance of the action. 3 All of the hopping terms with respect to the µ-direction are in L µ (ϕ n , ϕ n+μ ). This decomposition is actually not unique because the positions of the on-site interactions and some constants are free to choose. For our case in eq. (3.17), we find Note that βH n and βH m have the different signs for µ = 2.
The Boltzmann factor e −S B can be written as which is symmetric in the same sense as that of L µ . This symmetric property plays an important role in the subsequent discussion.

Discretization of non-compact field
The non-compactness of the variable ϕ is cumbersome in extracting the tensor structure from f µ (ϕ n , ϕ n+μ ) in practice. There are several possible ways to make the indices of the tensor. In our method, we first carry out a discretization of the variable ϕ itself, which automatically makes the partition function in eq. (3.15) into a discretized form. 3 We can express the action asSB = n∈Γ 2 µ=1 Kµ (ϕn, ϕ n+μ ) using a trial choice of Kµ. Actually, Kµ (ϕn, ϕ n+μ ) transforms to Kµ (ϕ−n, ϕ −n−μ ) by the PT-transformation, and the PT-invariance of the action tells us thatSB = n∈Γ 2 µ=1 Kµ (ϕ n+μ , ϕn). Thus the symmetric Lµ is always defined as To make the discussion of the discretization clearly understood, let us begin with a one-dimensional integral which converges for a given function f (x). We can formally approximate this integral with a discretized form for which I = lim K→∞ I (K) is simply assumed. K is a parameter to control the approximation of the integral by the sum. Now we suppose that S K is a set containing K numbers, x 1 , x 2 , . . . , x K , which are given by a discretization scheme (disc.), and that The Gauss-Hermite quadrature gives a concrete example of this abstract definition. The RHS of eq. (3.25) is then defined as follows: Here x i (i = 1, . . . , K) is the i-th root of the K-th Hermite polynomial, and w i are the weights given by the Hermite polynomial and x i . The RHS of eq. (3.26) has an extra exponential function because this quadrature is designed so that f (x) which has a damping factor e −x 2 is well approximated. In this case, we find that S K is a set of the roots and the weight w i e x 2 i is the ingredient of (GH) . For a well-behaved f (x), one can expect that I = lim K→∞ I (K).
With the prescriptions above, eq. (3.15) can be discretized as by replacing the measures for ϕ n by . 4 5 Note that we use the same discretization scheme for all components of ϕ n . Here eq. (3.22) is also used, and K is the number of discrete points. It is found that f µ (ϕ n , ϕ n+μ ) is a matrix whose indices are ϕ n and ϕ n+μ which take the K N discrete numbers in S N K . In this way, we now consider f µ as a matrix, and this fact provides a benefit for a numerical treatment; that is, one can use linear-algebra techniques instead of the functional analysis. The indices of the tensor will be naturally derived from this matrix structure of f µ as will be seen in section 3.2.4.

Construction of tensor
In order to derive the tensor network structure from eq. (3.27), one needs to separate ϕ n and ϕ n+μ in f µ . If this separation works, the original field ϕ n can be traced out at each n.
Since f µ is a symmetric matrix with complex entries in general, which is found in the previous sections 3.2.2 and 3.2.3, we carry out the Takagi factorization: for ϕ, ϕ ′ ∈ S N K , where U and V are unitary matrices, U T and V T are the transposes of U and V , respectively, and σ w and ρ s are non-negative. Note that this factorization depends on the discretization scheme, which determines the set S K . Instead of the Takagi factorization, we can also use the SVD as seen in the next section. We thus find that eq. (3.27) is written as where for all indices. One can verify eq. (3.30) from eq. (3.27) by applying the factorization in eqs. (3.28) and (3.29) to f µ (ϕ n , ϕ n+μ ) for µ = 1, 2 with the local indices w n , s n . Then the index w n (s n ) can be interpreted as a variable defined on the link which connects n and n +1 (n +2), so eq. (3.30) forms a tensor network on the two-dimensional lattice Γ as with the case of the fermion partition function in eq. (3.9). Here one finds the correspondence between the tensor indices and the hopping structure of the lattice action as in the fermion part. From this one can see that the tensor network structure is originated from the kinetic terms for both fermions and bosons. We expect that, in the large K limit, Z B (K) converges to Z B with an exact tensor network representation if we can find a proper discretization scheme so that T B (K) converges to T B in K → ∞. In practice one has to confirm that Z B (K) converges to Z B with increasing K in the choice of a discretization scheme. We will see this point in section 4.3.

Miscellaneous remarks
We give some miscellaneous remarks which may be important for future applications and deeper understanding of the symmetry of the tensor network. The tensor network representation in the form of eq. (3.32), which gives the boson partition function in eq. (3.14), is not uniquely determined. Let F and G be regular K × K matrices. We then find that eq. (3.32) also holds for another uniform tensorT B given bỹ for all indices. Furthermore, by using F n and G n that are regular matrices satisfying the periodic boundary conditions on the two-dimensional lattice Γ and transforming T Bw nsnw n−1 s n−2 by F n , G n , F −1 n−1 , and G −1 n−2 , Z B can also be written in terms of the non-uniform tensors. This means that the tensor network representation of the partition function is invariant under the gauge transformations for tensors.
The expression of eq. (3.30) is rather general in the sense that we can always find it for two-dimensional PT-invariant theories with the real scalars. It is very easy to generalize this result to more complicated cases, non PT-invariant actions which, for example, have only one of φ 2 n+μ φ n or φ 2 n−μ φ n terms or the theories with the complex scalars. For those theories, although f µ is not symmetric in general, we can use the SVD instead of the Takagi factorization. Then, U T and V T in eqs. (3.28) and (3.29) are replaced by other unitary matrices, and we can express the partition function by a similar construction of the tensor to eq. (3.31), where the second U and the second V are replaced with the other ones. An extension to the higher-dimensional theories is also straightforward.
We have much simpler expressions for the cases of S B,naive given in eq. (2.13) because the auxiliary fields are not needed (N = 1) and L µ is isotropic and given by a single L:

Equation (3.31) then becomes
because f µ ≡ f = e −L for µ = 1, 2 and because V = U and ρ i = σ i in eqs. (3.28) and (3.29). In this case, instead of the Takagi factorization, we can use the SVD: where O and P are real symmetric matrices. Then we have

Total tensor network
We have seen that the fermion and the boson partition functions can be expressed as the tensor networks in the previous two sections. By combining these results, we can also express the total partition function as a tensor network. Before presenting the total tensor, let us introduce combined indices X n , T n . X n is defined as X n = (u n , v n , w n ), where (u n , v n ) and w n are indices of the fermion and the boson tensors, respectively. T n is also defined as T n = (p n , q n , s n ), and the dimension of X n and T n is 2 × 2 × K N .
The total tensor is made by replacing e −S B (φ) in eq. (3.14) with e −S B (φ) Z F (φ) and repeating the same procedure for making the tensor network representation of the boson partition function. Additional contributions by Z F do not give any complexity. We find that the total tensor network representation is given by the boson one in eq. (3.30) multiplied by Z F (φ) from the right: where U , V , σ w , and ρ s are given by eqs. (3.28) and (3.29). The measure DΞ uvpq is given in eq. (3.10), ξ,ξ, χ,χ, η,η, ζ,ζ are one-component Grassmann numbers, and T F is the tensor for the fermion part defined in eq. (3.11). Note that T F (φ n ) includes only φ n which is a component of ϕ n . The total tensor T (K) XT X ′ T ′ is uniformly defined on the lattice. Now the original partition function Z is expressed as a tensor network Z (K). We have built it for a general superpotential by focusing on the hopping structure of the lattice action. Introduction of local interaction terms does not change our formulation but rather elements of tensor. Moreover, the same structure of the tensor network leads to the same order of computational complexity. 6 We will numerically verify that Z (K) indeed converges to Z by using the TRG as a coarse-graining scheme for the tensor network in the next section.

Numerical test in free theory
The partition function of the lattice N = 1 Wess-Zumino model has been expressed as a tensor network in eq. (3.38). In this section, we test the expressions in the free theory given by eq. (2.14) varying the mass for three lattice sizes V = 2 × 2, 8 × 8, 32 × 32 with the periodic boundary conditions. Numerical tests in the free theory are effective to study whether the tensor is correctly given by our new formulation because the tensor network structure is derived from the hopping terms in the action, i.e. the kinetic terms. The computation is performed with the value of the Wilson parameter r = 1/ √ 2 to reduce the computational cost because the auxiliary field G is decoupled as seen in section 3.2.1.

Some details
In sections 4.2 and 4.3, we compute Z F and Z B individually using the (Grassmann) TRG since they are independent with each other in the free theory. In section 4.3, the Witten index given by the total partition function Z is computed by the Grassmann TRG. Since the free theory is exactly solvable, we can compare an obtained result X TRG with the exact solution X exact by computing In what follows, we briefly describe the TRG while introducing D cut which defines the truncated dimension of tensors. The SVD allows us to express a tensor T ijkl (i, j, k, l = 1, 2, · · · , N ) of which the tensor network representation of a partition function Z is made as T ijkl = N 2 I=1 S ijI σ I (V † ) Ikl , where S and V are unitary matrices and σ I is the singular value of T ijkl . We assume that the singular values are sorted in descending order: σ 1 ≥ σ 2 ≥ σ 3 ≥ · · · σ N 2 ≥ 0. 7 In the TRG, T ijkl is approximately decomposed: where D cut , which is fixed throughout a computation, is used to truncate the dimension of the tensor indices if it is smaller than N 2 . If not so, the summation in eq. (4.2) is done up to N 2 without the truncation. A similar decomposition can be done with a different combination of the indices: The coarse-grained tensor T new IJKL with I, J, K, L = 1, . . . , min{D cut , N 2 } is then given by contracting the rank-three tensors √ σS, √ σ ′ S ′ , √ σV, √ σ ′ V ′ and forms a network again as with T ijkl . We can compute the partition function Z by repeating this procedure. Since the number of tensors decreases through the coarse-graining, Z is finally given by a single tensor for which the indices are contracted: Z = Dcut I,J=1 T new IJIJ . More details are shown in ref. [37], and appendix A is given for the Grassmann cases.
We employ the Gauss-Hermite quadrature (3.26) to discretize the integrals of φ and H in (3.22): for eq. (2.14) and r = 1/ √ 2. The two-dimensional variable ϕ n = (ϕ n,1 , ϕ n,2 ) = φ n / √ 2π, H n / √ 2π is again used for the notational simplicity. S K is a set of the roots of the K-th Hermite polynomial. We use the SVD to decompose f µ , which are K 2 ×K 2 real symmetric matrices, as For reducing the memory usage and the computational cost, we initially approximate the tensor network representation of eq. (4.4) by D init ≤ K 2 : Note that D init defines the bond dimension of the initial tensor. We will simply take D init = D cut for evaluating Z B in section 4.3 and D init = D cut /2 for the Witten index in section 4.4 because the bond dimension does not change after the coarse-graining steps under these choices.
Here we mention the computational costs for the coarse-graining of tensor networks and for the construction of tensors. Both of them are mainly consists of the SVD and the contraction of tensor indices. Since the cost of the numerical SVD for square matrices is proportional to the third power of the matrix dimension, the computational effort required for the numerical decomposition described in eq. (4.2) and in eqs. (4.6) and (4.7) are in proportion to N 6 and K 6 , respectively. A contraction of tensor indices is expressed as a summation of them, so the cost of the contraction depends on the number of the tensor indices. Then it is proportional to D cut 6 when contracting the rank-three tensors described around eqs. (4.2) and (4.3), and is proportional to K 2 × D init 4 when building the tensor in eq. (4.9). For the coarse-graining step, one can find that the volume-dependence of the cost is milder than D cut -dependence as follows. Since the TRG is a coarse-graining of spacetime, one can reach a large space-time volume by simply iterating the same local blocking procedures. More directly, the computational cost of the TRG is proportional to the logarithm of the space-time volume, i.e. the number of iterations. Summarizing the above, the computational cost for the coarse-graining of tensor networks is proportional to D cut 6 ×ln V , and that for the construction of tensors is proportional to max K 6 , K 2 × D cut 4 , where N = D init = D cut is simply assumed. The purple curves represent the exact solutions given by eq. (2.23). Three negative peaks at m = 0, − √ 2, −2 √ 2 correspond to the fermion zero modes, and the exact Pfaffian has the negative sign for −2 √ 2 < m < 0 as can be seen in eq. (2.23). In the top plot of figure 1, the green symbols (D cut = 8) around the peak at the center are rather deviated from the exact solution, and they even have the opposite sign. The deviation becomes smaller as D cut increases, and the yellow symbols (D cut = 16) have the correct sign and agree well with the exact one even near the peak. The situation is further improved by taking larger volumes even for the smallest D cut , and the numerical results fit well with the analytical curve in the center and the bottom figures.

Free Majorana-Wilson fermion
These observation can also be clearly understood in figure 2, which shows the relative errors δ (ln |Z F |) given by eq. (4.1). Note that the case for D cut = 16 on V = 2 × 2 have extremely small errors. This is because the maximal bond dimension of the coarse-grained tensors on V = 2 × 2 lattice is less than or equal to D cut . In other words, no truncation occurs in the TRG steps. This striking feature is only found in the pure fermion case. In contrast, the discretization error and the truncation error are inevitable in the boson case since the approximation already enters in deriving the tensor network representation of the boson partition function, and furthermore the tensor indices are truncated to carry out the numerical evaluation as seen in previous section. For all volumes used in the computation, the relative errors almost monotonically decreases as D cut increases.
Thus we can conclude that the Pfaffian with the correct sign is reproduced from the tensor network representation in eq. (3.9) with eq. (3.11) using the Grassmann TRG within tiny errors O 10 −3 for physically important parameters, |m| ≪ 1, and larger volumes.

Free Wilson boson
The boson partition function is given as a discretized form Z B (K) in eq. (4.4) by applying the Gauss-Hermite quadrature to the integrals of φ and H. Then K is the number of the discrete points. We prepare the initial tensor network approximately as eq. (4.8) and compute it using the TRG for m > 0 because the adopted quadrature does not effectively work for m < 0 (we will see this point later.). It is, however, sufficient to study the case of m > 0 because the boson action does not depend on the sign of m, but on m 2 , in the continuum theory.  32. Figure 5 shows that the results are systematically improved by increasing D cut as one expects. The exponential improvement may be explained as follows. Usually the singular values of the tensor are exponentially decaying; thus from a local point of view the truncation error gets exponentially smaller by increasing D cut . Since the free energy consists of the local tensors, it is likely that its error shows such a behavior as well.
The growth of the errors is observed near m = 0. Roughly speaking, this is because the massless theory has no damping factors in f µ of eq. (4.5). We can show that f µ is expressed as (4.10) One can see that the damping factors are actually provided for m > 0 with the damping rate m 2 but is not for −4 √ 2 < m < 0 on the line φ = −φ ′ , so the quadrature does not work for m < 0. For m > 0, we have to take K larger as m decreases so that the quadrature retains effective. That structure is encoded in the initial tensor in eq. (4.9) via the matrices O, P, S, T and the singular values σ w , ρ s in eqs. (4.6) and (4.7). The singular values of the initial tensor have unclear hierarchies for small masses as seen in figure 6. Thus we find that, if m approaches zero from the right, we have to take K and D cut as large as possible to obtain the precise result. 8 The K-dependence of the relative errors is investigated in figure 7. In order to purely see the discretization effect due to finite K, we set the maximum bond dimension of the tensor K 2 and choose the lattice size V = 2× 2 that allows us to carry out a full contraction for the computation of the partition function. Although there are no other systematic errors except for finite K, the value of K is practically restricted up to 10. Figure 7 shows that the errors decrease by increasing K. From this we can say that a simple discretization scheme such as the Gauss-Hermite quadrature well approximates the original integrals if K is sufficiently large, and that the tensor network representation reproduces the correct values of the boson partition function.

Witten index of the free N = 1 Wess-Zumino model
The Witten index computed by the Grassmann TRG is shown in figure 8. Figure 9 shows the relative error of the Witten index. As discussed in section 2.2, the fermion and the boson are decoupled from each other in the free case. In this section, however, we treat the free Wess-Zumino model as a combined system of fermions and bosons; thus we perform the Grassmann TRG for a single tensor network. One can see that the results tend to converge to the exact values by increasing D cut . The obtained indices with D cut = 64 (yellow symbols) take the values near one compared with those of D cut = 32 (green symbols).
Thus we can conclude that eq. (3.38) gives a correct tensor network representation of the two-dimensional lattice N = 1 Wess-Zumino model. Z F and Z B become extremely large and extremely small, respectively, for large space-time volume. For instance, Z F are of the order of O 10 400 at m = 1 on V = 32 × 32 lattice as seen in figure 1. Surprisingly,   O (1) values are obtained as the Witten index as seen in figure 8. Namely, the boson effect balancing huge Z F is correctly reproduced using the Grassmann TRG for the total tensor.  So we can say that the TRG is a very promising approach to study the supersymmetric field theories.

Summary and outlook
We have shown that the two-dimensional lattice N = 1 Wess-Zumino model is expressed as a tensor network. The known techniques of making a tensor were refined in the fermion sector and generalized in the boson sector in the sense that it is possible to define a tensor for any way of discretizing the integrals for scalar fields. We have also tested our formulation in the free theory by estimating the Witten index and comparing it with the exact solution.
The resulting indices reproduce the exact one as D cut , the dimension of the truncated tensor indices in the TRG, increases. Now we are tackling the issue on the supersymmetry breaking by estimating correlation functions from the tensor network. Before investigating the physical breaking effects, we have to show that the artificial ones by the lattice cut-off disappears in the continuum limit beyond the arguments of the perturbation theory. We will estimate the expectation value of the action, the supersymmetric Ward-Takahashi identity, and the mass spectra of fermions and bosons to show it. We will then see the supersymmetry breaking in the model with the double-well potential by estimating several physical quantities and study the phase structure in detail.
Although we have only dealt with the Wilson type discretization of derivatives, one may use another way such as the domain wall discretization. In that case, partition functions or Green's functions will be represented as three dimensional tensor networks. For such higher dimensional tensor networks, the higher order TRG was introduced in Ref. [38], and the Grassmann version was also proposed in Ref. [39]. In this way one can in principle go this direction; however, the computational cost could be severe. Therefore further improvements of the algorithm might be needed for the actual computation in higher dimensions.
We emphasize that the methodology of constructing the tensor is given for any superpotential, that is, any interacting case, in this paper. Since the Wess-Zumino model consists of various building blocks: the scalar field, the Majorana fermion, and their interactions such as the Yukawa-and the φ 4 -interactions, we expect that our method could be very useful in TRG studies of other theories.

A Coarse-graining step in Grassmann TRG
In this appendix, we describe the coarse-graining step in the Grassmann TRG for the current boson-fermion system. We basically follow ref. [34], which deals with a pure fermionic model (the N f = 1 Gross-Neveu model) and show the method in our notation making the difference that comes from the boson part clear.   We begin with the partition function that initially takes the following form: 9 Z = {X,T } n∈Γ T XnTnX n−1 T n−2 · n∈Γ dΞ uvpq n · n∈Γ ξ n+1 ξ n un χ n+1 χ n vn η n+2 η n pn ζ n+2 ζ n qn , (A.1)  and the tensor elements are not zeros only when u n + v n + p n + q n + u n−1 + v n−1 + p n−2 + q n−2 mod 2 = 0 (A. 3) holds. The tensor T is made of the fermionic one T F (φ n ) unvnpnqnu n−1 v n−1 p n−2 q n−2 in eq. (3.11) and the boson one T B (K) wnsnw n−1 s n−2 in eq. (3.31) as in eq. (3.39). The indices u n , v n , p n , q n take two values 0 or 1 while w n , s n run from 1 to D init as seen in section 4.1. The total indices X n and T n are given by X n = (u n , v n , w n ) and T n = (p n , q n , s n ), and they run from 1 to 2 × 2 × D init . As mentioned in section 4.1, we set D init = D cut /2 for the actual computations.
The coarse-graining of a tensor network mainly consists of three steps: the SVD of tensors, a decomposition of Grassmann measures, and a contraction of the indices and taking the integrals of Grassmann variables defined on Γ. The SVD and the decomposition for Grassmann measures are performed in a different manner for even and odd sites. We will see that the coarse-grained tensors take the same form as (A.1) with v = q = 0 and are defined on the coarse-grained lattice Γ ⋆ = n + 1 2 1 +2 n = (n 1 , n 2 ) ∈ Γ, where n 1 + n 2 is an even integer. .