Phase structure of the 1+1 dimensional massive Thirring model from matrix product states

Employing matrix product states as an ansatz, we study the non-thermal phase structure of the (1+1)-dimensional massive Thirring model in the sector of vanishing total fermion number with staggered regularization. In this paper, details of the implementation for this project are described. To depict the phase diagram of the model, we examine the entanglement entropy, the fermion bilinear condensate and two types of correlation functions. Our investigation shows the existence of two phases, with one of them being critical and the other gapped. An interesting feature of the phase structure is that the theory with non-zero fermion mass can be conformal. We also find clear numerical evidence that these phases are separated by a transition of the Berezinskii-Kosterlitz-Thouless type. Results presented in this paper establish the possibility of using the matrix product states for probing this type of phase transition in quantum field theories. They can provide information for further exploration of scaling behaviour, and serve as an important ingredient for controlling the continuum extrapolation of the model.


I. INTRODUCTION
Many quantum field theories of interest cannot be studied with perturbative methods. This concerns, most notably, quantum chromodynamics (QCD), but also seemingly simple models may have non-perturbative features. Since the seminal idea of Wilson [1], the method of choice for such systems is usually to formulate them on the lattice. In many cases, e.g. in QCD, this is the only way to obtain quantitative predictions directly from the Lagrangian. The path integral corresponding to a discretized system is finite-dimensional, but usually this dimension is very large, implying no alternative to Monte Carlo sampling. If such sampling is possible, the lattice provides an unambiguous way of reaching results with an arbitrary precision, with controllable total error. This led to spectacular successes, particularly in the most important theory studied with lattice methods, i.e. QCD. Several aspects of QCD were addressed with large-scale simulations, including hadron spectroscopy, hadron structure, thermal properties or fundamental parameters, such as the strong coupling constant and its scale dependence or the Cabibbo-Kobayashi-Maskawa matrix parameters. However, there are still areas in lattice field theories where the lattice has not offered quantitative or even qualitative insight due to an inherent shortcoming of the Monte Carlo approach. This shortcoming is that not all theories admit a positive-definite probability measure, a necessary ingredient for stochastic sampling. In the absence of positive-definite measure, a theory is said to have a sign problem, with a notable example of QCD at a non-zero chemical potential. Even though techniques to alleviate this issue have been invented, no real solution readily available for QCD exists. Another important example is real-time simulation, notoriously hard in quantum mechanics and quantum field theory in general.
Given these challenges, alternative approaches are being pursued. Methods tested in the context of lattice gauge theories include e.g. Lefschetz thimbles [2,3], complex Langevin simulations [4,5] and density of states methods [6]. In this paper, we concentrate on the Hamiltonian approach in the framework of tensor networks (TN). Tensor-network methods stem from quantum information theory and they have been heavily applied for description of, in particular, low-dimensional condensed matter systems. Their main feature is that they can represent quantum states, taking correctly into account the relevant entanglement properties of a system. In this way, they can define effective ansatzes for the wave function, which combined with efficient algorithms to utilize such ansatzes, can lead to a successful quasiexact description of a broad class of physical theories. A particularly efficient example of TN for one-dimensional systems is the so-called matrix product states (MPS) ansatz [7][8][9][10][11][12]. Generalizations to higher dimensions also exist, such as Projected Entangled Pair States (PEPS) [13]. For a pedagogical introduction to the TN approach, see e.g. Refs. [14][15][16][17].
The TN framework, mostly in combination with the Hamiltonian formulation, has also been used for solving lattice field theory models, such as the U (1) (the Schwinger model) and SU (2) gauge theories, the φ 4 -theory and the O(3) σ model. Among investigations that have been performed, one can name spectral computations [18][19][20][21][22][23][24][25], calculations of thermal properties [26][27][28][29][30], phase diagrams [31][32][33][34][35][36], entanglement properties [25,[37][38][39][40], non-zero chemical potential [41][42][43], studies of the θ-vacuum [44,45], and real-time evolution [19,37,38,[46][47][48]. The latter two are examples of computations hindered by a sign problem in traditional Monte Carlo simulations, which is, by construction, absent in the TN approach. From a more formal perspective, the interplay of symmetries with TN has been a fertile research topic [49][50][51] that has allowed, among others, a full classification of one dimensional gapped phases of matter [52][53][54] and the construction of topological states in two dimensional systems [55][56][57]. In the particular case of gauge symmetries, it is possible to explicitly incorporate the local gauge invariance into the ansatz, which can then be used to construct invariant models or perform numerical calculations [23,[58][59][60][61][62]. There are several aspects to these studies. From the computational perspective that occupies us in this work, their achievements are threefold. Firstly, it is important to check whether physics of quantum field theoretical models can be well described in the TN language 1 . Investigations by several groups led to unambiguously positive conclusions and the precision reached in these calculations in (1+1)-dimensional cases is in many cases unprecedentedly high. Secondly, TN methods have offered some new insight about well-known models, in particular about their entanglement properties, behaviour at non-zero particle density and real-time evolution. Finally, there is a close connection between TN formulations and proposals of quantum simulation for studying these theories in experiments with cold atoms, trapped ions or superconducting qubits [73][74][75][76][77][78]. Another related approach is the exploration of field theories employing full fledged quantum computers (for recent progress, see, e.g. Refs. [79][80][81]). Indeed, the connections among (discretized formulations of) quantum field theories, condensed matter models and quantum information techniques are promoting a fruitful research ground, where the interdisciplinary effort can shed light on various research areas [82][83][84].
In this work, we concentrate on the investigation of the (1+1)-dimensional massive Thirring model 2 , with the Minkowski-space action where m denotes the fermion mass, and g is the dimensionless four-fermion coupling constant. The current paper summarizes our result for exploring the non-thermal (zero-temperature) phase structure of the above theory using the staggered-fermion regularization [86,87]. This is the first step in a research programme on studying the Thirring model with TN methods. As described below in this section, the phase structure of the massive Thirring model can exhibit features such as infrared (IR) slavery, and the existence of a critical phase with a Berezinskii-Kosterlitz-Thouless (BKT) transition [88]. Understanding these features in lattice field theory using the TN approach is a worthy challenge. In particular, the connection between the numerical results and the continuum limit is a non-trivial aspect of the study, which nevertheless is of fundamental importance for further applications of TN methods to other quantum field theories. In addition, such exploration is essential for our future long-term work. Our final goal in this research programme includes the investigation of the model with chemical potential, and various aspects of its real-time dynamics [89].
The sector of vanishing total fermion number in the (1+1)-dimensional massive Thirring model is known to be S-dual to the sine-Gordon scalar field theory with zero total topological charge [90]. The sine-Gordon theory is described by the classical action, with two couplings, t and α. In Ref. [90], Coleman obtained the duality relations where Λ is a cut-off scale. In addition to the soliton solutions that can be related to fermions in the Thirring model [91], the sine-Gordon theory exhibits interesting scaling behaviour and phase structure, which can be understood from studying its renormalization group (RG) equations [92][93][94]. Since this aspect of the theory is important for the current work, here we describe the scenario in slightly more detail.
For convenience, let us define the following parameters, In terms oft and z, the RG equations (RGE's) of the sine-Gordon theory to O(z 3 ) are [93,94] where µ is the renormalization scale. From these RGE's, the crucial feature that can be seen immediately, in the limit z/Λ 2 1, is the following scaling behaviour. 2 The massless model was proposed by W. Thirring in 1958 [85] as an example of a solvable quantum field theory.
• In the regime wheret > ∼ 1/8π, i.e., t < ∼ 8π, the operator [cos φ(x)] in the sine-Gordon theory is relevant, resulting in the existence of solitons in the model.
• On the contrary, att < ∼ 1/8π, i.e., t > ∼ 8π, the operator [cos φ(x)] is irrelevant, and the model is a free bosonic theory at low energy. In this case, from Eq. (5), the coupling t will be scale-invariant in the IR.
These different scaling properties imply the existence of a phase transition at t ∼ 8π. To further understand the nature of this phase transition, we note that the sine-Gordon theory is known to be equivalent to the Coulomb-gas (vortex) sector of the two-dimensional O(2) σ model (the XY model), and the couplingt can be interpreted as the temperature in the latter 3 . Therefore, the above scaling behaviour can be shown to be closely related to the BKT phase transition [88]. In fact, it can be easily verified that RGE's for the sine-Gordon model, Eqs. (5) and (6), bear the same characteristics as the Kosterlitz scaling equations [96] for the XY model. This is also in accordance with the Mermin-Wagner-Coleman no-go theorem [97,98]. In a recent paper [99], the MPS approach has been implemented for the XY model to study the BKT transition in statistical physics. As we will demonstrate in this work, performing computations with the dual Thirring model using the MPS method enables us to probe this interesting phase structure. In addition, it also allows for future investigation of real-time dynamics of this BKT phase transition [89].
Employing the duality relations in Eq. (3), the sine-Gordon RGE's, Eqs. (5) and (6), can be rendered into with higher-order corrections in (m/Λ) n . These RGE's govern the scaling behaviour of the massive Thirring model in the limit m/Λ 1. Equation (7) implies that g always increases towards the IR limit as long as m/Λ = 0, although it does not signify asymptotic freedom. It also concurs with the fact that the massless Thirring model is a conformal field theory [85]. To understand the implication of Eq. (8), we first notice that, for the duality relation in Eq. (3) to be valid, the four-fermion coupling is restricted to be in the range Since we rely on the duality to understand the RG flow of the model and thus how to correctly approach the continuum, we need to restrict our exploration to this range of values. With this constraint, the RGE in Eq. (8) leads to the expectation that a corresponding non-thermal phase transition occurs in the massive Thirring model at in the regime where m/Λ 1, with the exact value of g * being (m/Λ)−dependent 4 . While the fermion mass, m, is a dimension-one coupling and remains relevant in the region g > g * , it becomes irrelevant at g < g * . Obviously, equations (7) and (8) This means that the "fixed line" on the g−m plane, m = 0 where both β g and β m vanish, is separated into two sectors, with g <ḡ * being stable and g >ḡ * being unstable. It is also straightforward to show that generically, from Eq. (8), g * (m/Λ) decreases with growing m/Λ under the condition m/Λ 1. Figure 1 shows qualitatively the above features for the RG flows of the Thirring model in the regime where −π < g and m/Λ < ∼ 0.01. Furthermore, Eqs. (7) and (8) predict that theψψ operator becomes irrelevant in the regime g < g * (m/Λ), although its classical dimension is one. This indicates that a large anomalous dimension is generated in the Thirring model. As will be discussed in Sec. V B, these features can guide the study of the continuum limit of the model.  7) and (8) in the regime where −π < g and m/Λ < ∼ 0.01. The arrows present the flows towards the IR limit. The line m = 0 is a fixed line under RG transformation. It is separated into two sectors, with g <ḡ * being stable and g >ḡ * being unstable.
It should be stressed that the above discussion is based on an expansion in terms of m/Λ. In this project, we carry out non-perturbative study for the non-thermal phase structure of the massive Thirring model through lattice simulations, employing the method of MPS. Our investigation can shed light on the scaling behaviour of the theory beyond perturbation theory. As mentioned earlier in this section, results reported in this paper are the first step in a research programme where we will further explore other aspects of this model, such as its properties out of equilibrium and relevant real-time dynamics. Although the massive Thirring model has been extensively examined in the past, our study presented in this article does reproduce known properties of the theory with novel ingredients. In addition to implementing a new strategy for research on lattice-regularized Thirring model, to our knowledge, the discussion in Sec. IV A is the first time that entanglement entropy is used for probing the zero-temperature phase structure of this quantum field theory.
The rest of this paper is organized in the following way. Section II contains the formalism of the theory that we simulate, and Sec. III gives details of the numerical implementation. In Sec. IV, we present main numerical computations in this project. The outcome of these computations is used in Sec. V for addressing the phase structure, the scaling behaviour, as well as the continuum limit of the massive Thirring model in 1+1 dimensions. We then conclude in Sec. VI. Preliminary results of this work were presented in our contributions to the proceedings for the Lattice conferences in Refs. [39,40].

II. LATTICE FORMULATION AND THE CORRESPONDING SPIN MODEL
In this section, we first describe subtleties in constructing the Hamiltonian at the operator level in the continuum, then discuss the latticization procedure of the system and the comparison with the quantum spin-chain model. In our numerical implementation, we use the XXZ-model Hamiltonian in Sec. II B.

A. The Hamiltonian operator at quantum level and the staggered regularization
To perform lattice simulations using the MPS approach, we first have to obtain the corresponding Hamiltonian of the classical action in Eq. (1). At the quantum level, the Hamiltonian operator cannot be related to the Lagrangian through a straightforward Legendre transform. The main subtlety arises from quantum effects that modify the current conservation laws, leading to an ambiguity in defining the vector current that appears in the four-fermion operator in Eq. (1) 5 . In the path integral formalism, these effects are easily understood via analysing anomalies that result from the fermionic measure in a field-redefinition procedure [100]. When working with the operator formalism, this ambiguity can be accounted for by employing a non-local definition of the currents [101,102]. As explained in the Introduction, we are interested in studying the sector of zero total fermion number in the Thirring model. To ensure that fermion number is conserved, we choose to maintain vector-current conservation in this work. Incorporating quantum effects in the energy-momentum tensor at the operator level, one derives the Hamiltonian [102], where Z ψ (g) is the wavefunction renormalization constant that can be found in Refs. [103,104]. Notice that the structure of the four-fermion operators and their couplings in H Th are different from that in the action of Eq. (1). This originates from the effects of quantum anomalies discussed above. Furthermore, in Eq. (12) we denote the fermion mass as m 0 , to emphasize that in the simulations the input mass is bare. As discussed below in Sec. II B, the input g can be related to the spectrum of a one-dimensional spin chain [105][106][107].
In this project, we work with the standard representation of the Dirac matrices in two dimensions, where σ i 's are the Pauli matrices. In this representation, the Hamiltonian can be expressed in terms of the two components of the Dirac spinor, ψ 1 and ψ 2 , whereg Having obtained the Hamiltonian operator in Eq. (14), we now proceed to discretize it, following the strategy as detailed in Refs. [86,87]. Upon latticizing the one-dimensional space, we implement the staggered regularization by keeping only ψ 1 on the even site, and ψ 2 on the odd site. That is, where a is the lattice spacing, and c i is the resulting dimensionless one-component fermion field on the i-th spatial lattice site. This procedure leads to the discretized theory with with N being the total number of sites. Notice that since we are only latticizing the spatial direction in 1+1 dimensions, the staggered regularization completely removes the fermion doubling problem. That is, H (latt) Th describes only one "taste" of fermion with the effective lattice spacing being 2a.

B. The XXZ spin chain
The Hamiltonian in Eq. (17) contains fermionic degrees of freedom, which TN can in general handle without causing a loss in efficiency [108][109][110]. But in the case of one spatial dimension, it is more convenient to map the fermions to spin matrices by applying a Jordan-Wigner transformation, where S ± n = S x n ±iS y n , and S α n = σ α /2 (α = x, y, z) with [S α n , S β m ] = 0 when n = m. Application of this transformation on H (latt) Th results in the quantum spin chain, In this work, we adopt open boundary conditions. This avoids the appearance of boundary operators that make the Jordan-Wigner transformation more complicated than Eq. (18) and, most importantly for our simulations, it allows us to use the most efficient variants of the MPS algorithms 6 .
The spin-chain Hamiltonian in Eq. (19) is precisely that of the XXZ model coupled to staggered and uniform external fields, in a regime in which the latter is equal to the anisotropy. But using exactly the form in Eq. (19) turns out not to be the adequate strategy, as the choice of parameters in terms of g does not account for the lattice effects, an issue that has been studied in the literature of condensed matter physics [105,112]. Based on a matching of the soliton bond-state spectrum between the exactly solvable field theory at m 0 = 0 and the Bethe-ansatz solution of the XXZ chain, the correct choice was found to be with where the parameter γ can be related to the four-fermion coupling, g, in the Thirring model by [105,107], Notice that the functions∆(g) andg(g) (in Eq. (15)) agree to O(g 2 ). The parameter ν(γ) accounts for the wavefunction renormalization.
For convenience, we define wherē This Hamiltonian,H sim , is what we actually implement for the simulations in this work. That is, the exploration of the phase structure of the (1+1)-dimensional massive Thirring model is performed via scanning the two spin-chain parameters, am 0 and ∆(g) in Eq. (24), that can be translated back to the two couplings in Eq. (12). It is obvious from Eq. (25) that we only explore the regime −1 < ∆(g) ≤ 1. This corresponds to −π < g ≤ π, which encompasses the expected phase-transition point, g ∼ −π/2, discussed in Sec. I, and keeps ∆(g) single-valued.
The XXZ Hamiltonian with a uniform magnetic field has been profusely studied in the literature. The simultaneous presence of a staggered component in the same direction affects the phase diagram, but has been much less explored, beyond the initial work by Alcaraz et al. [107], and some recent results for particular values of the parameters [113]. Furthermore, as discussed in the next section, our simulation restricts the system to be in the sector of vanishing total z−component of the spin. This is necessary in order to employ duality properties to understand numerical results.

III. STRATEGY FOR NUMERICAL SIMULATIONS
A. Matrix product states In this paper, we are interested in the phase structure of the Thirring model at zero temperature. Hence, we need to search for the ground state of the system for different sets of parameters. For small system sizes, it is possible to find the ground state numerically, either by diagonalizing the full Hamiltonian or by targeting only the very lowest state, e.g. with the Lanczos algorithm. However, the main difficulty in this approach is that when one enlarges the system size, N , the dimension of the Hilbert space increases exponentially with N . This exponential explosion can sometimes be evaded by using tensor network methods [14,16,17,[114][115][116]. In particular, the MPS ansatz [10,12,[117][118][119] is especially adequate for one dimensional quantum many-body systems. For a one-dimensional lattice system with N sites, each of them with a d−dimensional local Hilbert space generated by {|σ i } d σi=1 , a generic quantum state can be written In a MPS, the basis coefficients c σ1···σ N adopt the particular form, where each M σ k is a D × D matrix. With open boundaries, M σ 1 and M σ N are instead D−dimensional vectors, and the trace reduces to the simple product of matrices.
The parameter D is called the bond dimension and determines the maximum amount of entanglement in the state. With respect to a bipartition of the chain in two regions A|B, the entanglement entropy is defined as the von Neumann entropy of the reduced state, where ρ A = tr B |Ψ Ψ| is the reduced density matrix for subsystem A (and analogously for B, ρ B = tr A |Ψ Ψ|) 7 . In the generic case, the entropy of A can be as large as the number of sites it contains (its volume). If the state is a MPS and A is a subchain, it is easy to see that S A|B ≤ 2 log D, i.e. it is bounded from above by a constant independent of the system size. The MPS ansatz hence satisfies by construction an area law [15].
Any state |Ψ can be exactly written in the MPS form with bond dimension D ≤ d N/2 . Alternatively, a MPS approximation with a maximum bond dimension D cut can in principle be found by performing a sequence of singular value decompositions (SVD) 8 across each bond of the chain, and retaining only the largest D cut singular values for each of them [12,118]. The distance to the true state can be bounded by the discarded weights, which thus give an estimation of the truncation error.
For the purpose of practical calculations, the representation is useful only if D is relatively small. This is in fact the case for states which satisfy certain entropic area law [121], and, in particular for ground states of local Hamiltonians with a gap [122][123][124]. Even in the gapless case, where the entropy presents logarithmic corrections to the area law [125], a MPS approximation is possible with bond dimension that scales only polynomially (not exponentially) with the system size [122].
These properties, together with the existence of several efficient algorithms for finding MPS approximations to the ground states of one-dimensional problems [14,16,17,115,116], make the MPS ansatz one of the most precise numerical methods for strongly correlated models in one spatial dimension.
Here we use a variational search to optimize the MPS that minimizes the energy with respect to the Hamiltonian (24). The algorithm, closely related to the density matrix renormalization group 9 (DMRG) [126,127], has a computational cost that scales as N D 3 (for open boundaries) 10 . The method proceeds by optimizing the d × D × D dimensional tensors M one by one, in a sequential manner, sweeping iteratively over the chain back and forth until the desired level of convergence has been reached in the energy. Namely, we choose to stop the iteration when the relative change in energy is smaller than the tolerance = 10 −7 . Because the optimization of a local tensor can only decrease the energy, the algorithm is guaranteed to converge (although it could be to a local minimum). This provides an approximation to the ground state of the model, which can be systematically improved by increasing the bond dimension.

B. Matrix product operator
The form of the Hamiltonian most suitable for numerical calculations is given in Eq. (24), representing a spin-1/2 system with a local Hilbert space of dimension two. Since the duality between the massive Thirring model and the sine-Gordon model is only valid in the zero-charge sector, the simulations should be performed in this sector, corresponding in the language of the spin model to S z tot = n S z n = 0. Within the MPS framework, it is possible to construct states that satisfy this constraint explicitly, by imposing the proper structure in the tensors M . In this work, instead, we use generic tensors, and introduce a penalty term in the Hamiltonian to ensure that the ground state is in the desired S z tot sector, namelȳ where the magnitude of λ should be chosen large enough to ensure that all the undesired states have energy above the lowest state of the targeted sector.
For efficient contraction of MPSs with an operator, such as the Hamiltonian, the operator has to be expressed also in the tensor network language, i.e. as the so-called Matrix Product Operator (MPO) [129,130]. We find that the MPO representation of Eq. (29) is given by 9 To our knowledge, there has not been extensive exploration of quantum field theories, as well as the XXZ spin model in Eq. (24) with the presence of a staggered magnetic field, using the original DMRG approach. 10 Variants of the algorithm exist also for periodic boundary conditions, in which case the cost scales with N D 5 [111,128].
for the tensors at the boundaries and in the bulk, where

C. Simulation details
We describe now our simulation strategy. The variational search begins with a randomly-initialized MPS with bond dimension D = 50. To have reliable results, we aim to observe convergence in D. With several different values of D, one can investigate the truncation error systematically, and extrapolate the physical quantities to the infinite-D limit (see Sec. IV B). Having results with D = 50, we gradually increase the bond dimension to 100, 200, 300, 400, 500 and finally 600. To do so, the size of the optimized tensors is increased to the desired D value, and the additional components are initialized to zero or a small random number. This MPS is used as initial guess for the variational procedure, which is run again until convergence. This is repeated, successively increasing the bond dimension, until our final D = 600 is reached.
Similarly, we study finite size effects, using four system sizes, N = 400, 600, 800, 1000, which allows us to perform an infinite volume extrapolation. We cover a wide parameter range to study the phase structure. The coupling ∆(g) is chosen from the range −0.9 ≤ ∆ ≤ 1.0, with five different masses,m 0 a = 0.0, 0.1, 0.2, 0.3, 0.4. To study the mass dependence in more detail, for some values of the coupling, we simulate also additional masses,m 0 a = 0.005, 0.01, 0.02, 0.03, 0.04, 0.06, 0.08, 0.13, 0.16. We set the parameter of the penalty term, λ, to 100, and target the zero-charge sector (S target = 0).
In performing the search of the ground state using the variational method, we observe that the convergence of the algorithm is slower in some regions of the parameter space, namely for m 0 a = 0 and, in the massive case, for large negative ∆(g). This is consistent with a regime where the theory may become critical. Figure 2 shows examples of these fast-and slow-convergence cases. For the slow cases, not only it takes more sweeps for the algorithm to converge, but also iterations of the Jacobi-Davidson solver used to solve for the local tensors are also significantly more time-consuming.

IV. NUMERICAL RESULTS FOR PROBING THE PHASE STRUCTURE
This section describes numerical computations for quantities that can be employed to probe the non-thermal phase structure of the Thirring model. As indicated by the perturbative RGE's, Eq. (7) and (8) in Sec. I, it is expected that there are at least two phases in the massive Thirring model, with theψψ operator in Eq. (1) being relevant in one of them and irrelevant in the other. Since thisψψ operator is dual to the cos(φ) term in the sine-Gordon theory in Eq.
(2), one envisages that in the regime whereψψ is irrelevant in the Thirring model, the corresponding bosonic theory is free. Furthermore, since we are investigating two-dimensional systems, and the sine-Gordon model is closely related to the XY model [95], it is foreseen that the phase transition in the Thirring model is of BKT type. Below, we show that these features in the phase structure can be observed using the ground state obtained with the technique of MPS described in Sec. III.
In the following, we demonstrate calculations for the entanglement entropy (Sec. IV A), the fermion bilinear condensate, ψ ψ (Sec. IV B), as well as density-density and fermion-antifermion correlators (Sec. IV C). Our results of these four objects can be employed to obtain knowledge of the phase structure that we summarize in the next section.

A. Entanglement entropy
To investigate the non-thermal phase structure of the massive Thirring model in 1+1 dimensions, we first study the von Neumann entanglement entropy for the ground state of the XXZ spin chain with finite spatial extent of size N (cf. Eq. (28)). Cutting the n-th link of the chain divides the chain in two parts, containing n, and N − n sites, for n = 0, 1, 2, . . . , N − 1. If the state is a MPS, the corresponding Schmidt decomposition gives at most D non-zero values, {s α (n)}, which can be easily recovered from the MPS representation [12]. The entanglement entropy with respect to this cut can thus be computed as Since the XXZ spin chain studied in this work is equivalent to the massive Thirring model, this entanglement entropy, computed in the spin model, is a useful tool to probe properties of the ground state in the latter. In particular, it can be employed to identify critical points in the field theory. As demonstrated by Calabrese and Cardy [125,131], at criticality, S N (n) exhibits the scaling behaviour, where c is the central charge, and k is a constant that can be argued to reflect the boundary effects. Figure 3 shows examples of the S N (n) from our analysis atm 0 a = 0. In both plots displayed in this figure, the critical Calabrese-Cardy scaling behaviour is manifest when the bond dimension is large enough. This property is present at all values of ∆(g) for vanishingm 0 a = 0, signalling that the system is at criticality. In Fig. 4, we exhibit results from performing the same scaling test form 0 a = 0.2 at three choices of ∆(g). In this case, as in results from all other simulations with non-vanishingm 0 a, the Calabrese-Cardy scaling is observed only for ∆(g) smaller than a certain value, ∆ * , which depends onm 0 a [∆ * = ∆ * (m 0 a)]. This is indicated by the left panel of Fig. 4. In the middle panel of the same figure, we show a case where ∆(g) > ∼ ∆ * . It is clear that the scaling behaviour is no longer present. Finally, when ∆(g) ∆ * (right panel of Fig. 4), the entanglement entropy, S N (n), is almost independent of n and its value is small. This gives strong hints that the theory is in a gapped phase at ∆ > ∆ * whenm 0 a = 0 11 . Furthermore, our data show a clear trend that ∆ * decreases with increasingm 0 a. In the regime of very smallm 0 a ( < ∼ 0.04), we see that ∆ * ∼ −1/ √ 2. This is what one expects from the small-mass RGE's, Eqs. (7) and (8), as discussed in Sec. I. The above features of the phase structure are confirmed by our analysis of fermion correlators, which will be discussed in Sec. IV C.
In the cases where the Calabrese-Cardy scaling behaviour is observed, one can use Eq. (34) to determine the central charge, c. For this purpose, we study the entanglement entropy, S N (n), as a function of X, where In Fig. 5, we display such an example for ∆(g) = −0.88 andm 0 a = 0.2, where the Calabrese-Cardy scaling is manifest. From this plot, it is obvious that at large enough bond dimension, one can read off the central charge and the result is c ∼ 1. Throughout the entire regime where the system is found to be critical, we extract the central charge by fitting S N (n) computed at the largest bond dimension (D = 600) in this work to Eq. (34). Such analysis finds that c is always at most ∼ 1% away from one. This is consistent with a free bosonic field theory at criticality, and is in accordance with the prediction of the RGE's, Eq. (5) and (6). That is, there exists a phase where the cos[φ(x)] operator in the the dual sine-Gordon model Eq.
(2) is irrelevant, so is theψψ operator in the dual massive Thirring model.
To summarize, we find strong evidence for the existence of two phases in the massive Thirring model. In one of the phases the Calabrese-Cardy scaling is valid, while in the other this scaling behaviour is absent. Through the detailed studies of the correlators, presented in Sec. IV C, it will be established that the theory is at criticality in the former, and the latter phase is gapped.

B. Fermion bilinear condensate
In order to obtain further information for the nature of the phase transition discussed in Sec. IV A, we investigate the chiral condensate,χ = aχ = a ψ ψ .
Under the Jordan-Wigner transformation, That is, the chiral condensate in the (1+1)-dimensional Thirring model corresponds to the staggered magnetization in the XXZ spin chain. In the spin chain model, the magnitude of anisotropy never exceeds one according to Eq. (25). This implies that a Néel state can only arise when the staggered magnetic field, am 0 = 0, is applied. For the corresponding quantum field theory, the Thirring model, this means that the chiral condensate is expected to be zero in the massless limit. Such a feature is consistent with the fact that the massless Thirring model in 1+1 dimensions is a conformal field theory. Furthermore, due to the presence of a uniform magnetic field at non-zero ∆(g) in the spin model, it is obvious that the value of the staggered magnetization (chiral condensate in the Thirring model) will exhibit ∆−dependence.
It is worth noticing that the Mermin-Wagner-Coleman theorem [97,98] is applicable for the theory under investigation. As also elaborated in Ref. [132], this means that the condensate defined in Eq. (36) is not a viable order parameter for phase transitions in 1+1 dimensions, and can be non-vanishing when the theory is in the critical phase.
We compute the chiral condensate at all values of ∆(g) and am 0 where numerical calculations are performed. At every choice of [∆(g), am 0 ], we first estimateχ in the limit of infinite bond dimension. Although there are no errors on the "raw data" forχ at finite D, we employ the following procedure to assign errors on infinite−D results of the condensate.
1. Using the raw data ofχ at D = 500 and 600,χ 500 andχ 600 , a "linear extrapolation" to the D → ∞ limit is first performed. Result for the chiral condensate from this step is denoted byχ 2. Estimate the central value of the infinite−D condensate,χ ∞ , by computinĝ . Notice that errors on the data points and the extrapolated result for the right panel are too small to be discernible on the plot.

The (symmetric) error onχ in the limit of infinite bond dimension is evaluated by
4. If the error, δχ ∞ , obtained in the previous step is smaller than the chosen precision of the variational algorithm in this work, = 10 −7 , we replace the numerical value of δχ ∞ by .
It is worth noticing that the above procedure is a conservative approach for assigning errors toχ in the infinite−D limit,χ ∞ . This is because our data show that the chiral condensate converges very well at D ≥ 400.
Results ofχ in the infinite−D limit are then used for the procedure of taking the thermodynamic limit, N → ∞. For this extrapolation in N , we use data points at N = 400, 600, 800 and 1000. Figure 6 shows examples of extrapolations in D and N for [∆(g), am 0 ] = [−0.9, 0.01]. In this figure, results of the condensate obtained at D = 200, 300, 400, 500, 600 and N = 400, 600, 800, 1000 are exhibited. From these plots, one can easily observe thatχ depends very mildly on D for D ≥ 400, and it also extrapolates smoothly to N → ∞ using data at N ≥ 400. Such a feature is in fact observed for all choices of [∆(g), am 0 ] in this work. As mentioned above, this means that our method of assigning errors toχ ∞ is a conservative approach.
We also remark that finite volume effects (FVE) in the Thirring model were investigated a long time ago in Ref. [133].
The authors compared an exact solution for the infinite-volume mass spectrum with results of Monte Carlo simulations at finite volume and concluded that FVE are suppressed to the level of a few percent when N > ∼ 50. In our case, despite using different observables, we observe O(1%) FVE at N ∈ [400, 1000], which is qualitatively consistent with the results from Ref. [133]. According to results of the entanglement entropy discussed in Sec. IV A, the theory is in the gapped (massive) phase at am 0 = 0 for ∆ > ∆ * where ∆ * is mass-dependent, while it can be in the critical phase at ∆ < ∆ * . This feature will be further confirmed by the investigation of fermion correlators, presented in Sec. IV C. We notice that, in Fig. 7,χ is non-vanishing at ∆(g) = −0.8 for all data points with am 0 = 0. Through the study of the entanglement entropy reported in Sec. IV A, we also know that for ∆(g) = −0.8, the theory can be in two different phases separated at am 0 ∼ 0.2. In other words, the chiral condensate can be non-zero in both phases. This means that χ is not an order parameter for the observed phase transition, providing further evidence that this transition can be of the BKT-type [132]. In Fig. 7, it can be seen thatχ extrapolates smoothly to zero at vanishing am 0 . As mentioned above, this is in accordance with the fact that the massless Thirring model in 1+1 dimensions is a conformal field theory. Furthermore, we find that the chiral condensate computed directly at am 0 = 0 is zero for all values of ∆(g). Given that all simulations that lead to results in Fig. 7 are performed at finite system sizes, we carry out checks for am 0 = 0 calculations with infinite-size simulations by employing the variational uniform MPS (VuMPS) method [134]. These checks confirm thatχ obtained from simulations at am 0 = 0 indeed vanishes. Results of the VuMPS approach will be published in a separate article where we will report our study of real-time dynamics associate with the BKT phase transition in the massive Thirring model [89].

C. Correlation functions
We now proceed to present results for the correlation functions in the Thirring model. We will consider two types of correlators: density-density (which in the continuum is ψ (x 1 )ψ(x 1 )ψ(x 2 )ψ(x 2 ) , and in terms of the discrete fermions corresponds to c † n c n c † m c m /a 2 ) and fermion-antifermion ( ψ (x 1 )ψ(x 2 ) or c † n c m /a 2 ). The connected part of the former can be expressed in the spin language as: The sum over n indicates that, for a given distance x, we average over all pairs of spins within the 200-site subchain in the middle of the lattice, to avoid boundary effects. Despite using open boundary conditions, translational invariance is realized almost ideally in this region. We take x to be odd, i.e. we look at the correlator between even-odd or odd-even sites. The number of pairs corresponding to a given distance x = 1, . . . , 199 is denoted by N x and is equal to 200 − x. To calculate the connected part of the correlator, we subtract the product of single-site expectation values.
The fermion-antifermion correlation function is expressed in terms of the spin operators as the following string correlator: with the same subchain averaging procedure as for the distance-dependent part of C zz (x). This correlator can be shown to correspond to the soliton-soliton correlation function in the dual sine-Gordon theory [91].
We first concentrate on the short-distance behaviour of both correlation functions. The decay of C zz (x) and C string (x) is expected to be power-law in the critical phase and power-exponential in the gapped phase. To test this behaviour, for all parameter values and for both types of correlators, we perform fits using the power-law ansatz: with the fitting parameters α, β, and the constant C allowed only for C string (x). Moreover, we fit three types of exponential ansatzes: with fitting parameters A, B, η or A i , B i (we assume A 1 > A 2 > ...), and the constant C allowed in the case of the string correlator. Even if the true behaviour is power-law or mixed power-exponential decay, a MPS with finite bond dimension can only approximate it by a finite sum of exponentially decaying terms. Their decay rate, in a translationally invariant case, is determined by the eigenvalues {λ i } of the transfer matrix E ≡ σ M σ * ⊗ M σ , more concretely by the ratios λ i /λ 1 , where λ 1 is the largest one. Only in the limit of infinite D would it be able to reproduce a true power law at arbitrarily long distance. Thus, we also test how well a multi-exponential ansatz can describe the observed data. The parameters A or A i will be related to the aforementioned ratios of transfer matrix eigenvaluesin particular A is the ratio of the second largest to the largest eigenvalue, λ 2 /λ 1 . Thus, in the critical phase (strictly conformal in the massless case and in the continuum), A = 1 and the power-exponential fitting ansatz becomes the power-law one. A discussion is in place also for the constant C in the fermion-antifermion correlator. Such a constant cannot emerge in the continuum theory. However, on the lattice, discretization effects may cause the constant to appear. In the critical phase, all lattice effects should vanish in the infinite volume limit, since the theory is conformal in this regime (the mass operator is irrelevant, as discussed above in the context of RG equations). Thus, the constant should be zero, which we check explicitly. In the gapped phase, however, the mass operator becomes relevant and an energy scale appears. This breaks the conformality, and the continuum limit requiresm 0 a → 0. In such a situation, a constant may be generated in the string correlator. This constant is expected to increase the further away from the continuum limit one simulates -i.e. when the fermion mass or ∆(g) is increased. At large mass and deep in the gapped phase, a string order configuration may be favoured as an artefact of the staggered discretization. In spin language, the system is subject to a staggered magnetic field, leading to a Néel-type phase when this field is strong. Inspecting the definition of this correlator, Eq. (41), it becomes clear that this renders a constant at large distances.
In Fig. 8, we show example fits for the density-density correlator, for a system size of 1000 sites, coupling ∆(g) = −0.7 and a small fermion massm 0 a = 0.02. All of the results of this subsection are obtained in the limit of infinite bond dimension, using a procedure analogous to the one described in Sec. IV B. The left panel presents the power-law fit of Eq. (42) and the right panel the three exponential-type fits, Eqs. (43)- (45). In both cases, the fitting range is distances between x = 21 and x = 39, i.e. 10 consecutive data points. For the case depicted in these plots, we expect that the system is close to the point of the BKT transition and hence C zz should be well described by the power-law fitting ansatz. Indeed, as the left panel of Fig. 8 demonstrates, such a fit provides an excellent description of the correlator decay, following a straight line on a log-log plot. Obviously, the power-exponential fitting ansatz works equally well, since it boils down to the power-law functional form when A = 1, as reproduced by the fit. A single exponential is not enough to describe the data, as clear from visual inspection of the right panel -there is clear curvature visible with only the y-axis in logarithmic scale. However, two-exponential and three-exponential fits provide good description of the data. Nevertheless, they do not reproduce the expected power-law behaviour of the critical phase, i.e. the fitting parameter A 1 is different from one (although it increases towards one when more exponentials are included). The rate of the power decay of the correlator in the critical phase is described by the parameter α, equal to η for A = 1. The expected value of α is −2, well reproduced by the fit [135,136].
To illustrate how the gapped phase manifests itself in the considered correlator, we show the fits also for ∆(g) = 0.4 (other parameters unchanged), see Fig. 9. Even though the decay looks approximately linear in the log-log plot (left panel), our precision is good enough to conclude that power-law fit does not describe the data well, indicated by the large value of χ 2 /dof. The mixed power-exponential ansatz works well, with A ≈ 0.923 and η ≈ −1.54. Thus, the spectrum of the transfer matrix is gapped and the exponent η significantly deviates from −2. We note that this picture is consistent with the one from entanglement entropy, i.e. the phase where power-law fits are proper is the one where the central charge is found close to 1, while the region where the exponential correction to the correlator  decay is important manifests itself by flat behaviour of the entanglement entropy, i.e. the Calabrese-Cardy scaling is not observed.
The correlators are expected to follow a given type of behaviour asymptotically, i.e. at large enough distances. Hence, a proper choice of the fitting interval has to be made. To analyse the dependence of the fitting parameters on the coupling ∆(g), we avoid the arbitrary choice of the fitting interval by adopting a systematic procedure, similar to the one used e.g. in Ref. [18] (see the Appendix of this reference). We consider all possible fits in the interval x ∈ [5, 49] encompassing a minimum of 10 consecutive distances 12 . Each fit is weighted with exp(−χ 2 /dof) and we build histograms of the fitting parameters for each fitting ansatz. The central value for each fitting parameter in a given physical setup (same system size, fermion mass and coupling) is extracted as the median of this distribution and the error as half of the interval in which 68.3% of the weighted fits around the median are contained (corresponding to a 1-σ deviation in the case of an ideally Gaussian distribution). We note the obtained distributions are approximately Gaussian and the thus extracted systematic error is in most cases a factor 5-10 larger than the error obtained from typical fits. Finally, we add this systematic error to the one of the fit that best describes the data, defined as the one with the smallest error among the fits with χ 2 /dof ≤ 1.
The result of applying this procedure to the density-density correlator is shown in Fig. 10, again for fermion mass m 0 a = 0.02. In the left panel, we show the A parameters extracted for different couplings. We note the A parameter of the power-exponential fit becomes compatible with 1 somewhere between ∆(g) = −0.4 and −0.6. The expected location of the BKT crossover is at ∆(g) ≈ −0.7, however the smooth transition between the functional forms of the power-law and power-exponential type of behaviour makes it impossible to locate the BKT point at the current level of precision. We expect that close to the critical point, there is only a small admixture of the exponential factor to the power-law term, impossible to disentangle without much better precision. Further into the gapped phase, at ∆(g) > ∼ −0.4, the exponential term becomes clearly visible and A is no longer consistent with 1. The value of A drops when ∆(g) is increased and the exponent η of the power-law factor in the fitting ansatz increases towards less negative values. Thus, the exponential decay becomes relatively more important deeper in the gapped phase. This is also indicated by the smaller difference of the parameter A and the parameter A 1 of the 3-exponential fit for positive ∆(g), which would agree in the limit of purely exponential behaviour.
In Fig. 11, we show the same kind of plot for a larger fermion mass,m 0 a = 0.3. In this case, the dependence of the parameter A of the power-exponential fit is much steeper and A becomes compatible with 1 between ∆(g) = −0.82 and −0.84. This signals that the BKT transition moves towards more negative values of the coupling with increasing fermion mass. For ∆(g) > −0.82, the system is clearly in the gapped phase, which is indicated also by χ 2 /dof 1 for the pure power-law fits. In contrast, such fits in the small fermion mass case are still reasonable until ∆(g) ≈ 0, as a consequence of our rather conservative error estimate procedure. We therefore conclude that the BKT crossover is more pronounced for larger fermion masses. However, interestingly, the value of the exponent η of the powerexponential fit is consistent with −2 for all couplings. A comparison of the coupling dependence of the parameter A for different fermion masses,m 0 a = 0.005, 0.02, 0.08, 0.3 is shown in Fig. 12. It offers a rather clear picture -while for the smallest mass the parameter A is consistent with 1 up to ∆(g) ≈ 0, increasing the fermion mass confines the location of the BKT transition to smaller and smaller ranges of the coupling. Nevertheless, only in the largest mass case in this plot, it is possible to determine the transition point with a precision of a few percent. It is very likely that the critical coupling, expected to be around g * = −π/2 (∆(g * ) ≈ −0.7), moves towards more negative values as the mass is increased, but the numerical evidence for this is convincing only for the largest mass, at our level of precision in this correlator. For smaller fermion masses, their effect can be interpreted as a perturbation of the massless conformal theory. The mass perturbation opens the gap at the transition point, but the almost power-law decay of the correlator for small fermion masses is a remnant of the conformal phase that extends considerably into the gapped phase.
We also consider the fermion-antifermion (string) correlator. In Figs. 13, 14, we show its decay for the same parameters as in Figs. 8, 9 (N = 1000,m 0 a = 0.02 and two couplings, ∆(g) = −0.7 and ∆(g) = 0.4). We illustrate our fits again with examples for the fitting range of x ∈ [21, 39] and we reach conclusions that are consistent with the ones from the density-density correlator. For coupling close to the BKT transition point, the decay is power-like, indicated by the good χ 2 /dof of the power-law and power-exponential fits, the latter having A close to 1 and thus being dominated by the algebraic term. The constant C of the power-exponential fit is small, but non-zero, indicating that the actual transition may occur for slightly more negative ∆(g). As we concluded above for the density-density correlator, it is challenging to locate the transition point more precisely and it would require much more accurate data. When the coupling is increased and the system is clearly in the gapped phase, the pure power-law fit does not describe well the data and the proper fitting ansatz is the power-exponential one, that can be mimicked by a few exponentials. A rather large value of the constant C is observed and, similarly as for the density-density correlator, the value of A is much smaller than 1 and decreases towards larger ∆(g), as can be seen in Fig. 15. The exponent η increases from below −1 in the critical phase towards less negative values deep in the gapped phase, indicating that the fits are dominated by the exponential term in the latter. An important difference with respect to the C zz correlator is that the exponent α = η in the critical phase does not have a universal value. In the density-density correlator, α = η = −2 for all couplings in the gapless phase, while in the string correlator, there is rather strong dependence of this exponent when moving deep into this phase [88]. For larger fermion masses, the qualitative picture is the same as for them 0 a = 0.02 case. We illustrate the g-dependence of the fitting parameter A of the power-exponential fit in the left panel of Fig.  16. Comparing with the analogous plot for the density-density correlator, Fig. 12, we observe that the value of A for a given pair (m 0 a, ∆(g)) is basically independent of the correlator type. The right panel of Fig. 16 shows the coupling  dependence of the constant C. Similarly to the case of the parameter A, there is a pronounced dependence of C for the largest fermion mass -its value is zero in the critical phase and steeply increases when moving deeper into the gapped phase. For smaller masses, this dependence is milder, again as for the A parameter. However, in general, the relative systematic error of C is smaller than the one of A, allowing to pinpoint the transition point a bit more precisely.

V. DISCUSSION
In this section, we will use the numerical results reported in Sec. IV to discuss the phase structure of the massive Thirring model in 1+1 dimensions. Implication for the scaling behaviour and the continuum limit of the model is also addressed. From the investigation of the entanglement entropy and two types of correlators, as presented in Secs. IV A and IV C, we observe concrete evidence for the existence of two phases in the massive Thirring model in 1+1 dimensions. In one of these two phases we find power law decaying correlations and critical scaling of the entanglement entropy, indicating a conformal phase with central charge equal to unity. The other phase is a massive (gapped) phase: it exhibits exponentially decaying correlations and bounded entropy. Combining this with the results from Sec. IV B, where we showed that the chiral condensate does not vanish in either phase when m 0 a = 0, and thus cannot serve as order parameter, it can be concluded that the relevant phase transition is of the BKT-type. Although this phase structure is expected from duality properties of the model, our current work adds important new ingredients to this research direction. Our results demonstrate that the MPS methods are still applicable to study the lattice version of a model with this feature, despite the intrinsic difficulty of locating BKT transitions numerically. This leads to possibilities for further understanding physics of the BKT phase transitions, since the MPS formulation can naturally be implemented for studying real-time dynamics [89].
Our numerical calculation confirms that the (1+1)-dimensional Thirring model is a conformal field theory in the massless limit. At non-vanishing values of the bare mass in lattice units,m 0 a, both the massive and the critical phases are seen. To depict the phase diagram, we resort to the resultant central value of the constant term, C, in fitting the fermion-antifermion correlator defined in Eq. (41) to the ansatz in Eq. (43) which describes the powerexponential analysis of the data. It is reported in Sec. IV C that this ansatz always leads to good reduced χ 2 in the fitting procedure. At criticality, its results also converge well to those given by the power-law analysis using Eq. (42), with η compatible with α, and the parameter A in Eq. (43) being consistent with one. As already discussed in Sec. IV C, the parameter C is expected to be non-zero in the gapped phase where the string order in the corresponding spin chain can emerge, while it vanishes in the critical phase. The reason for choosing C as the main probe of the phase structure is because it is the best determined quantity amongst all the observables computed in this work. Figure 17 shows results for the central value of C at all choices of [∆(g),m 0 a] in our simulations. It is clear from this plot that on the ∆(g)−m 0 a plane, there is a region where this parameter is well consistent with zero. Here we also remind the reader that typical error on C is of percentage or subpercentage level.
We then use the information presented in Fig. 17 to extract the non-thermal phase structure of the (1+1)-dimensional massive Thirring model. While C is the most accurately computed quantity in our analysis procedure, it still contains systematic error. For cases where this parameter is small but non-vanishing in the fermion-antifermion correlator, it is challenging to decide whether the theory is in the conformal or gapped phase. In fact, we discover that when the central value of C is between 0.001 and 0.01, it is impossible to use our data for a clear judgement on this issue, as all the fit ansatzes in Eqs. (42), (43), (44) and (45) lead to acceptable reduced χ 2 . That is, power-law, power-exponential and pure-exponential functions all describe our data well when 0.001 ≤ C ≤ 0.01. For this reason, we identify the cases in this range as being "undetermined". It may be possible to reduce the size of the undetermined region by means of better (e.g. with larger bond dimension) simulations. Nevertheless, such a high-precision determination of the phase boundary is beyond the scope of this work.
From the above discussion, we classify our data points into three groups according the to the central value of the parameter C in fitting the data of fermion-antifermion correlator to Eq. (43).
1. The critical regime is where C < 0.001.
2. The gapped phase is identified as where C > 0.01.
Using this categorization, we plot the non-thermal phase structure of the Thirring model in Fig. 18. In this figure, there is a grey area. This area indicates qualitatively the regime where we find the undetermined points described above. It is obvious that the BKT transition occurs in this grey region, with the phase boundary described by a function where g * (m 0 a) is the value of g for the BKT phase transition in the massive Thirring model with bare mass m 0 and lattice spacing a. At a given value ofm 0 a, the theory is conformal when ∆ is less than ∆ * (m 0 a), and becomes gapped otherwise.
With the definition, FIG. 18: Non-thermal phase structure of the massive Thirring model from our numerical investigation. In addition to the data points that can be identified to be in the gapped phase (blue stars) or at criticality (red circles), there are points (black squares) where our simulations cannot determine which phase the theory is in. The grey area indicates the regime where we find these "undetermined" point. The BKT phase transition must occur within this grey area. Figure 18 shows that∆ * ∼ −0.7 .
As discussed in Secs. I and II, we work with the value of the four-fermion coupling, g, in the range between −π and π such that the S-duality is valid. In this range of the coupling, ∆ is singled-valued and increases monotonically with g, as can be seen from Eq. (25). Therefore, Eq. (48) implies that which is consistent with the result, Eq. (11), extracted from the continuum RGE's. It can also be inferred from Fig. 18 that ∆ * [g * (m 0 a)] and g * (m 0 a) decrease as the value of the fermion mass grows. This qualitative feature agrees with the that of the phase boundary presented in Fig. 1. Here we stress that Eq. (11) and Fig. 1 are obtained from a perturbative expansion in m/Λ, while results presented in this section are extracted non-perturbatively.

B. Implication for scaling behaviour and the continuum limit
Numerical results presented in this paper can be employed to infer scaling behaviour of the massive Thirring model in a non-perturbative fashion. Such investigation is particularly interesting for the critical phase. As already discussed in Sec. V A, with the phase boundary introduced in Eq. (46), it is found that the theory with non-zero bare mass is conformal in the regime ∆(g) < ∆ * (m 0 a) at a given value ofm 0 a. Because of chiral symmetry, the bare mass is proportional to the renormalized mass in the Thirring model 13 . For a massive theory to be conformal, the only possibility is that the fermion mass is an irrelevant coupling. This is consistent with the scaling behaviour in the small-mass limit as predicted by Eq. (8), and one can use the RG-flow diagram of Fig. 1 to illustrate this feature.
In this diagram, there is a stable "fixed half-line" at (g <ḡ * , m = 0) whereḡ * = −π/2. Let us denote the phase boundary (the dotted blue curve) in Fig. 1 by a function g * (m). The whole region spanned by g < g * (m) at each value of positive m is then the critical surface of this fixed half-line, where there is no relevant scaling variable.
Corresponding to the above scenario, theψψ operator is irrelevant in the critical phase of the massive Thirring model. That means, although this operator is of dimension-one at the classical level, quantum effects from the four-fermion interaction must render it into being irrelevant (with dimension greater than two) in the critical phase. In other words, a large anomalous dimension has to be generated through interactions. In Refs. [104,137] the scaling dimension of this operator was studied using the operator product expansion (OPE), and was found to increase with decreasing coupling strength, g. It is allowed to be larger than two when g is small enough. In this regard, our numerical investigation shows the same feature as that from these previous studies. It would also be possible to have further detailed numerical examination for the scaling properties of theψψ operator using the MPS method. Nevertheless, this calculation requires a dedicated project, and is beyond the scope of our current work.
In the gapped phase, theψψ operator is relevant, and a spectrum of excited-state masses can be generated. We denote the set of these masses by {M i }, and assume the hierarchy These masses can be computed using variational methods in the framework of MPS [138], with the tools developed in Ref. [18] for the study of the Schwinger model. With the fermion mass approaching zero, m → 0, all M i must vanish because the massless Thirring model is scale-invariant. Nevertheless, we observe that in the regime where the fermion mass m is well below the cut-off scale, the gapped phase can be interpreted as a mass-deformed conformal field theory, and conformality can be restored at any value of the four-fermion coupling on the unstable fixed half-line, (g >ḡ * , m = 0), in Fig. 1. In this case the scaling behaviour of M i can be studied by employing the hyperscaling techniques detailed in Ref. [139]. To proceed, we first choose a conformal point, (g =ḡ + , m = 0) whereḡ + is greater thanḡ * = −π/2. This conformal field theory is then deformed by introducing a small mass perturbation (ḡ + , m deform ), where m deform in lattice units is much smaller than unity. In other words, this mass-deformed theory is in the vicinity of the fixed point (g =ḡ + , m = 0). With the condition that the four-fermion coupling is not changing under a RG transformation, it is then straightforward to show that [139] where γ(ḡ + ) is the anomalous dimension of theψψ operator evaluated atḡ + , and c i is function ofḡ + . In other words, when the model is very close to a particular fixed point, all the excited-state masses will scale to zero with the same exponent. This interesting behaviour can be tested in future numerical computations for the spectrum of the model [138]. It also offers an approach to determine the anomalous dimension, γ(ḡ + ), in the gapped phase near the fixed half-line. We stress again that Eq. (51) is obtained with the assumption that g takes the value ofḡ + and does not change under a RG transformation.
The above scaling properties can serve as a guide to understanding the continuum limit of the massive Thirring model in 1+1 dimensions. By regarding the field theory as a statistical mechanics system, the latticized model approaches its continuum counterpart at criticality where the correlation, ξ, diverges, One can then associate ξ with typical scale at low energy, e.g., the fermion mass, m. This immediately leads to the conclusion that the theory is in the continuum limit on the linem 0 a = 0 in Fig. 18. However, because of the BKT phase transition, complication arises in the Thirring model. This is because the line,m 0 a = 0, is divided into two sectors which are associated with different scaling properties discussed above. For the critical area with non-vanishing mass in Fig. 18, the model describes the same physics as that on the stable fixed half-line, (g <ḡ * ,m 0 a = 0). This whole critical regime withm 0 a = 0 is then a continuum limit where the theory remains gapless/massless in the IR.
In the gapped phase, approaching the continuum limit, the unstable fixed half-line described by (g >ḡ * ,m 0 a = 0), is less straightforward. Below we address the related issues.
Since the extrapolation to the continuum limit requires the knowledge of how the theory scales when changing the

VI. CONCLUSION AND OUTLOOK
In this paper, we report our study for the non-thermal phase structure of the massive Thirring model in 1+1 dimensions employing the MPS approach. By restricting the model in the sector of vanishing total fermion number, it is expected to be S-dual to the sine-Gordon theory where classical soliton solutions are known. Our construction of the Hamiltonian makes use of the staggered-fermion discretization in the spatial direction. Through the application of the JW transformation, this Hamiltonian can be shown to be equivalent to the XXZ spin chain coupled to uniform and staggered external magnetic fields. In addition, a "penalty term" is included in the Hamiltonian [Eq. (29)] to lift the energies of other sectors to the cut-off scale. This ensures that the variational optimization of the MPS with the DMRG method will result in the ground state with zero total fermion number 15 . Numerical simulations are then performed at fourteen values of bare fermion mass, ranging from 0 to 0.4 in lattice units, and at twenty four choices of the four-fermion coupling, g, straddling the regime where phase transitions are expected to occur. Furthermore, we carry out computations at seven different bond dimensions (D = 50, 100, 200, 300, 400, 500, 600), and four system sizes (N = 400, 600, 800, 1000). This facilitates the extrapolation to the limit of infinite D, as well as to the thermodynamic limit.
Results of this work clearly identify the existence of two phases, one critical and the other gapped, in the (1+1)dimensional massive Thirring model. We also demonstrate that in the critical phase, the model is equivalent to the free bosonic field theory with the central charge being unity. Through the study of the entanglement entropy, theψψ condensate, as well as the density-density and the fermion-antifermion correlators, we observe unambiguous numerical evidence that the two phases are separated by a BKT transition. The theory in the chiral limit is found to be in the critical phase. This is in accordance with the expectation that the two-dimensional massless Thirring model is a conformal field theory. For the case of non-vanishing fermion mass, m, both the critical and the gapped phases appear. The phase boundary, as depicted in Fig. 18, can be specified with a function, g * (m), that takes decreasing value against growing m. The critical phase is in the regime g < g * (m) at a given m. Furthermore, our calculation shows thatḡ * ≡ lim which is consistent with previous studies.
The investigation presented in this article provides information for the scaling behaviour and the control of the continuum extrapolation in the massive Thirring model. One intriguing feature of the theory with non-zero bare fermion mass is that it can be in the critical phase, as indicated in Fig. 18, where theψψ operator becomes irrelevant although its classical dimension is one. This implies the generation of a large anomalous dimension through the fourfermion interaction. Detailed numerical examination of the scaling dimension can be carried out, although it requires a dedicated implementation and is beyond the scope of current work. Another interesting character of the massive Thirring model is that there are two distinct types of continuum limits. The first kind is defined on the stable fixed half-line of g <ḡ * on the m 0 a = 0 axis of the g−m 0 a plane, plus the entire region of the critical phase with m 0 a = 0. The second kind of continuum limit is the unstable fixed half-line, (g >ḡ * , m 0 a = 0). In Sec. V B of this article, we argue that this half-line contains an infinite number of possible continuum limits which can be distinguished by ratios of the excited-state masses extrapolated to the limit m 0 a → 0 in the gapped phase.
Our work described in this paper establishes the possibility of applying the MPS method for probing the BKT phase transition in quantum field theories. It is also the first step towards further exploration of the massive Thirring model using this approach. As mentioned above, it is interesting to perform detailed studies of the scaling behaviour and the continuum limit of the model, which requires the computation of excited-state masses [138]. In addition, the MPS offers another handle of the scaling properties, through examining how the theory responds to the change of the bond dimension [141][142][143][144]. Finally, it is also essential to understand the dynamical aspects of the BKT transition in the massive Thirring model, since it can shed light on the dynamics related to the process of engineering a field theory into a conformal phase. In this regard, the MPS technique provides a natural way to implement real-time evolution with a quench across the phase boundary, and this is one of our next steps in this research programme [89].
Extending the study to higher dimensions is also possible in principle. The natural generalization of MPS to two (or more) spatial dimensions is the family of projected entangled pair states (PEPS), introduced in Ref. [13]. Variational algorithms can also be employed to find a PEPS approximation to the ground state of a given Hamiltonian (see e.g. recent results in Ref. [145]). However, the computational cost involved is still considerably higher than in the one-dimensional case. Therefore, carrying out such an exhaustive exploration of the phase diagram in 2+1 dimensions is presently a longer-term goal that will require optimizing or adapting the algorithms.