Grassmann higher-order tensor renormalization group approach for two-dimensional strong-coupling QCD

We present a tensor-network approach for two-dimensional strong-coupling QCD with staggered quarks at nonzero chemical potential. After integrating out the gauge fields at infinite coupling, the partition function can be written as a full contraction of a tensor network consisting of coupled local numeric and Grassmann tensors. To evaluate the partition function and to compute observables, we develop a Grassmann higher-order tensor renormalization group method, specifically tailored for this model. During the coarsening procedure, the blocking of adjacent Grassmann tensors is performed analytically, and the total number of Grassmann variables in the tensor network is reduced by a factor of two at each coarsening step. The coarse-site numeric tensors are truncated using higher-order singular value decompositions. The method is validated by comparing the partition function, the chiral condensate and the baryon density computed with the tensor method with exact analytical results on small lattices up to volumes of $4\times4$. For larger volumes, we present first tensor results for the chiral condensate as a function of the mass and volume, and observe that the chiral symmetry is not broken dynamically in two dimensions. We also present tensor results for the number density as a function of the chemical potential, which hint at a first-order phase transition.

a complex action, i.e., with a sign problem, were successfully studied, as for example the three-dimensional O(2) model with a chemical potential [10]. For systems with fermions, which are represented by Grassmann variables in the partition function, the Grassmann HOTRG (GHOTRG) was recently developed for cases where the Grassmann variables cannot be integrated out locally [11,12,13].
The aim of the current paper is to demonstrate the applicability of tensor-network methods to strongcoupling QCD with staggered quarks. The HOTRG method cannot be applied as such to strong-coupling QCD as non-local sign factors occur in the meson-baryon-loop representation of the partition function [1,2]. Clearly, this property is impossible to encode in a local tensor. To resolve this problem, we do not integrate out all Grassmann variables as in the meson-baryon-loop representation, but keep the baryonic combinations of Grassmann variables in a Grassmann tensor. Then, the partition function can be written as a full contraction of a tensor network with local numeric and Grassmann tensors. To evaluate the partition function we then apply an iterative blocking procedure, which uses ideas of the original GHOTRG method, but is specifically tailored for strong-coupling QCD.
We validate our Grassmann tensor-network method by comparing its results for the partition function, the chiral condensate and the number density with exact analytical results computed on small lattices of sizes up to 4 × 4. Then, we apply the method to larger lattices to compute the chiral condensate as a function of mass and volume, and observe that the chiral symmetry is not broken dynamically in the twodimensional case. Furthermore, we compute the number density at nonzero chemical potential, which hints at a first-order phase transition. We also briefly discuss the convergence of the tensor-network results with increasing bond dimension.
The paper is structured as follows. In Sec. 2 we reformulate the partition function of strong-coupling QCD as a tensor network of numeric and Grassmann tensors. In Sec. 3 we introduce auxiliary Grassmann variables in order to decouple the nearest-neighbor interaction terms in the different directions. Then, we discuss how the lattice can be coarsened by blocking adjacent local tensors. We present our numerical results in Sec. 4 and our conclusions in Sec. 5.

Strong-coupling QCD and its tensor formulation
In the strong-coupling limit (β → 0) of QCD, the gauge action vanishes and only the fermion action survives. For a single staggered quark field 1 with mass m, the lattice action is where x ∈ {1, . . . , V } enumerates the sites on a lattice with temporal extent L 1 , spatial extent L 2 , and volume V = L 1 L 2 . For tensor-network studies L 1 and L 2 are taken to be powers of 2. The SU(3) matrices U x,ν are defined on the links of the lattice, ψ x andψ x are 3-dimensional vectors of Grassmann variables, representing the colored quark and antiquark fields on the site x. The staggered phases are η x,1 = 1 and η x,2 = (−1) x1 , where x 1 is the time coordinate of site x. The quark chemical potential µ and an anisotropy factor γ are introduced for the Euclidean time direction. 2 To describe the system in thermal equilibrium, we use antiperiodic boundary conditions in the time direction and periodic boundary conditions in the space direction for the Grassmann variables. In the infinite-coupling limit, the SU(3) gauge fields can be exactly integrated out [1,2], giving rise to a system of mesons and non-intersecting baryon loops. For each configuration contributing to the partition function, each lattice site is assigned either to a baryon loop or to a mesonic contribution, as all Grassmann variables must be saturated, i.e., each site has to contain 3 quarks and 3 anti-quarks, in order to contribute to the partition function. The partition function can then be written as [3] where the differentials are defined as dψ x dψ x = dψ x,3 dψ x,2 dψ x,1 dψ x,1 dψ x,2 dψ x,3 = dψ x,1 dψ x,1 dψ x,2 dψ x,2 dψ x,3 dψ x,3 and 3 with mesonic combinations M x =ψ x ψ x , baryonic combinations B x = 1 3! i1i2i3 ψ x,i1 ψ x,i2 ψ x,i3 , antibaryonic combinationsB x = 1 3! i1i2i3ψx,i3ψx,i2ψx,i1 and ζ ν = γ 3 exp(±3µ) for ν = ±1 , 1 else .
To integrate out the Grassmann variables, we first expand the exponential in the mass and write the product of sums in (2) as a sum of products. Looking at a single term in the sum, i.e., a specific configuration, we observe that, for nonzero contributions, the sites have to be either baryonic or mesonic due to the Grassmann nature of the variables. Therefore the product over directions ν = 1, 2 cannot mix baryonic and mesonic contributions on a single site.
To apply Monte Carlo simulations to strong-coupling QCD, the Grassmann variables in both the mesonic and baryonic combinations are integrated out. The partition function then consists of configurations of closed, non-intersecting baryon loops with remaining sites saturated by meson contributions (including mass terms) [2,3]. A particularity of this representation is that each baryon loop contributes a multiplicative factor of (−1) to the weight of the configuration to which it belongs, coming from a reordering of the Grassmann variables along the loop when performing the Grassmann integration.
Although the mesonic part of the partition function is easily converted into a consistent tensor-network formulation, as was already shown for U(N ) [15], the baryonic contributions introduce a new problem as the baryon-loop sign factors are of a global nature and can therefore not be included in a local tensor without further ado. Therefore, HOTRG cannot be applied as such on this model. However, it turns out that this problem can be resolved using a variant of GHOTRG [11,12,13], which we specifically develop for this model.
In the following we explicitly integrate out the Grassmann variables in the mesonic combinations, but leave the baryonic ones unintegrated to avoid the generation of non-local sign factors. After integration of the mesonic Grassmann combinations, the baryonic Grassmann variables B x andB x can be regarded as fundamental (non-composite) Grassmann variables that are integrated over. This results in 4 l −1 0 1 l + 0 0 1 l − 1 0 0 Table 1: Relation between the occupation numbers l, l + and l − given by l ± = l(l ± 1)/2 and l ∈ {−1, 0, 1}.
where the set of configurations on a two-dimensional lattice of volume V is the set of all tuples k = (k 1,1 , . . . , k V,2 ) and l = (l 1,1 , . . . , l V,2 ) of mesonic and net baryonic link occupation numbers k x,ν ∈ {0, 1, 2, 3} and l x,ν ∈ {−1, 0, 1}, respectively. The occupation numbers l ± x,ν for baryons and antibaryons are mutually exclusive in (6), see also (4), and can therefore be written as functions of the net occupation number l x,ν with l ± x,ν = l x,ν (l x,ν ± 1)/2 ∈ {0, 1}, see also Table 1. The weight functions in (6) are where we introduced the notation k x,−ν ≡ k x−ν,ν and l x,−ν ≡ l x−ν,ν . For every configuration that yields a nonzero contribution to the partition function, each site x is either baryonic or mesonic: (a) baryonic site x ∈ B: All surrounding links must have k = 0, i.e., ν (k x,−ν + k x,ν ) = 0. In order to yield a nonzero contribution to the partition function when all Grassmann variables are integrated out, each baryonic site x must be occupied by exactly one factor of B x and one factor ofB x . This means that for each baryonic site x we require l − x,1 + l − x,2 + l + x,−1 + l + x,−2 = 1 and l + x,1 + l + x,2 + l − x,−1 + l − x,−2 = 1. We represent the baryon condition by (b) mesonic site x ∈ M: All surrounding links must have l = 0 and the SU(3) condition requires n x ≥ 0, see (9). This is represented by the meson condition where we use the convention Θ(0) = 1 for the Heaviside-theta function.
The configurations (index combinations) for which any single site x is neither baryonic nor mesonic, i.e., with indices such that δ x∈B = 0 and δ x∈M = 0, are not contributing to the partition function.
Note that in the partition function (6), for each contributing index configuration (k, l), the Grassmann variables appearing in the mesonic combinations have already been integrated out, and the remaining integrals only apply to the baryonic terms. Because of the baryonic condition (10), the Grassmann differentials for such a configuration can be rewritten as For each configuration contributing to the partition function (6), any link is either mesonic (k = 0, l = 0), baryonic (k = 0, l = 0) or empty (k = l = 0). Hence we can combine the mesonic and baryonic combined index occupation numbers k x,ν and l x,ν into a single combined index j x,ν of dimension 6, to reduce the total number of configurations in the partition function (i.e., we effectively reduce the number of configurations with zero weights). The relation between the combined index 0 ≤ j ≤ 5 and the mesonic and baryonic occupation numbers k and l is given in Table 2.
The partition function (6) can be written as a full contraction of a tensor network where the local tensors have a numeric and a Grassmann part, where each configuration is now characterized by its 2V indices j = (j 1,1 , . . . , j V,2 ) and we again introduce the notation j x,−ν ≡ j x−ν,ν . The local numeric tensors S (x) and the local Grassmann tensors G (x) have the following entries: where the indices k x,ν and l x,ν are implicitly defined as functions of j x,ν as in Table 2, and we recall that l ± x,ν = l x,ν (l x,ν ± 1)/2 ∈ {0, 1}, see Table 1.
Note that each configuration j in the partition function (13) selects one entry for each local tensor, and due to (14), the nonzero tensor entries of S can be classified as either "mesonic" or "baryonic". For the mesonic entries, the corresponding entries of the Grassmann tensor are all equal to 1, as all l + and l − are zero around a mesonic site. For the baryonic entries, the corresponding entries of the Grassmann tensor are non-trivial and will play a crucial role in GHOTRG. 5 As all the interaction terms in G (x) involve an even number of Grassmann variables, they are mutually commuting and so their order does not matter. Furthermore, due to the factors δ x∈B and δ x∈M in S (x) , the Grassmann tensors G (x) can be considered to be commuting with each other in (13), since for every nonzero entry of S (x) the corresponding entry of G (x) is Grassmann-even.
The local numeric tensor S (x) only depends on the site x through the staggered phase η x,ν . As η x,1 = 1 and η x,2 = (−1) x1 , there are only two different realizations of S (x) , for sites with odd and even time coordinates x 1 , respectively. 6 The bond dimension of the initial local tensor is D initial = 6, corresponding to the dimension of the index j in Table 2. During the iterative blocking procedure, the bond dimensions of the coarse-lattice tensors, which would in principle grow exponentially, are truncated to a chosen value D using HOSVD approximations [9].
In order to validate our version of the GHOTRG method, we will also investigate a simplified partition function containing only baryons. The sum over j in (13) is then restricted such that all sites are baryonic, i.e., x ∈ B for all x.

Grassmann HOTRG for strong-coupling QCD
We now explain how to evaluate the partition function given by the full contraction of the Grassmann tensor network in (13). The method can be summarized as being an iterative blocking procedure where each blocking step consists of two parts: First, new Grassmann tensors are generated on the coarse lattice, which reduces the number of Grassmann variables by a factor of two and gives rise to local sign factors. Then, an HOSVD approximation is applied to the contraction of two adjacent numeric tensors. In this process, the local sign factors are absorbed in the new numeric tensors on the coarse lattice.
The peculiarities of the strong-coupling QCD model, i.e., the use of staggered quarks and the existence of both, mesonic and baryonic contributions, require the development of a tailor-made GHOTRG.

Decoupling the Grassmann interaction terms through auxiliary variables
In order to integrate out the Grassmann variables B x andB x in the partition function (13) for one particular configuration j, satisfying x ∈ B or x ∈ M for all x, we decouple the interaction terms in different directions. This is achieved by introducing auxiliary Grassmann variables c and inserting identities of the form The interaction terms (B x B x+ν ) l + x,ν and (B xBx+ν ) l − x,ν are mutually exclusive, i.e., l + x,ν and l − x,ν cannot simultaneously be equal to one, see Table 1. Therefore, we can use the same auxiliary variable c x,ν to rewrite the two interaction terms as where each interaction term is split into two commuting factors. These identities can be applied for all x and ν = 1, 2 independently. The order of the Grassmann variables on the right hand side is chosen to facilitate the integration of B x andB x below. For later convenience we will also introduce the notation c x,−ν ≡ c x−ν,ν . After introducing the auxiliary Grassmann variables using (17), all original Grassmann variables B x and B x can be integrated out independently for different sites. To this end, we gather all eight interaction terms involving B x orB x for one particular site x, together with the differentials contained in G (x) (the commuting pairs are reordered to gather the contributions in B andB separately, such that the Grassmann integrations can be performed without generating additional sign factors). For x ∈ B or x ∈ M, satisfying the conditions (10) and (11), respectively, we obtain with B and B indicating that we only integrate over the Grassmann variables B andB. If the chosen site x is baryonic, then H (x) always contains exactly two Grassmann variables due to (10), one from each product in square brackets. On the other hand, if x is mesonic, H (x) = 1 since all l ± x,±ν = 0. In both cases H (x) is Grassmann even.
Note that to collect all interaction terms involving B x andB x in Z, resulting in (18), one also needs (17) with x → x −ν. When x is on the lower edge of the lattice in theν-direction, the interaction terms between x −ν and x will wrap around the lower edge of the lattice. These terms actually stem from the interaction terms in (13) which wrap around the lattice at its upper edge in that direction (because (13) only contains interactions from x to x +ν for x ∈ {1, . . . , V }). For these interaction terms, the variables B x andB x will be subjected to the boundary conditions in the directionν. However, we always havē B y B y =B y±Lνν B y ±Lνν , for all y, y , ν since we use antiperiodic (for ν = 1) and periodic (for ν = 2) boundary conditions. In order to preserve this property for products of original and auxiliary variables, appearing on the right hand side of (17), and to avoid explicit sign factors in the partition function, we choose the boundary conditions of the new auxiliary variables such that we always have where a y , b y are place holders for any of the original or auxiliary variables (or differentials). The conditions above are automatically satisfied by requiring the auxiliary variables to satisfy the same boundary conditions as the original Grassmann variables B andB. After integrating out B x andB x , according to (18), for all sites x in the partition function (13), we are left with 2V new Grassmann variables c x,ν and their differentials, and the partition function can be written as with numerical tensors S (x) and Grassmann tensors H (x) given in (14) and (18), respectively. Recall that the indices l x,ν of H (x) are implicit functions of j x,ν as given in Table 2. For each configuration j, the integral in (21) applies to all auxiliary fields having differentials with unit exponent. The main difference compared to the original formulation (13) for Z is that the new Grassmann variables live on the links, while the original Grassmann variables where defined on the sites of the lattice. This is crucial for deriving a consistent blocking procedure, as will be shown in the next sections.
As explained above, the Grassmann tensors H (x) in the partition function can be considered to be commuting since every entry H (x) lx,−1lx,1lx,−2lx,2 is accompanied by a factor δ x∈B or δ x∈M in the numeric tensor S (x) .
To facilitate the further manipulations, we reorder the factors in H (x) in a canonical order (chosen such that the Grassmann integrations can be more easily performed when blocking two tensors as described in further sections), with sign factor 7 ω lx,−1lx,1lx,−2lx,2 = (−1) l + which will eventually be absorbed in the numeric tensor. As l ± x,ν = l x,ν (l x,ν ± 1)/2, the products in (22) can be rewritten as such that the Grassmann tensor H (x) becomes We now introduce the notation After defining a new Grassmann tensor and absorbing the sign factor ω in a new numeric tensor the partition function (21) becomes Note that the Grassmann tensors K (x) can be considered to be commuting (Grassmann-even) in (28), since the baryonic condition δ x∈B included in (11). Hence, the Grassmann tensors can be reordered in the partition function when performing the coarsening steps discussed below, without generating additional sign factors.
In the next sections we will describe the renormalization group (RG) steps, which coarsen the lattice iteratively and halve the number of lattice sites at each iteration. The blocking of two adjacent Grassmann tensors K (x) and K (x+ν) will produce a new tensor on the coarse lattice with a Grassmann structure identical to that of the original tensors, and a sign factor that can be absorbed in the coarse-lattice numeric tensor. In the following we will call the index f x,ν ≡ f x,ν (j x,ν ) ∈ {0, 1} the Grassmann parity of the index j x,ν . The Grassmann parity f x,ν is the exponent of the Grassmann variable living on the link between x and x +ν, see (26). In the original local Grassmann tensor K (x) fx,−1fx,1fx,−2fx,2 , the indices f = f (j) are related to j by Table 2 and f = l 2 . However, in general, the index f is a function of the index j, which will be updated at each step of the blocking procedure, as will be explained in detail below.
We will see that after each RG step, the partition function will always have the shape (28), albeit with an updated numeric tensor on the coarse lattice. The RG steps are repeated until the tensor network has been reduced to a single tensor. The sum over the remaining indices of that tensor then yields the partition function Z.

Coarsening the time direction 3.2.1. Blocking adjacent tensors in the time direction
As part of GHOTRG we now discuss an RG step in the1-direction, which consists of (identical) contractions of all V /2 pairs of adjacent local tensors in that direction, where (x, x +1) denotes the pair of sites which will eventually be fused in a new coarse-grained site, and T (x,x+1) and K (x,x+1) are the new numeric and Grassmann tensors, respectively, on the coarse lattice. The integral only represents an integration over the Grassmann variable c x,1 , which is defined on the link that connects the two sites. Note that the summation variable j x,1 also appears in f x,1 = f x,1 (j x,1 ).
We first consider the Grassmann part of this contraction, which is the product The order of the two factors in the product is irrelevant as the Grassmann tensors can be considered to be commuting (as explained above), and we place them such that the Grassmann integration over the shared link can be directly performed, Note that K (x,x+1) does not depend on f x,1 due to the integration formula (16). Therefore the sum over j x,1 in (29) actually only applies to the numeric tensors, such that In terms of the blocked tensors T and K, the partition function is given by where j now only contains all remaining indices, and X = (x, x +1) represents the sites on the coarse lattice. Note that a tensor on the coarse lattice is connected to each neighbor in the contraction direction1 by a single shared index j and to each neighbor in the perpendicular direction2 by two such indices, which form "fat indices" (in the following we therefore call the links in this direction "fat links"). In (32) we denote the fat indices of T by the pairs (j x,−2 , j x+1,−2 ) and (j x,2 , j x+1,2 ). In the following we want to apply the ideas of HOTRG to the blocked partition function (33) and reduce the bond dimension of the fat indices from D 2 back to D, the bond dimension of the original indices. 8 Note that we cannot apply the HOSVD procedure as such to the coarse numeric tensor (32) because the coarse Grassmann tensor (31) depends on the same indices through f (j).
By choosing the1-direction as the first contraction direction, we always combine T (x) on a site with odd time coordinate with T (x+1) on a site with even time coordinate, such that the contributions of the staggered phases are the same for all T on the coarse lattice. Therefore, the new numeric tensor T (X) is identical for all X on the coarse lattice. The new K (X) can again be considered to be commuting in Z, as is explained in Appendix A.

Reducing the number of Grassmann variables in the space direction
To reduce the number of Grassmann variables in the blocked Grassmann tensor, we will integrate out the Grassmann variables in the direction perpendicular to the contraction direction in (31). However, the Grassmann variables and their corresponding differentials belong to tensors K on different coarse sites in the partition function (33). The differentials belonging to the fields c x,2 and c x+1,2 in K (X) can be found in K (X+2) . Therefore, we want to reshuffle Grassmann differentials between all K in the partition function (33) to be able to integrate out the Grassmann variables in the2-direction. To do so, the differentials will be moved from the coarse site X to X −2, and will be replaced by the differentials which are moved in from site X +2 to X. This applies to all X on the coarse lattice. This reshuffling of Grassmann differentials would however introduce non-local sign factors, and the partition function would no longer have the form of a tensor network.
To resolve this problem we define new auxiliary Grassmann variablesc on the fat links of the coarse lattice by introducing a factor Note thatf X,−2 is not an independent variable, but just an alias for the expression in (37), which we call the Grassmann parity of the fat index (j x,−2 , j x+1,−2 ). This definition guarantees that the sum f x,−2 + f x+1,−2 +f X,−2 is even, such that the product is commuting. After introducing (36) in (31) and reordering the differentials and fields, we find with a sign factor The partition function (33) can now be written as with a modified numeric tensor and a new Grassmann tensor Note that the tensor K has the same six indices as K sincef is defined by (37). The integral overc X,−2 in (39) is now part of the integral in Z.
We are now able to move the commuting combination (38) from the coarse site X to X −2, for all X, without generating any sign factors. Hence, in K (X) this combination is replaced by the commuting expression which is moved in from site X +2 to X. This is done for all X on the coarse lattice. 9 The partition function (41) can now be written as with a new Grassmann tensor, where we moved the commuting combination (44) to the appropriate position to perform the Grassmann integrations over c x,2 and c x+1,2 without generating additional sign factors.
The new Grassmann tensors K (X) can always be considered to be commuting, as the entries of the corresponding numeric tensors T (X) are nonzero only when (see Appendix A)

HOSVD of the numeric tensors
In the following we will show how to apply an HOSVD approximation to reduce the dimension of the coarse-lattice numeric tensor T , by truncating its fat indices (j x,−2 , j x+1,−2 ) and (j x,2 , j x+1,2 ). As these indices also occur in K, it may seem as if this procedure cannot be applied. However, after the integration in (46), the new Grassmann tensor K (X) only depends on j x,−2 and j x+1,−2 through the sum of their Grassmann parities inf X,−2 , see (37). Similarly, it only depends on j x,2 and j x+1,2 through the sum of their Grassmann paritiesf X,2 . Therefore, truncations of T (X) are now possible if we separately truncate subspaces with even and odd Grassmann paritiesf X,2 andf X,−2 . Let us first analyze the HOSVD of the numerical tensor T . The HOSVD procedure requires the computation of the left singular vectors of the matrizations M of the coarse-lattice tensor T of (42) with respect to its fat indices. For a contraction in the1-direction, the matrization with respect to the backward2-direction yields the matrix The matrix entries of M − are just a reordering of the tensor entries of T . From (47) we see that the entries of M − are nonzero only when the Grassmann parities of its indices satisfỹ The columns of the semi-orthogonal matrix U are the D leading left singular vectors u (j X,−2 ) of M − . Therefore the D 2 × D dimensional matrix U is also block diagonal and its columns can be assigned a definite Grassmann parityg X,−2 .
Nevertheless, the matrices M − and M + generically have different singular values and singular vectors. However, for HOTRG we have to project the vector spaces of dimensions D 2 belonging to the backward and forward directions on the same D-dimensional subspace. Hence, a common D 2 × D semi-orthogonal truncation matrix U has to be constructed from M − and M + . The standard approach [8] consists of constructing U with the D leading left singular vectors of M − or M + , depending on which one yields the smallest truncation error, i.e., the largest value for the sum of their D largest singular values. However, we have developed a so-called SuperQ method [16], which reduces the combined local approximation error, defined below in (51) Therefore, in T of (42) we can truncate the fat indices (j x,−2 , j x+1,−2 ) and (j x,2 , j x+1,2 ) with dimension D 2 to new thin indicesj X,−2 andj X,2 of dimension D with Grassmann paritiesg X,−2 ≡g X,−2 (j X,−2 ) and g X,2 ≡g X,2 (j X,2 ) and construct a new tensor The error of the combined truncations of the fat indices can be quantified by the Frobenius norm Note that in (50) the same semi-orthogonal matrix U is used for backward and forward directions, such that the tensor network on the coarse lattice can be written in terms of T , with bond dimensions D on all links. 13

Applying the truncation matrices
As the indices of the numeric tensor also appear in the Grassmann tensor, we need to determine the effect of applying the truncation matrix U to the product of the two tensors. We first consider the application of U to truncate the fat index for the backward direction, Due to the block-diagonal nature of U , we observe that for thosej X,−2 which have Grassmann paritỹ g X,−2 = 0, only fat indices (j x,−2 , j x+1,−2 ) with Grassmann parityf X,−2 = 0 contribute to the sum. Similarly, forj X,−2 withg X,−2 = 1, only fat indices (j x,−2 , j x+1,−2 ) withf X,−2 = 1 result in nonzero contributions. Therefore we can replace the indexf X,−2 of K withg X,−2 , such that (52) becomes Similarly, for the truncation in the forward direction,f X,2 can be replaced byg X,2 in K. After these replacements the truncation matrix U only acts on the numeric tensor T , which leads to the truncated numeric tensor T of (50). The Grassmann paritiesg of the indicesj of T become the new coarse site Grassmann indices in the Grassmann tensor K.
As explained above, the coarse tensor T is identical for all sites X on the coarsened lattice. Therefore the truncation procedure is identical for all these sites. The links (x, −1) and (x +1, 1) on the original lattice become (X, −1) and (X, 1) on the coarse lattice, see Fig. 2, such that the new coarse local tensor T has the following entries for the coarse site X = (x, x +1), wherej are the new indices introduced in the truncation procedure (50). Using the same change of notation, the Grassmann tensor (46) on the coarse lattice becomes with f X,±1 ≡ f X,±1 (j X,±1 ) andg X,±2 ≡g X,±2 (j X,±2 ). As is explained in Appendix A, the tensors K can still be considered to be commuting. The partition function on the coarse lattice then reads where V = V /2. We can now rename and finally After this change of notation, the partition function (56) has the exact same form as the original Z in (28), albeit now on the coarsened lattice of half the volume. This means that the coarsening procedure detailed above is self-reproducing. Below we will show that the same holds for contractions in the2direction. Therefore, the blocking steps in either direction can be repeated iteratively using the exact same manipulations until the complete lattice has been reduced to a single site. Note that all sign factors generated by Grassmann manipulations only depend on local indices, i.e., indices connected to the sites being contracted. This property, which allows us to absorb the sign factor in the new numeric tensor on the coarse lattice, is crucial for the application of the iterative renormalization group procedure.
The implementation of the GHOTRG algorithm is described in Appendix B. To improve its efficiency, the truncated coarse-lattice tensor T is constructed without explicitly computing the full coarse-lattice tensors T of (32) nor T of (42).
After defining a new Grassmann tensor and a new numeric tensor the partition function (28) can be rewritten as We observe that the partition function (63) has a structure identical to (28), up to an exchange of the directions 1 ↔ 2 and a renaming of T → T and K → K. From here on, everything derived in Sec. 3.2 for a contraction in the1-direction can be applied to a contraction in the2-direction by just exchanging 1 ↔ 2 everywhere. This means that we again integrate out the Grassmann field along the contracted link, introduce new Grassmann variablesc X,1 on the fat links perpendicular to the contraction direction, move the differentials one site backward in the1-direction, and integrate out the old Grassmann variables in that direction. The coarse numeric tensor is truncated using HOSVD which yields the partition function (56) with directions1 and2 exchanged. Similarly to (57) we rename and again X → x and V → V . Then, the partition function on the coarse lattice is identical to (63), albeit on a lattice of half the volume. We now convert the Grassmann tensor back to its canonical form (26), such that further blockings in either direction can be applied, using the procedures detailed in Secs. 3.2 and 3.3. After a contraction in the2-direction, the coarse Grassmann tensor reads with sign factorσ given in (60). After defining a new Grassmann tensor which is in the canonical form, and a new numeric tensor the partition function on the coarse lattice again has its original form (28). The use of a canonical order for the variables in the Grassmann tensor conveniently allows for a flexible order of contraction directions.
The boundary conditions on the Grassmann variables are easily applied in our version of the GHOTRG procedure. As we have shown in the sections above, the boundary conditions are automatically transferred to the coarse site Grassmann variables c x,ν at each coarsening step. Therefore, the antiperiodic boundary conditions in time are implemented by imposing c x,−1 = −c x,1 in (69), while the periodic boundary conditions in space are given by c x,−2 = c x,2 . This results in and the partition function is thus given by In analogy to matrices, the sums in (71) are often referred to as tensor traces in the corresponding directions.

ASAP-tracing
As an alternative to the procedure described in Sec. 3.4, we can also apply "ASAP-tracing". As soon as the lattice has been reduced to a single one-dimensional slice, the tensor can be traced out in the perpendicular direction. This slightly improves the accuracy of the GHOTRG method since it avoids unnecessary truncations that would otherwise arise in the further coarsening steps of the remaining direction. When performing this ASAP-tracing, one first integrates out the Grassmann variables in the perpendicular direction taking into account the boundary conditions. We start from the Grassmann tensor in its canonical form, and consider a tracing in either the1-or2-direction below.

Tracing the time direction
For a slice in the2-direction, we reorder the Grassmann variables in (72) to integrate out the variables in the1-direction with f x,1 = f x,−1 and apply the antiperiodic boundary conditions in time by setting c x,−1 = −c x,1 . This leads to We now define and such that the partition function is Only the entries of the matrix T (1d) with even f x,−2 + f x,2 are nonzero, see Appendix A, such that K (1d) can be considered to be commuting. In the subsequent spatial tensor contractions, the Grassmann tensor is given by and after taking X → x and K → K the Grassmann tensor is identical to (74), albeit with x now on the coarse lattice. This means that the Grassmann tensor is self-reproducing in the one-dimensional coarsening procedure. The contraction of the two adjacent numeric tensors then yields (the Grassmann tensor no longer depends on f x,2 ), which corresponds to a matrix multiplication. When taking X → x and T → T , the partition function again looks like (76), albeit on the coarsened lattice.
After contracting the remaining sites in the2-direction until only one site is left, the final trace with periodic boundary conditions in the spatial direction yields

(80)
In this case we define the Grassmann and numeric tensors on the remaining one-dimensional lattice as and such that the partition function can be written as When blocking sites in the remaining1-direction, the product of Grassmann tensors is self-reproducing, as which after taking X → x and K → K again yields the structure of (81), albeit with x on the coarsened lattice. The contraction of the two adjacent numeric tensors then yields (the Grassmann tensor no longer depends on f x,1 ), which is a matrix multiplication. When taking X → x and T → T , the partition function again has the form of (83), albeit on the coarsened lattice.
After contracting the remaining sites in the1-direction, the final trace with antiperiodic boundary conditions in the time direction yields

Results
In this section we report about the application of our GHOTRG method for two-dimensional strongcoupling QCD with staggered quarks, where we set the anisotropy factor γ = 1.
We first consider a baryon-only version of the model to validate the Grassmann blocking without being affected by possibly large mesonic contributions. For small lattices, the numerical results are verified with exact analytic computations.
Next we report about the application of the GHOTRG method to the full strong-coupling meson-baryon system and again compare with exact results for small lattices. Furthermore, we investigate the convergence of log(Z)/V with the bond dimension D. Besides the partition function itself, we also compute the chiral condensate and the quark number density For large volumes, we study the behavior of the chiral condensate as a function of the mass and the volume at zero chemical potential, in order to investigate the chiral symmetry of the model. Finally we present results for the quark number density and the chiral condensate as a function of the chemical potential and obtain some evidence for a first-order phase transition. We implemented our version of the GHOTRG procedure as an extension to our already existing C++ HOTRG library. Some specifics of our implementation are described in Appendix B.

Computing observables with stabilized finite differences
In tensor studies, observables are often computed using finite differences of log Z or using an impurity method. To overcome the drawbacks of both methods, we developed a stabilized finite-difference (SFD) method [10]. To motivate the method, it is useful to describe the problem encountered with the traditional finite-difference computations in tensor methods. Numerical finite differences only work properly for functions that are sufficiently smooth. However, the very nature of tensor-network methods is that discrete truncations are applied during the blocking procedure, and these truncations very easily break the required smoothness property of log Z with varying parameter values. The problem occurs when the computed log Z jumps between close-by parameter values required for the evaluation of finite differences. In tensor methods, such jumps are typically caused by (almost-)degenerate singular values and/or level crossings of singular values, which lead to discontinuous changes of the vector subspaces used to truncate the coarse-lattice tensors. This problem can in principle only be resolved by taking the bond dimension D so large that the systematical error on log Z is much smaller than the difference between the exact values of log Z for two different parameter values. If such a bond dimension cannot be achieved, as is often the case, the computed finite differences will have large errors.
A solution to this problem, which we developed with the SFD method, is to modify the HOSVD truncations in order to improve the smoothness properties of the computed log Z, required for the application of the finite-difference method. The stabilization uses a heuristic approach that operates on the singular vectors of HOSVD to maximize the overlap between the truncated vector spaces constructed for adjacent parameters values. This is achieved by considering almost-degenerate singular values for both parameter values, and introducing separate basis changes in the respective subspaces. The method uses the fact that small variations of the parameter values generically lead to small rotations of these subspaces and allows for the use of very small step sizes in the finite-difference formula.
Note that observables can also be computed using the impurity method. Although this method yields smoother data (which does not necessarily mean more accurate) than the non-stabilized finite-difference method, it has its own systematic error because the same singular vectors are used to truncate the pure and impure tensors. We therefore use the SFD method as method of choice to compute observables. The SFD method was also used successfully to stabilize second-order finite differences in the computation of susceptibilities, e.g., the specific heat of the three-dimensional O(2) model [10].

Baryon-only partition function
As a first validation of our GHOTRG method for strong-coupling QCD, we discard all mesonic contributions in (28), i.e., we replace δ x∈M → 0 in (14). Then, the resulting partition function is independent of the mass and all sites of contributing configurations are baryonic.
We computed log(Z)/V as a function of the chemical potential µ on lattices of sizes 2 × 2, 2 × 4, 4 × 2, 4×4, 4×8 and 8×4 using GHOTRG with D = 32 and compared these results with the analytical predictions given in Appendix C. We find very good agreement between the numerical and analytical results. This is illustrated in Fig. 3 for the 2 × 2 and 4 × 4 cases. These results confirm that the global minus signs that appear in the standard baryon-loop formulation of the model [2] are correctly taken into account by the GHOTRG procedure.

Meson-baryon partition function
Next we consider the full meson-baryon system of strong-coupling QCD in two dimensions, described by the tensor network (28). In Fig. 4 we show log(Z)/V as a function of the chemical potential for a 2 × To verify the accuracy of the GHOTRG results, we show their relative deviation from the exact result on a 4 × 4 lattice for m = 0 and m = 0.2 as a function of µ for various bond dimensions D in Fig. 5. As expected, the accuracy typically improves with increasing D, but the behavior does not hold for all µ and D, which is related to the small size of the lattice, see below. For m = 0 (left plot) the error is about a factor of 10 larger than for the nonzero mass m = 0.2 (right plot). This shows that the tensor method is more accurate for larger masses, as is the case with most other simulation methods. Nevertheless, even in the chiral limit (m = 0), the tensor method gives very satisfying results for this two-dimensional system. We also compute the mass dependence of log(Z)/V at zero chemical potential and verify our results with the analytic expression (C.10) on a 4 × 4 lattice, see Fig. 6. From the relative deviation , shown in the bottom row of the figure, we observe that the accuracy improves as D becomes larger and the results converge to the exact values. However, larger values of D are required to get accurate results for smaller masses.
In Fig. 7 we show a convergence study of log(Z)/V with respect to the bond dimension D for m = µ = 0 on 4 × 4 and 1024 × 1024 lattices. The convergence behavior is quite erratic on the small lattice. Even though the accuracy is very good when D is sufficiently large, the convergence is far from being monotonous. For the large lattice, the convergence is much more stable, and a quadratic fit in 1/D allows us to make an extrapolation to D → ∞.

Chiral condensate at zero chemical potential
After validating the GHOTRG method for small lattices, where analytical results are available, we now consider larger lattices of size L × L. First we compute the chiral condensate (87) as a function of the mass and lattice volume at zero chemical potential. The aim is to investigate if the chiral symmetry is dynamically broken in this two-dimensional theory. To this end we look at the zero-mass and infinite-volume limit of the chiral condensate where the order of the limits is crucial. In Fig. 8 we show the evolution of the chiral condensate as a function of the mass for various lattice sizes at fixed D = 64. Although these results are computed at fixed D, they already illustrate how the chiral condensate converges to its infinite volume limit for the different mass values. As the mass gets smaller, larger volumes are needed to approach this limit.
We now perform a detailed analysis of the chiral condensate by extrapolating to D → ∞ for each mass and volume, and then extrapolating this infinite-D result to V → ∞ for each mass value. An example for such an extrapolation, together with its error estimate, is shown in Fig. 9 for m = 10 −5 . We observe that the lattice size needed to obtain an estimate for the V → ∞ limit increases with decreasing mass. The results of these extrapolations as a function of the mass are shown in Fig. 10. For small m < 0.005, the results lie on a straight line in a log-log plot and are thus well fitted by lim V →∞ ψ ψ = am b . For the fit shown in Fig. 10 (left), the fit parameters are given by a = 2.77, b = 0.0414. This shows that the chiral symmetry is not dynamically broken in two-dimensional strong-coupling QCD with (two tastes of) staggered quarks.
For large masses, the chiral condensate is asymptotically given by 3/m at leading order, which can easily be derived from the partition function (6). We therefore fit the infinite-volume limit of the chiral condensate over the full mass range by the empirical formula f (m) = (am b +cm)/(1+dm+(c/3)m 2 ), which interpolates between the asymptotic behaviors, see Fig. 10 (right). The fitted parameter values are a = 2.77, b = 0.0409, c = 1.05, d = 0.770.

Particle number density and chiral condensate at nonzero chemical potential
Finally, we use the GHOTRG method to investigate the behavior of the model at nonzero chemical potential. For m = 0.1 we study the quark number density (88) and the chiral condensate (87) as a function of the chemical potential, for lattice sizes up to L = 128. The GHOTRG results for D = 64 are shown in Fig. 11. For small lattices (2 × 2 and 4 × 4), they agree well with the exact values obtained from the analytic formulas of Appendix C. For larger lattices, the results quickly converge to the V → ∞ limit. There appears to be a first-order phase transition around µ c ≈ 0.3508. Above this critical chemical potential, we observe both a nonzero quark number density ρ and a restoration of the chiral symmetry, i.e., ψ ψ = 0.

Conclusions
In this paper we developed a tensor-network renormalization group framework, based on the GHOTRG method, specifically tailored for strongly-coupled two-dimensional QCD with staggered quarks. In its dual formulation, the partition function is composed of mesonic and baryonic degrees of freedom. The Grassmann variables in the baryonic contributions to the partition function cannot be integrated out without introducing    non-local sign factors. Therefore, the partition function cannot be written as a full contraction of a numeric tensor network and the standard HOTRG method cannot be applied. However, this problem can be resolved by constructing a tensor network consisting of local numeric and Grassmann tensors. When the lattice is then coarsened during the renormalization group procedure, the blocking of two adjacent sites now consists of two steps: First the Grassmann tensors on the two adjacent sites are blocked, yielding a new Grassmann tensor on the coarse lattice. This procedure generates a local sign factor, which is absorbed in the contraction of the numeric tensors on the two adjacent sites. Just as in standard HOTRG, the coarse numeric tensor is then subjected to an HOSVD approximation to avoid an exponential increase of its dimensionality. After each renormalization group step, the partition function is represented by a coarsened tensor network of local tensors that are again products of a numeric and a Grassmann tensor. At each blocking step, the number of Grassmann variables is reduced by a factor of two and the HOSVD procedure reduces the dimensions of the fat indices of the numeric tensor back from D 2 to D. This procedure is repeated until the whole lattice has been reduced to a single site and the partition function can be computed, taking into account the boundary conditions. Our version of the GHOTRG procedure allows for a tensor-network computation of the partition function with a computational cost that is similar to that of standard HOTRG. This can be achieved since the Gram matrices, used in the construction of the truncation matrices, are block diagonal in the Grassmann parity. Without this block-diagonal structure, truncations from dimension 4D 2 to D would be required at each blocking step to keep the dimensions of the coarse local tensors under control, as is the case in applications of the GHOTRG method for some other fermionic problems [13]. This would make the GHOTRG method substantially more expensive.
We have validated our version of the GHOTRG method by comparing with exact results on small lattices. On large lattices, we have studied the chiral condensate as a function of the mass and the volume at zero chemical potential, and showed that the chiral symmetry is not dynamically broken in this two-dimensional model in the chiral limit (m → 0). At nonzero chemical potential we computed both the quark number density and the chiral condensate and found some evidence for a first-order phase transition to a phase with nonzero density where the chiral symmetry is restored (for nonzero mass).
In future work we will apply the method to strong-coupling QCD in higher dimensions and also extend it beyond the infinite-coupling limit.
For a contraction in the1-direction, the matrization with respect to the backward2-direction leads to the Gram matrix Note that the Grassmann parities f in the sign factors σ are functions of the corresponding indices j.
To improve the efficiency of the computation of Q − and to reduce the required storage, we would like to reshuffle the factors in the previous expression, such that the tensors T at the same positions are contracted first, as is usually done in standard HOTRG. The additional couplings between the tensors, caused by the sign factors, complicate the reordering of the product.
As the nonzero entries of M − satisfy (49), the nonzero entries of Q − have Grassmann parities satisfying such that Q − is block diagonal with the nonzero blocks being either even-even or odd-odd blocks in f We now reorder the sums in (B.1) such that the tensors on equal sites can be contracted first.
The sums over the indices that only appear in the two factors of T (x) yields where the indices j (0) and j (1) in the sum denote the indices j with even and odd Grassmann parities, respectively. The storage and computational costs scale as Mem ∝ 2D 4 and Comp ∝ D 6 , respectively. Note that we cannot add the f x,2 = 0 and f x,2 = 1 contributions when constructing the auxiliary tensor A − since f x,2 also appears in σ fx,2f x+1,−2 f x+1,−2 . Analogously, for the two factors T (x+1) , we construct with Mem ∝ D 4 and Comp ∝ D 6 . Finally we compute with Mem ∝ 2D 4 and Comp ∝ 2D 6 , by first contracting A − and B − and then summing over f x,2 .

Appendix B.3. Staggered phase
To avoid multiple definitions of the initial local tensor that would just differ in the staggered phase, we decided not to include the latter in the tensor (14), but instead to modify the construction of U and the truncation T of the coarse tensor, described in the sections above, to take into account the staggered phase. If we choose the very first contraction to be in the1-direction, these modifications only have to be applied in this contraction. In subsequent contractions, all local tensors on the coarse lattice are identical, and the constructions of U and T are performed as described above, without any modifications. Therefore, to take into account the staggered phase, we only need to modify the equations for the first contraction, by replacing T (x) by η fx,2 x,2 T (x) and T (x+1) by η f x+1,2 x+1,2 T (x+1) , see also (14). From the definition of the staggered phase, these two staggered phases will have opposite signs. We choose η x,2 = 1 and η x+1,2 = −1, such that only operations involving T (x+1) will be affected.

Appendix C. Analytical results
We developed a code to generate all configurations with baryon and meson loops automatically, from which an exact analytic formula for lattices up to 8 × 4 can be computed. The code can be restricted to baryon-loop only configurations or applied to the full meson-baryon system, with chemical potential µ, mass m, anisotropy parameter γ, and arbitrary combinations of periodic and antiperiodic boundary conditions. The results presented below use periodic boundary conditions in space and antiperiodic boundary conditions in time. Although we derived the analytical results for arbitrary γ, we only give the formulas for γ = 1 for conciseness.