Tensor Network study of the (1+1)-dimensional Thirring Model

Tensor Network methods have been established as a powerful technique for simulating low dimensional strongly-correlated systems for over two decades. Employing the formalism of Matrix Product States, we investigate the phase diagram of the massive Thirring model. We also show the possibility of studying soliton dynamics and topological phase transition via the Thirring model.


Introduction
Tensor Network (TN) methods [1][2][3] provide a non-perturbative tool for studying (lattice) quantum many-body systems. A common ingredient of TN algorithms is an entanglement-based ansatz for the quantum many-body wave function, which allows an efficient description of the system in spite of the exponential growth of the Hilbert space dimension with the number of lattice sites. Tensor network methods have been extremely successful in one-dimensional problems, and are becoming competitive for higher dimensional ones. They are free of the sign problem, and can be used to approximate ground states, low-lying excitations, thermal states, and (to some extent) real-time evolution. In the last years they have been applied to the Hamiltonian formulation of lattice gauge theories .
In this article, we consider the zero-charge sector of the massive Thirring model, where the model is dual to the sine-Gordon theory [27]. In fact, early explorations of the Thirring model suggested that the massless case is equivalent to a free scalar theory in two dimensions [28][29][30]. Without the presence of mass, a fermion-antifermion pair created from the vacuum does not separate and forms a bound state in two dimensions. Although, in general, this would not be possible if the fermions have mass, Coleman illustrated that we can still have a bosonic description for the zero-charge sector of the massive Thirring model [27]. Since in Euclidean space the sine-Gordon theory contains topological structures, this connection allows one to simulate soliton dynamics through investigating the massive Thirring model.
On the other hand, the sine-Gordon theory in Euclidean space is also dual to the vortex sector of the classical XY model. This model undergoes a phase transition without the appearance of long-range order in two dimensions, which is known as the Berezinskii-Kosterlitz-Thouless (BKT) transition [31, Speaker, e-mail: tanlin2013.py04g@nctu.edu.tw 32]. It has been found that the topological excitations, vortices and anti-vortices, play an important role in this system. At low temperatures, a vortex forms a bound state with an anti-vortex. The vortexanti-vortex pairs are condensed in the quasi-long-range order phase, where the correlation functions decay with the distance as a power law. At high temperatures, the system is in a disordered phase and the correlators exhibit the behavior of exponential decay. The BKT phase transition happens without breaking the O(2)-symmetry of the XY model. Instead, these two phases are distinguished by different scaling laws of the correlators.
This article is organized as follows. We begin by presenting the definition of the models and their relationships in Sec. 2. We explain our implementation for lattice simulations in Sec. 3. In Sec. 4, we formulate the Matrix Product State (MPS) ansatz and the algorithm we adopt. Finally, we present our preliminary numerical results in Sec. 5, and conclude in Sec. 6.

Models and Dualities
The massive Thirring model is a theory of a single Dirac field on a (1+1)-dimensional space-time, with the action defined by where m 0 denotes the bare mass, and g is a dimensionless coupling constant.
In [27], the equivalence between the zero-charge sector of the massive Thirring model and the sine-Gordon model was established. The sine-Gordon model is a theory of a single scalar field on a (1+1)-dimensional space-time, with the action defined by where α 0 is a dimension-two bare coupling, and t is dimensionless. The correspondence between operators and couplings is summarized bȳ where Λ is the ultra-violet cutoff. In particular, the first formula of Eq. It is well known that the sine-Gordon model is also dual to the vortex sector of the classical XY model described by the Hamiltonian where θ is an angle defined on the two-dimensional lattice, and K is a coupling constant. Below the critical temperature, T c (T < T c ), all angles are almost aligned. Yet, the long-range order is destroyed by thermal fluctuations. Despite this, the system still has a quasi-long-range order phase, with the correlation function decaying algebraically. At high temperatures (T > T c ), the system is in a disordered phase and the correlator decays exponentially. To summarize, where ξ is the correlation length. In addition, we can describe the vortex sector of the XY model, which is equivalent to the Coulomb gas, by the grand canonical partition function [33], where β T = 1 T and G L is the Green's function for the two-dimensional Laplacian operator. We also include the vortex density n r and the core energy of each vortex c . One can find the sine-Gordon representation of Eq. (6) as, where β t = 1 t , and the sine Gordon action Eq. (2) is now written in Euclidean space (and thus the sign in front of α 0 is flipped). The correspondence of parameters is identified by where A is a constant, and e −β T c is known as the fugacity. Finally, we summarize the correspondence of the parameters of these three models in Table 1.

Preliminaries of the Lattice Calculation
In our approach, the system is described by the quantum many-body wave function. To start with, we employ staggered fermions as proposed by Kogut and Susskind [34]. In this formulation the lattice Hamiltonian derived from the action of Eq. (1) reads where c † n and c n are fermionic creation and annihilation operators obeying the anti-commutation relations, {c † n , c m } = δ nm , {c n , c m } = 0, {c † n , c † m } = 0, and a is the lattice spacing. Furthermore, the fermion variables can be mapped onto spins by using the Jordan-Wigner transformation,  with S + n and S − n being the ladder operators of spin-1/2. We note that this step is not strictly required, since TN can be applied to fermionic degrees of freedom. However, it is convenient to construct the full Hamiltonian using the Pauli matrices. Therefore, the spin Hamiltonian we use in the simulation is given by We will explain how to solve the Hamiltonian numerically in Sec. 5. In addition, the total S z polarization S z tot = n S z n is a conserved quantity, related to the charge in the fermionic language. We can thus explore the zero-charge sector by requiring S z tot = 0. Further details will be illustrated in the next section.

Tensor Network Methods
The concept of entanglement is central to the Tensor Network states and the algorithms based on them. Any pure state |Ψ of a composite Hilbert space where s i 's, known as Schmidt coefficients, are positive, monotonically decreasing real numbers which determine the entanglement with respect to the bipartition A|B, and r ≤ min(d A , d B ). We can compute the entropy of entanglement where ρ A is the reduced density matrix for the subsystem A. If the system is weakly entangled, the Schmidt coefficients often decrease very fast. This provides us with a way to approximate the state by introducing an entanglement cutoff which truncates the number of the Schmidt coefficients from r to r < r.

Matrix Product States
For a N-partite system, with Hilbert space H = H 0 ⊗ H 1 ⊗ · · · ⊗ H N−1 , H n = C d for every n, a MPS is a state of the form [3] |Ψ = a 0 ,a 1 ,··· where M[n] σ n are D × D matrices for 0 < n < N − 1, and vectors of dimension D for n = 0, N − 1 (see Fig. 1(a)). The parameter D is called the bond dimension. Any state of the Hilbert space can be written as a MPS for D ≤ d N/2 , while restricting D to a small value plays the role of a cutoff in the Schmidt rank.

Density Matrix Renormalization Group
Our objective is to find a MPS approximation to the ground state (and excitations) of the model. This can be achieved by a variational minimization of the energy over the set of MPS with fixed bond dimension D, similar to the Density Matrix Renormalization Group (DMRG) algorithm [3]. The algorithm is simplified when the Hamiltonian is written as a Matrix Product Operator (MPO), i.e. a MPS in the operator vector space [1,2], see Fig. 1(b).
where W[n] are χ × χ matrices (χ dimensional vectors for the edges, since we adopt open boundary conditions). We start with expressing the functional in terms of a Tensor Network as Fig. 1(c). If the M[n] tensors for all but one site, k, are fixed, Eq. (16) can be optimized over M[k] by solving a generalized eigenvalue equation [3]. This can be repeated for all sites, from left to right and back, until sufficient convergence of E/N is achieved. The Hamiltonian in Eq. (11) conserves the total S z polarization S z tot = n S z n . It is thus convenient to restrict the search to a sector of specific S z tot = S target . To this end we add a penalty term [5] to the Hamiltonian where the magnitude of λ should be chosen to be large enough to ensure that all the states with Ψ| N−1 n=0 S z n |Ψ S target have energy above the lowest state in the desired sector. The MPO representation of the Hamiltonian of Eq. (17) is therefore given by for the tensors at the boundaries and for the remaining sites, where

Results
In this section, we present our results obtained in the presence of the penalty term with S target = 0. Figure 2 illustrates our investigation of the ground state (GS) energy, extracted using the DMRG algorithm as described in Sec. 4, for system size N = 40 and bond dimension D = 100. The plot in Fig. 2(a) shows a scan on the (g, m 0 L) plane, with L denoting the spatial volume of the lattice. In the regime g < 0, it is evident that the GS energy decreases when g turns more negative. This is in accordance with Coleman's argument that the GS energy is always negative and scales as λ 1/(1+g/π) , where λ is a dilatation parameter [27]. Figure 2(b) exhibits the GS energy per site (E/N) as a function of g. We observe that when g −1.5, E/N is almost independent of the mass parameter m 0 . To understand this behavior, we notice that the renormalization group (RG) analysis of the sine-Gordon theory concludes that the (cosφ) operator is irrelevant when t > 8π. This corresponds to the region g < −π/2. In other words, the Thirring model is indeed describing a free bosonic theory in this regime. This can lead to the observed independence of m 0 .
To further examine the convergence of the algorithm, we also check S z tot and the bipartite entanglement entropy computed using the GS extracted from our DMRG simulations. To calculate the bipartite entanglement entropy, we divide the system into two equal-size subsystems, A and B, and then determine the entanglement entropy S A|B using Eq. (13). Figure 3 displays results from this study. We notice that in Fig. 3(a) the GS does not always converge to the desired S z tot sector, even in the presence of the penalty term with S target = 0. This happens when g −1. In addition, it is observed in Fig. 3(b) that S A|B shows unstable, scattering behavior. In order to understand this instability, we perform DMRG simulations starting from a GS with S z tot = 0 at g < −1, and slowly changing g. Namely, we use the ground state with S z tot = 0 obtained at a value of g as our initial MPS  in the DMRG computation for another choice of g. With this procedure, we find that the instability demonstrated in Fig. 3 disappears. Presently we are carrying out further investigation of this issue.

Conclusion and Outlook
In this project, we implemented the MPS method for the massive Thirring model on the lattice, employing staggered fermions. This article presents our exploratory numerical results from this work. Our initial findings show that the DMRG computation for this model can be performed. Using the dualities described in Sec. 2, we will investigate the topological phase transition and the soliton dynamics in the XY model and in the sine-Gordon theory.
It is well known that there is a BKT-type phase transition at T ∼ πK/2 in the XY model. From Table 1, it means that this transition can be observed at g ∼ −π/2 in the massive Thirring model. The issue is to use the dualities and the Jordan-Wigner transformation to write the XY-model vortex-antivortex correlators in the format of MPO and MPS. This aspect of the project is now under investigation, and will appear in a separate publication in the near future [35].