Dense $\textrm{QCD}_2$ with matrix product states

We study one-flavor $\mathrm{SU}(2)$ and $\mathrm{SU}(3)$ lattice QCD in ($1+1$) dimensions at zero temperature and finite density using matrix product states and the density matrix renormalization group. We compute physical observables such as the equation of state, chiral condensate, and quark distribution function as functions of the baryon number density. As a physical implication, we discuss the inhomogeneous phase at nonzero baryon density, where the chiral condensate is inhomogeneous, and baryons form a crystal. We also discuss how the dynamical degrees of freedom change from hadrons to quarks through the formation of quark Fermi seas.


Introduction
The study of the phase structure of cold dense QCD is of great importance to understand the physics inside neutron stars as well as to deepen our understanding of QCD as a part of the fundamental theory of our universe [1,2].In dense QCD matter, the physics of confinement and colored Fermi seas of quarks play an essential role in determining the phase structures.Indeed, various interesting phases have been conjectured so far such as the quarkyonic matter [3][4][5], and color superconductivity [6] (see e.g., refs [7][8][9] for reviews of the QCD phase diagram).However, the direct study of cold dense QCD is very difficult except at large chemical potential, where interactions between quarks become weak due to asymptotic freedom, and perturbative QCD is applicable.In the presence of the baryon chemical potential, which is not high enough to accommodate color superconductivity, the ab-initio study of QCD is demanded, but conventional Monte Carlo simulations based on lattice QCD are no longer valid due to the notorious sign problem; the fermion determinant in the path integral becomes complex at nonzero baryon chemical potential, so that the importance sampling breaks down [10].Several methods have been proposed to circumvent the sign problem in the path integral formulation [11,12].
Recently, the Hamiltonian lattice gauge theory is reconsidered with the help of tensor networks and quantum computing to attack the problems suffering from the sign problem such as finding groundstate at finite-density or in the presence of the θ term, and real-time problems [13,14].In particular, numerical simulations based on a tensor network known as matrix product state (MPS) [15] are found to be effective to solve (1+1)-dimensional QED, which is also known as the Schwinger model [16][17][18][19][20][21][22][23][24] (see e.g., refs.[25,26] for applications of the tensor network in the path integral formulation of the Schwinger model).Variational computations based on MPSs have also been adopted to (1 + 1)-dimensional QCD, which we refer to as QCD 2 for SU(2) [27][28][29][30], and SU(3) [31,32].However, to our knowledge, comprehensive study of QCD 2 at finite density has not been done yet.
In this paper, we study one-flavor SU(2) and SU(3) QCD 2 at zero-temperature and finitedensity based on Hamiltonian lattice QCD and matrix product states [15].We adopt the formulation developed in ref. [30] to decouple nonabelian gauge fields explicitly in open boundary conditions.We use a variational approximation for the wave function of the ground state of quark many-body systems interacting with color coulomb force based on density matrix renormalization group (DMRG) [33,34], and compute physical observables such as the equation of state, chiral condensate, and quark distribution function as a function of the baryon number density, where the conventional lattice QCD simulation suffers from the sign problem.As a physical implication, we discuss the inhomogeneous phase at nonzero baryon density, where the chiral condensate is inhomogeneous, and baryons form a crystal.We also discuss how dynamical degrees of freedom change from hadrons to quarks through the formation of quark Fermi seas.We find that the study of dense QCD 2 is a good playground to test theoretical tools developed for ab-initio study of cold-dense states in QCD or QCD-like theory.
The remainder of this paper is organized as follows.In section 2, we review the Hamiltonian formulation of SU(N c ) lattice QCD in (1 + 1) dimensions.We eliminate the gauge field from the Hamiltonian by a unitary transformation under open boundary conditions in section 2.2.We map our system to the spin system through the Wigner-Jordan transformation [35] in section 2.4.In section 3, we introduce a matrix product state and matrix product operator, and explain the DMRG technique briefly.In section 4, we show the numerical results.We used iTensor [36] for numerical calculations.Section 5 is devoted to the summary.In Appendix A, we show several results of the free theory with open boundary conditions to compare with QCD 2 .
2 QCD 2 on a lattice

Hamiltonian formulation
We consider SU(N c ) gauge theory with N f flavors in (1+1) dimensions, although our numerical calculations mainly deal with the cases of N c = 2, 3 and N f = 1.We employ the Kogut-Susskind Hamiltonian with staggered fermions [37], which is composed of three parts: electric H E , hopping H hop , and mass H m parts: H = H E + H hop + H m , where Here, g 0 , a, m 0 , and N are the bare gauge coupling, lattice spacing, bare mass, and number of sites, respectively.χ and χ † are annihilation and creation operators of N c × N f component fermions, U is the link variable with the fundamental representation that is a N c × N c matrix, and E 2 i (n) is the square of the electric field.The index "i" denotes the color of adjoint representation (i = 1, • • • N 2 c −1), and repeated indices are summed unless otherwise noted.For convenience in later numerical calculations, we introduce adimensional Hamiltonians defined by Here, J = ag 0 /2, w = 1/(2ag 0 ), and m = m 0 /g 0 are the dimensionless couplings.Hereafter, we use the unit system with g 0 = 1 and omit g 0 where it does not lead to confusion.The fermion fields satisfy the anticommutation relation: where c and f represent the color and flavor indices, respectively.In the Hamiltonian formulation on a lattice, there are two types of electric fields, R i (n) and L i (n), which generate the gauge transformation of U from the right and left, respectively.Their commutation relations are as follows: ) ) ) Here, T i are generators of fundamental representation that satisfy the commutation relation, with the structure constant f k ij .R i (n) and L i (n) are not independent but related by a parallel transport: where [U adj (n)] j i is the link variable with the adjoint representation defined by which satisfies [U † adj (n)] j i = [U adj (n)] i j .From eq. (2.13), the square of R i (n) and Therefore, H E is independent of the choice of electric fields.
A physical state |Ψ⟩ is gauge invariant, which needs to satisfy the Gauss-law constraint: where G i (n) is defined by with the color-charge density operator, (2.17) Q i (n) satisfy the following commutation relation as do electric fields: Note that the Hamiltonian is gauge invariant, so that the Hamiltonian commutes with G i (n), [H, G i (n)] = 0.
(2. 19) This can be accomplished if the unitary operator Θ satisfies ) The construction of Θ is as follows: We can express the link variable U (n) by introducing a gauge field A i (n) as Using this gauge field A i (n), we introduce the following operator: which satisfies The unitary operator Θ satisfying eqs.(2.20) and (2.21) is constructed from V as Using Θ, we eliminate the link variables from the Hopping term H hop , The consequence of this transformation is that the electric part of the Hamiltonian, H E , becomes nonlocal and more complex.However, the Gauss-law constraint becomes simpler and can be easily solved under open boundary conditions.The right-electric field R i (n) and the color-charge density Q i (n) transform under V as and which lead to (2.29) The transformation of the left-electric field L i (n) can be obtained from eqs. (2.29) and (2.13) as (2.30) Using these transformations, the Gauss-law constraint G i (n) for n > 1 reduces to On the other hand, for n = 1, it is Here, L i (0) represents the incoming flux to the system.In the following, we impose no incoming flux condition, L i (0) = 0; then we have The Gauss-law constraint is now a simple recurrence relation with the initial value (2.33), which can be analytically solved: (2.34) Substituting eq.(2.34) into eq.(2.29), we obtain (2.35) Therefore, the electric part of the Hamiltonian in the rotated basis becomes which is expressed only by the fermion fields.Note that transformation of the mass term is trivial: Hm = H m .Consequently, the Hamiltonian in the rotated basis becomes with ) ) Thus far, the outgoing flux R i (N ) is not restricted, although we imposed the no incoming flux condition L i (0) = 0.If we impose no outgoing flux condition, R i (N ) = 0, we may add a penalty term, into the Hamiltonian, H + H λ .We will consider this setup in the numerical calculations.

Observables
We are interested in local observables, correlation functions, and thermodynamic quantities.
Here, we summarize these observables represented on a lattice.

Local observables
The two-component fermions ψ = (ψ L , ψ S ) in the Dirac representation correspond to and gamma matrices are γ 0 = σ z , γ 1 = iσ y , γ 5 = γ 0 γ 1 = σ x (See Appendix A for the detailed correspondence in the case of the free theory).The physical coordinate n runs from 1 to L := N/2.Bilinear local operators, the quark-number density, the current, scalar, and pseudo-scalar operators in the rotated basis, are respectively.The expectation value of an operator O is defined as ⟨O⟩ = ⟨Ω| O |Ω⟩, where |Ω⟩ is the ground state wave function.

Correlation and quark distribution functions
Another observable in the rotated basis is a two-point function, where s, s ′ ∈ {L, S}.We note that eq.(2.48) is gauge invariant.This is because χ ) in the original basis, which is manifestly gauge invariant.We define the Wigner transform of s S < s,s (n, n ′ ) by Because there is a boundary, the possible momentum is restricted, unlike the infinite system.For numerical calculations, we will use the center position, x = L/(2w), which allows the maximum number of points for the momentum.We have introduced the floor function to restrict wx + 1/2 ± ℓ/2 to an integer.The 1/2 in the argument is a shift term introduced to ensure the distance between two points changes by one when ℓ is increased by one.One could use the ceiling function instead of the floor function, whose difference disappears in the large volume limit L → ∞.Using eq.(2.49), we define one particle distribution function n(p) by The vacuum contribution is subtracted so that the number density vanishes at the vacuum.
In the case of the free theory at zero temperature and finite density in the continuum a → 0 and large volume L → ∞ limits, it corresponds to where θ(z) is the Heaviside step function, and E p = p 2 + m 2 is the energy, µ is the quark chemical potential.There is the Fermi surface at p = µ 2 − m 2 .If one turns on an attractive interaction, the Fermi surface will disappear, and the distribution function n(p) will become a smooth function.We will see this is the case for both two and three colors in section 4.

Thermodynamic quantities
We are interested in a finite-density system at zero temperature.To this end, we minimize instead of the Hamiltonian, where N B = N q /N c is the baryon number operator, and N q is the quark number operator, Here, we subtract the constant LN c N f to ensure the baryon number vanishes at the vacuum.
Similarly, the Hamiltonian is used with the vacuum energy implicitly subtracted.µ B is the baryon chemical potential that relates the quark chemical potential µ through µ B = N c µ.
The pressure P is given by where we defined the baryon number n B and energy density ε as with V = L/w being the physical volume.

Mapping to spin system
For numerical calculations, it is more convenient to use spin degrees of freedom than fermionic ones.Using the Wigner-Jordan transformation [35], we can map our fermionic system to a (N N f N c )-sites spin system: where σ ± l = (σ x l ± iσ y l )/2 with the Pauli matrices σ i (i = x, y, z), and we introduced the mapping of sites from the original system (n, f, c) to the spin system, ℓ(n, f, c) Using the spin degrees of freedom, the hopping term (2.39) is expressed as where Σ f,c (n) is defined as Similarly, the mass term is For the electric part of the Hamiltonian (2.38), we need to express the color-charge densities by spins.For this purpose, we introduce the following basis of generators: for i ̸ = j, and Let us now map Q (ij) (n) to the spin degree of freedom: for i > j, and for the diagonal parts, where ) ) ) (2.72) Similarly, the two-point function is where we defined (2.74) 3 Matrix Product state and Density matrix renormalization group We utilize the MPS under open boundary conditions as a variational ansatz, which is given explicitly as where |i 1 , i 2 , • • • i N ⟩ form the 2 N -dimensional Hilbert space of the N -site spin chain, and ,N ] are D × D complex matrices [D-dimensional complex vectors], respectively.α i are called bond indices, and D is the bond dimension.MPS is graphically expressed e.g., for N = 8 as where shapes and lines represent tensors and their legs, respectively.The connected lines imply contractions of connected indices of tensors.Similarly, we express the spin Hamiltonian H in the matrix product operator (MPO) form: Graphically, it is expressed as Our purpose is to find the wave function (3.1) that minimizes the expectation value of the Hamiltonian.To achieve this, we adopt the MPO version of the DMRG (two-site update) [33,34].The procedure of the algorithm is as follows: Using the singular value decomposition, we first prepare an initial MPS that is expressed into an orthogonal form (canonical form), , where the orthogonality center is set to the site q U y g 9 K 9 9 c r h W o 1 u l + W 4 B b c h r t 4 h z y E K j z F 8 9 r A + q / g E 3 y G L w s / 0 o v p 5 X R 2 H J q a i 5 j r 8 F 9 L 3 / w H R S 8 q 1 Q = = < / l a t e x i t > i 2 . (3.7) To achieve this, we consider an effective Hamiltonian obtained by contracting the remaining tensor with physical indices i 3 , i 4 , • • • i N and bond index α 2 .Due to the orthogonal property (partial isometry), MPS tensors such as can be understood as the basis transformation from the physical indices (i, j, k) to the virtual bond indices α.Thus, we obtain the effective Hamiltonian by transforming H into the (i 1 , i 2 , α 2 ) basis as explicitly given as Since the Hamiltonian is given by MPO, the contractions of tensors can be done sequentially and efficiently by defining and computing R j tensors as We note that the R j tensors can be recycled during the sweep of DMRG.In order to update the tensor (3.7), our task is now to find a new bond tensor B 12 , which is defined in the (i 1 , i 2 , α 2 ) basis and minimizes the effective Hamiltonian (3.10).If abstract space spanned by α 2 well approximates the subspace on the (i 3 , • • • , i N ) space, in which the full eigenvector has support, the full Hamiltonian can be minimized by optimizing the bond tensors.The DMRG algorithm tries to improve this abstract space by sequentially updating the bond tensors.
After obtaining the optimized bond tensor B 12 , which is done by e.g., the Lanczos algorithms, we need to restore the bond tensor to the MPS form.This can be done by using the singular value decomposition: which is used as the update of eq.(3.7).Here, we keep only the m largest singular values of B 12 and corresponding columns (rows) of U 1 (V 2 ).Alternatively, we can retrieve the MPS tensor by diagonalizing the density matrix: < l a t e x i t s h a 1 _ b a s e 6 4 = " U l c P 0 9 q w Y 9 < l a t e x i t s h a 1 _ b a s e 6 4 = " C v A P y 3 c w U B t g where we truncate U 1 according to the eigenvalue (S 12 ) 2 , and M 2 is given as U † 1 •B 12 .Although these methods are formally equivalent, the latter has several advantages.One of them is the noise term perturbation [39], which is used to accelerate the convergence of DMRG.The MPS tensor U 1 obtained from these methods is left orthogonal, i.e., = . (3.14) We represent the left-orthogonal tensors by solid square shapes.
Next, we move on to the next site and update the MPS tensor with indices (α 1 , i 2 , i 3 , α 3 ).We define L 1 tensor from U 1 and the i 1 component of the MPO tensor of H as Using this, H is transformed as  which is the effective Hamiltonian in the (α 1 , i 2 , i 3 , α 3 ) basis.By minimizing the transformed Hamiltonian with the Lanczos algorithms, we obtain the optimized bond tensor B 23 .Using the truncated singular value decomposition (or diagonalization of the density matrix), the bond tensor is restored to the MPS form: < l a t e x i t s h a 1 _ b a s e 6 4 = " / V 9 E S b 5 O o z q t G j w i B w 1 5 q m Y 0 s y V 5 3 i q h f g w n P z P Y d S C r u h K S V J t B e 0 n k y 8 j P O o H k V a S / Y / X V l X 9 y j S e j L T Q p 8 2 P i N u 3 F u p M 7 Z S a a 5 P z z S Z n l f q M z 9 Z j 5 q 2 J F O q E v K E s z Q n L 0 S x M 1 R 6 r h 5 2 d n g i q U 5 L 9 N l P r a f y E T E 5 q U y g 9 K 9 9 c r h W o 1 u l + W 4 B b c h r t 4 h z y E K j z F 8 9 r A + q / g E 3 y G L w s / 0 o v p 5 X R 2 H J q a i 5 j r 8 F 9 L 3 / w H R S 8 q 1 Q = = < / l a t e x i t > i 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " y x h b Y 8 V w Y t / w E P / g C A e K 3 j x 4 J v J p o 0 N 3 e 0 s u / v e m / e 9 N / P e M K Z n c z 8 g 5 H g q d W X a u H p t 5 n r 6 x u z N W 3 O Z 7 P y 2 7 3 a l x W q W a 7 t y < l a t e x i t s h a 1 _ b a s e 6 4 = " e r 2 H o 7 x 3 7 h L N 6 y I / Y V 5 y E H v 6 3 3 6 7 z + h r J L Z A 5 p z g 7 V l 7 j K I e o D z N j D z K P I 0 2 s 6 U e K 5 2 g x X u w A X n F v v G a w l s J t E a U u h n i c / n n g e 1 d F 7 F J I X H 3 8 6 s h 7 t U U h + P N P G m A 4 + Y 2 6 y t 4 o q t h N p Q a d n l k y u q + j M T 6 9 < l a t e x i t s h a 1 _ b a s e 6 4 = " X X y 9 e q l U 9 9 f J t f 7 b h y p / 1 x y H q 0 R y H 5 y U w b Y z r 4 j L j x 3 i q q 2 E 6 l B Z 2 e a T K 9 r q I z P 7 k e d < l a t e x i t s h a 1 _ b a s e 6 4 = " e r 2 H o 7 x 3 7 h L N 6 y I / Y V 5 y E H v 6 3 3 6 7 z + h r J L Z A 5 p z g 7 V l 7 j K I e o D z N j D z K P I 0 2 s 6 U e K 5 2 g x X u w A X n F v v G a w l s J t E a U u h n i c / n n g e 1 d F 7 F J I X H 3 8 6 s h 7 t U U h + P N P G m A 4 + Y 2 6 y t 4 o q t h N p Q a d n l k y u q + j M T 6 9 which is used as the update of the MPS tensor.Similarly, by defining L 2 tensor from L 1 and U 2 as we next optimize the bond tensor B 34 defined on the (α 2 , i 3 , i 4 , α 4 ) basis.We continue this procedure until we reach the end of the chain.This is half of a single DMRG sweep.Next, starting from the bond located between i = N −1 and N , we continue the procedure in reverse order back to the first bond.After reaching the first bond, one full sweep of DMRG is finished.We iterate the sweep of DMRG with gradually increasing the accuracy of truncation of the singular value decomposition until the desired accuracy of the groundstate energy is reached.

Numerical results
In this section, we show the numerical results using DMRG with the ITensor Library [36].

SU(2)
DMRG was performed using iTensor [36] with maximum bond dimension 200, truncation error cutoff 10 −8 , and noise strength 10 −8 for the final DMRG sweeps.We need typically a few thousand DMRG sweeps for convergence.

Baryon size
First, let us look at the single-baryon state to see the typical size of baryons.Concretely, we prepare the vacuum state by minimizing the Hamiltonian, by using the quantum-number conserving DMRG algorithm.We employ parameters, N = 150, J = 1/20, and w = 5.The physical spatial volume is V = N/(2w) = 15, which is sufficiently large to contain the single baryon, as will be seen numerically below.Upon preparation of the vacuum state, we apply a baryon operator, ψ † 1 (x)ψ † 2 (x) at the center position to create the N B =1 state.Then, the DMRG optimization is performed to minimize its energy with fixed baryon number N B = 1.This process is stopped at a finite number of sweeps, around 2000.The baryon remains localized at the center position in the resultant state that is not a true minimized energy state; however, the energy is low enough to see the distribution of the single baryon, where the energies are 1.59, and 2.57 for m = 0.5 and 1.0, respectively.Note that the baryon is localized at the boundary in the true minimum state.The baryon size can be inferred from the expectation value of the baryon number density, shown in figure 1.We can see that the baryon is localized at the center of the position.In addition, there are localized modes at the boundaries.This is due to the open boundary condition.In calculations at finite density, particularly at low densities, boundary effects are expected to be unavoidable; thus, the volume should be as large as possible with keeping densities.We fit the distribution by a Gaussian function A exp −(x − x 0 ) 2 /(2r 2 ) with parameters A, x 0 , and r, using the data points 3.4 < x < 11.8.The fitting results are also shown in figure 1.From these fitting results, the baryon sizes for m = 0.5 and 1.0 are estimated to be r = 1.07 and 0.88, respectively.
In figure 2, we show the quark-distribution function defined in eq.(2.50) for the single baryon state.Both for m = 0.5 and m = 1.0, the distributions rapidly decay as p increases.The distribution function decays of order 1 = g 0 , reflecting the fact that the baryon size is of order 1/g 0 .The decay rate of the distribution function does not seem to depend significantly on mass.Because the quarks are confined within the baryons, no Fermi surface or Fermi sea of quarks is formed in the single baryon state.

Thermodynamic quantities
Next, we examine the thermodynamic quantities such as the pressure and the number density.At finite density, we minimize H − µ B N B (2.52) using the DMRG algorithm without quantum number conservation.We employ parameters N = 160, J = 1/8, and w = 2, where the physical volume V = 40 is much larger than the typical size of baryons ∼ 1.The left panel of figure 3 shows the pressure as a function of µ B .Due to a finite volume, the pressure is continuous but not smooth.This fact can be seen from the baryon number density, which is the derivative of the pressure, shown in the right panel of figure 3. The baryon number density shows a step-like behavior.Due to the confining energy, the pressure and baryon number density start to increase at a higher point than the threshold of the free theory µ B = N c m.The threshold values of the chemical potential are 1.26 and 2.27 for m = 0.5 and m = 1.0, respectively.The baryon number density first rises sharply and then increases linearly.The behavior can be understood as follows: At high densities, g 0 /µ ≪ 1, and m/µ ≪ 1, the contributions from the interaction and mass become negligible.As a result, the system is approximated by free-quark fields, and then the baryon number density behaves as n B = N c µ/π = µ B /π, which is a linear function of µ B .The behavior at low density is more nontrivial, and we discuss the physical meaning along with other observables below.
When we compare physical quantities with different masses, it is better to consider them as functions of the baryon number density (or energy density) rather than the chemical potential.For the energy, this can be determined unambiguously; however, subtleties arise when considering the pressure.This is caused by the ambiguity of the chemical potential; the baryon number density n B is not a continuous function of µ B , and therefore, the inverse function does not uniquely exist.Given that the energy is uniquely determined as a function of the baryon number, the chemical potential can be defined as the change in the energy when one more baryon is added: where E(N B ) = ⟨H⟩ is the total energy for a given baryon number.The choice between a forward or backward difference causes ambiguity.Note that when we compute a physical quantity as a function of the chemical potential, the number and energy densities remain constant within the interval µ − B (n B ) < µ B < µ + B (n B ).Alternatively, we can introduce µ B (n B ) as the average of µ (4.4) In the large volume limit V → ∞, all values of µ ± B (n B ) and µ(n B ) converge to the same value.Since the averaged chemical potential µ B (n B ) has an improved volume dependence compared to µ ± B (n B ) at a large volume, we define the pressure as a function of n B by using µ B : The numerical results are presented in the left panel in figure 4. We also plot the energy density as a function of n B in the right panel in figure 4. For comparison, we plot the pressure and energy density of the free theory (See appendix A for details of the free theory).The To see the difference between QCD 2 and the free theory in more detail, let us focus on behaviors of the energy per quark and averaged chemical potential.The left and right panels in figure 5 show the energy per unit quark number ε/(N c n B ) and quark chemical potential µ = µ B /N c , respectively.These quantities will help us understand whether the system behaves like quarks or like baryons, by comparing their behavior in the theory.In regions of a large baryon number density, ε/(N c n B ) behaves like free quarks.On the other hand, at low densities, the energy per quark and averaged chemical potential are higher than those of the free theory due to the effect of confinement energy.The change in density is more gradual than that in the free theory, and the behavior seems to change to that of the free theory around n B ≈ 0.2.
One of the important characteristics of two color QCD 2 is that the baryons are bosons due to the even number of colors.Therefore, if we assume interactions are negligible at low density, the baryons will degenerate to the lowest energy state.Consequently, ε/(N c n B ) and µ are expected to be constant at low densities, which is consistent with the behaviors in Fig. 5.However, we must be careful with this consideration, since interactions are generally not negligible in (1 + 1)-dimensional systems.
Lastly, in this subsection, let us look at the averaged chiral condensate and sound velocity.The averaged chiral condensate is given as We need a renormalization of Σ because it diverges in the continuum limit.We are interested in the change of Σ at finite density, so we introduce Here, Σ vac is the unrenormalized chiral condensate in the vacuum.The left panel of figure 6 shows ∆Σ as a function of n B together with the results of the free theory.Both cases m = 1 and m = 0.5 behave similarly to the free theory, with the difference being more significant for the m = 0.5 case.This may be because the lighter masses are more sensitive to chiral symmetry breaking.However, we do not discuss the spontaneous breaking of chiral symmetry and its restoration for several reasons.Firstly, the masses we employ are relatively heavy.Secondly, in the Hamiltonian formalism, staggered fermions do not possess continuous chiral symmetry.Lastly, the open boundary conditions explicitly break chiral symmetry.
In the continuum theory, the sound velocity c s is given by which represents the response of the pressure to the change of energy.At zero temperature, we can express it by using the chain rule as We evaluate the sound velocity by replacing the derivative with the central difference, which can be expressed by µ(N B ) as The right panel of figure 6 shows the squared sound velocity as a function of n B together with the results of the free theory.The overall behavior is the same as before, behaving like a free theory.At low densities, the sound velocity is suppressed compared to the free theory.This is consistent with the free baryon picture, though interactions might play a significant role.
In the context of neutron star physics, the possibility of a peak in the sound velocity has been discussed [40,41], and may exist in two color QCD [42,43].This is due to the fact that in a (d + 1)-dimensional system, c 2 s asymptotically approaches 1/d at high densities.Consequently, there may exist an intermediate density region where c 2 s > 1/d.This is not the case in our d = 1 situation.Assuming that the sound velocity does not exceed the speed of light, there is no peak satisfying c 2 s > 1/d

Size dependence
While the purpose of this paper is to qualitatively understand QCD 2 at finite density, it is worthwhile to look at the lattice spacing and volume dependence.For this purpose, we focus on the energy density at a fixed baryon number density with w = 2. Figure 7 shows the volume dependence of the energy density for 2wV = N = 80, 120, 160, 200, 240.The volume dependence can be approximated by a linear function, with slope and intercept for m = 0.5 being −0.24 and 0.72, and for m = 1.0 being −0.21 and 1.04.Therefore, the energy densities in the infinite volume limit are 0.72 for m = 0.5 and 1.04 for m = 1.0, respectively.The results at N = 160 corresponding to V = 40 coincide with those in the infinite volume limit with an error of less than 1%.We also checked that the energy densities for w = 2 (a = 0.25) and w = 4 (a = 0.125) coincide with an accuracy of approximately 1% in the range 0 ≤ n B < 0.5.

Inhomogeneous phase in a finite volume
In Gross-Neveu, chiral Gross-Neveu models, and massless QCD in (1 + 1) dimensions, inhomogeneous phases have been considered to exist within the mean-field approximation [46][47][48].It is an interesting question whether such an inhomogeneous phase can be realized in QCD 2 in the presence of quark mass.In this paper, we discuss the realization of the inhomogeneous phase in our method.We, however, need to be careful about the possibility of an inhomogeneous phase.In (1+1) dimensions, the spontaneous breaking of continuous symmetry  which wavenumber the modulation occurs.Figure 9 shows the wavenumber with the largest amplitude for QCD There are two ways of looking at this periodic structure.One is based on the baryonic picture.
Assuming repulsive interactions between baryons, if a box of length V is packed with equally spaced baryons, the spacing is V /N B = 1/n B , and thus, the wave number is k = 2πn B .The other is based on the quark picture.If there is no interaction, the quarks form a Fermi surface at p = p F = πn B .If one turns on an attractive interaction between quarks, the Fermi surface becomes unstable due to the Peierls instability and condensation of particle and hole pair occurs [52].In this condensation, the relative momentum of the particle and hole is zero, but the total momentum is 2p F .Therefore, the wave number is k = 2p F = 2πn B .The same modulation is obtained regardless of whether the quark or baryon picture holds.This observation implies that the modulation is independent of the number of colors.As will be discussed in section 4.3, the same behavior will be found in the SU(3) case.We shall see that the amplitude of this modulation is significantly larger than the free case at any density.Figure 10 shows the largest amplitude in the Fourier decomposition of the baryon number density.The figure shows that the amplitude in QCD 2 is roughly twice that of the free theory.In the case of QCD 2 , there is a point where the amplitude behaves discontinuously, but this is presumably due to the finite volume and not a phase transition.From this numerical calculation for QCD 2 , it is challenging to ascertain whether the amplitude with modulation persists in larger volumes.Instead, let us look at the distribution function of quarks and argue that the occurrence of condensation is plausible, given in the absence of Fermi surfaces.

Quark distribution function
It is interesting to see how quarks are distributed in the inhomogeneous phases.We defined the gauge-invariant quark distribution function n(p) in eq.(2.50).In order to observe the transition from hadronic matter to quark matter, we calculate n(p) at N = 240 and w = 2.0 varying the baryon density (n B = 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5).We employ the larger site number compared with N = 160 used in the previous sections, to reduce the finite volume effect.This is because n(p) around p = 0 is sensitive to the volume size.The results for m = 0.5 and m = 1.0 are shown in the left and right panels of figure 11, respectively.As will be discussed below, these figures indicate that condensation is caused by the instability of the Fermi surfaces.
First, the two figures show that the difference in mass does not seem to affect the qualitative behavior of the distribution function.Next, let us look at the density dependence of the distribution function.At low densities, the structure has a peak at p = 0 and decays at large p.This shape is similar to the quark distribution in the single baryon shown in figure 2. This would mean that the degrees of freedom of the system are baryonic.As the baryon number density increases, the maximum of n(p) at p = 0 increases.After the maximum reaches 1, n(p) forms the Fermi sea.This behavior is consistent with the argument based on the quarkyonic picture in (3 + 1) dimensions discussed in ref. [41].Numerical results show a depression at p = 0 when the density is high.We expect that the depression near p = 0 is a lattice artifact because such depression also appears in the free theory.Both panels of the figure 11 show that the Fermi sea begins to form near n B = 0.2.We expect that the behavior of the system changes from a baryonic matter to a quark matter.This is consistent with the behavior of the energy per unit baryon number and the averaged chemical potential shown in figure 5, which changes to the free-quark behavior around n B = 0.2.
We note that such a behavior can also be observed in ultracold atomic gases, which exhibit the crossover from Bardeen-Cooper-Schrieffer (BCS) superfluids to molecular Bose-Einstein condensates (BEC), tuning the interaction using the Feshbach resonance.Although the distribution function of the ultracold atom fermion does not reach the maximum value, 1, at the origin in the BEC (bosonic) region, it forms the Fermi surface as it changes to the BCS (fermionic) region [53,54].
In BCS-type condensation, the half-density point coincides with that of the free theory because the degrees of freedom near the Fermi surface contribute symmetrically to the condensation.Our results in figure 11 show that this is realized for n B > 0.2.This is again consistent with the transition from baryon to quark degrees of freedom around n B = 0.2.Precisely speaking, it is an inhomogeneous phase with density waves, so it is different from BCS-type condensation, but the mechanism of condensation due to the instability of the Fermi surface itself is the same.

SU(3)
In the previous sections, we have numerically studied QCD 2 for N c = 2 at finite density.In this section, we generalize it to the case of N c = 3.The numerical calculation for N c = 3 is more computationally demanding than that for N c = 2.This is due to the increased color degrees of freedom as well as the enhanced nonlocality of the Hamiltonian represented by the spin system.Therefore, we present simulation results for a relatively smaller size with parameters N = 48, w = 2, J = 1/8, which corresponds to the physical volume size V = 12 in g 0 = 1 unit.DMRG was performed using iTensor [36] with maximum bond dimension 500, truncation error cutoff 10 −8 , and noise strength 10 −8 for the final DMRG sweeps.We need a few thousand DMRG sweeps for convergence.
We show graphs corresponding to figures 3-11 (except figure 7) in figures 12-18 for m = 1.As in the case of N c = 2, the pressure and baryon number density start to increase from µ B > N c m due to the confining energy (see figure 12).The critical value of the baryon chemical potential is given as µ B = 3.58 (> 3m).The overall behavior is similar to the case of N c = 2, even though the properties should be different since baryons are fermions and not bosons.Thermodynamic quantities behave like free quarks at high densities.As in the case of N c = 2, the inhomogeneous phases are realized, shown in figure 16, and the wavenumber of oscillations with the maximum amplitude is given by k = 2πn B (see figure 17).The physical interpretation is the same as for N c = 2, which is independent of the number of colors.
It is not easy to see the transition from baryons to quarks from behaviors of the energy per quark and the averaged chemical potential in figure 14.One of the reasons is that the volume is not large, so the number of points is not enough.On the other hand, the behavior of the distribution function shows the crossover transition from baryons to quarks around n B = 0.3 in figure 18.
Finally, let us consider the low-density behavior from the point of view of large N c counting.When comparing free quarks with free baryons, we find that both quarks and baryons have the same Fermi momentum p F = πn B .This is a unique property in (1 + 1)-dimensional theory.The pressure of baryons and quarks for a given baryon number density n B are respectively.Here, m B is the baryon mass of order N c , and m q is the (constituent) quark mass of order N 0 c .Since m B ∼ N c m q , the pressures for the free baryons and the free quarks are different by an order N 2 c .We plot the pressures of the free baryons in figure 13 for comparison.It can be seen that the pressure of free baryons is strongly suppressed compared to that of free quarks.The result from the numerical calculations (red dots) is closer to that for the free quarks rather than for the free baryons.It implies that the contribution from the interactions of order N c plays an essential role.Note that the fact that the interaction between baryons is of order N c is consistent with the large-N c counting [55].to QCD 2 .Not only is this computational method efficient, but also it is free of sign problems.We have calculated thermodynamic quantities, such as the pressure, baryon number density, and energy density.We have also calculated the quark distribution function from the Fourier transform of the Wigner function.Our study may pave the way toward ab-initio study of cold-dense QCD on the basis of tensor network methods.
We have several future directions, which make the analysis more realistic by introducing gauge fields, but may be reachable with current computational resources (although those are still far from real QCD).First, if we employ periodic boundary conditions in QCD 2 , one link corresponding to the spatial Wilson loop is no longer integrated out by gauge transformations.This would be the most simple setup to simulate a nontrivial nonabelian gauge field, although it is not dynamical.Furthermore, we may be able to study the effects of dynamical gauge fields by considering quasi-one-dimensional lattice such as a two-leg ladder geometry, and formulating gauge fields on the basis of the recently developed Hamiltonian Yang-Mills theory [56][57][58][59][60].

k a 1 e 5 b t 9 7 7 3
P u / d L 5 3 p 2 c I P G D u a S 8 0 v p C 8 t L l 3 O X L l 6 b X k l m 1 v d 8 t 2 + s n j D c m 1 X 7 Z i G z 2 0 h e S M Q g c 1 3 P M U N x 7 T 5 t r n / R P d v D 7 j y h S s 3 g 6 H H 9 x y j H H s s Q / y b c w 9 c Z Z + B e a + X l j Y 0 / i a y D / w B 6 p 1 T s y A P o w C M a s c A Z e K T o u V h j f v D i 9 f H G 4 3 o x v M P e s V 8 4 i 7 f s i H 3 F e c j B b + v 9 O q + / o e w S m U O a s 0 P 1 J a 5 y i P o A M / Y w 8 y j y 9 J p O l H i u N s P V L s A F 5 9 Z 7 B m s J 7 C Z R 2 l K o 5 8 m P J 5 5 H d f Q e h e T F x 5 + O r E d 7 F J I f z 7 Q x p o P

7 l m 2 3
/ v e + 7 x 3 v 3 S W 7 4 g g Z O x g L n N u 3 j i / s H g h e / H S 5 a X l X P 7 K Z u D 1 l c 0 b t u d 4 a t s y A + 4 I y R u h C B 2 + 7 S t u u p b D t 6 y 9 R 7 p / a 8 B V I D y 5 E Q 5 9 v u u a X S k 6 w j Z D l F q 5 6 H s E p l 9 m r N L 9 S W u c o T 6 A D P 2 M P M o 9 v S a T p R k r j b D 1 c 7 A h a f W e w o r K e w G U d p S q B f I T y a e x 3 X 0 H k X k J c c f j 6 x H e x S R n 8 y 0 M a a D z 5 i b 7 K 2 i i u 1 U W t D p m S X T 6 y o 6 8 9 P r U S c l m d I r o Y 4 4 e b 4 m o X 4 O S 5 9 Z 7 B a g a 7 q S l l C d R L 2 k 8 n D p M 6 a o 8 i 7 a X H n 4 6 s p / c o 0 n 4 6 0 8 a Y D j 4 j b r y 3 Q l d s Z 9 J M n 5 5 p M r u u 0 4 H l t Y P 7 3 c A z f 4 P v 8 S L u v F b S l s W t m L m b u w n 9 D e / w P h i I x / g = = < / l a t e x i t > (S 12 ) 2 , n m X 7 v e + 9 z 3 v 3 S 2 f 5 j g h C x g 7 n M h f m j Y s L i 5 e y l 6 9 c X V r O 5 a 9 t B V 5 f 2 b x h e 4 6 n d i w z 4 I 6 Q v B G K 0 O E 7 v u K m a z l 8 2 9 p / r P u 3 B 1 w F w p O b 4 d D n e 6 7 Z l a I j b D N E q Z U r N S 3 e F T L i z 6 S p l D

4 ,
7 B b b i L d 8 g D q M I T P K 8 N r P 8 K P s F n + D L / w 1 g w l o z c O D Q z F z P X 4 b 9 m 3 P w H 3 i k q w A = = < / l a t e x i t > R 4 8 p m n T F p o a Z j e T Z O j u 7 D q 7 S Y 1 L / k B P v X n w p O B B v H t V 8 O I f 8 N A / I I j H C l 4 8 + G a y a W N D d z v L 7 r 7 3 5 n 3 v z b w 3 j O U 7 P A g J O Z r J X J o 1 L l + Z u 5 q 9 N n / 9 x k I u f 3 M z 8 H r S Z n X b c z y 5 b d G A O V y w e s h D h 2 3 7 k l H X c t i W t b e q 5 r f 6 T A b c E x v h w G e 7 L u 0

k a 1 e 5 b t 9 7 7 3
P u / d L 5 3 p 2 c I P G D u a S 8 0 v p C 8 t L l 3 O X L l 6 b X k l m 1 v d 8 t 2 + s n j D c m 1 X 7 Z i G z 2 0 h e S M Q g c 1 3 P M U N x 7 T 5 t r n / R P d v D 7 j y h S s 3 g 6 H H 9 x y j H H s s Q / y b c w 9 c Z Z + B e a + X l j Y 0 / i a y D / w B 6 p 1 T s y A P o w C M a s c A Z e K T o u V h j f v D i 9 f H G 4 3 o x v M P e s V 8 4 i 7 f s i H 3 F e c j B b + v 9 O q + / o e w S m U O a s 0 P 1 J a 5 y i P o A M / Y w 8 y j y 9 J p O l H i u N s P V L s A F 5 9 Z 7 B m s J 7 C Z R 2 l K o 5 8 m P J 5 5 H d f Q e h e T F x 5 + O r E d 7 F J I f z 7 Q x p o P P m t a F S 7 Z 9 l + 7 3 v v 8 9 7 9 0 l m + I 4 K Q s Y O 5 z L l 5 4 / z C 4 o X s x U u X l 5 Z z + S u b g d d X N m / Y n u O p b c s M u C M k b 4 Q i d P i 2 r 7 j p W g 7 f s v Y e 6 / 6 t A V e B 8 O R G O P T 5 r m t 2 p e g I 2 w b c w d c d Z + B + a / n l t f U / q a y L / x B 6 x 1 T i y E P o w E M a s c A Z + K T o u d h j f v D i 9 e H 6 o 3 o p u s X e s V 8 4 i 7 f s g H 3 F e c j B b / v 9 G q + / o e w S m X 2 a s 0 v 1 J a 5 y h P o A M / Y w 8 y j 2 9 J p O l G S u N s P

k e b 6 9 3 ≃
E y S y X m l P v P j 9 a h q S z y l K i F P O V N z 8 l I U O 0 c l 5 + p i Z w e n k u q 0 R J / D x H o q H x G R I 2 2 M x l u l d P 4 O m R S 2 l 4 u l t e L q x n K + X I 7 u l x m 4 C / f h I d 4 h j 6 A M T / G 8 1 j D / K / g E n + H L 9 A 9 j 1 l g w b g 9 d U 1 M R s w D / D e P e P + o F L Q E = < / l a t e x i t > ↵ < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 l s B N x 0 p N m d R 4 b w a e s Q W 3 f Q s

k a 1 e 5 b t 9 7 7 3
P u / d L 5 3 p 2 c I P G D u a S 8 0 v p C 8 t L l 3 O X L l 6 b X k l m 1 v d 8 t 2 + s n j D c m 1 X 7 Z i G z 2 0 h e S M Q g c 1 3 P M U N x 7 T 5 t r n / R P d v D 7 j y h S s 3 g 6 H H 9 x y j

9 7 9
0 t n e r b w A 8 Y O Z j L n Z r P n L 8 z N 5 y 5 e u n x l I V 9 Y 3 P T d v r J 4 w 3 J t V 2 2 b h s 9 t I X k j E I H N t z 3 F D c e 0 + Z a 5 9 0 j 3 b w 2 4 8 o U r N 4 K h x 3 c d o y

k a 1 e 5 b t 9 7 7 3
P u / d L 5 3 p 2 c I P G D u a S 8 0 v p C 8 t L l 3 O X L l 6 b X k l m 1 v d 8 t 2 + s n j D c m 1 X 7 Z i G z 2 0 h e S M Q g c 1 3 P M U N x 7 T 5 t r n / R P d v D 7 j y h S s 3 g 6 H H 9 x y j

Figure 1 .Figure 2 .
Figure 1.Spatial distribution of the baryon number N B = 1 for m = 0.5 (left panel) and m = 1.0 (right panel) with N = 150, w = 5.0.The black dots are numerical results, and the red lines are Gaussian fits of the black points.

Figure 3 .
Figure 3. Pressure P (left) and baryon number density n B (right) for SU(2) as functions of µ B for m = 0.5 and 1.0 with w = 2, and N = 160.The pressure and number density begin to increase at µ B = 1.26 for m = 0.5 and at µ B = 2.27 for m = 1.0, respectively.

Figure 6 .
Figure 6.Chiral condensate subtracting the contribution from the vacuum ∆Σ = Σ − Σ vac (left) and squared sound velocity (right) for SU(2) as functions of n B for m = 0.5 and m = 1.0 with w = 2.0 and N = 160.

Figure 7 .
Figure 7. Size dependence of the energy density fixed w = 2 and n B = 0.4 for 2wV = N = 80, 120, 160, 200, 240.The dashed vertical line corresponds to the point of N = 160.The energy densities in the infinite volume limit, obtained by linear fitting, are 0.72 for m = 0.5 and 1.04 for m = 1.0, respectively.

Figure 9 .
Figure 9. Wave number of n B (x) with the largest amplitude as a function of n B for m = 0.5 (left) and m = 1.0 (right) with w = 2.0 and N = 160.

Figure 10 .
Figure 10.Largest amplitude in the Fourier expansion of n B (x) as a function of n B for m = 0.5 (left) and m = 1.0 (right) with w = 2.0 and N = 160.

Figure 16 .
Figure 16.Spatial dependence of the baryon number density n B (x) and the chiral condensate subtracting the contribution from the vacuum ∆Σ(x) = Σ(x) − Σ vac (x) with m = 1.0, w = 2 and N = 48.

Figure 17 .
Figure 17.Wave number of n B (x) with the largest amplitude (left) and amplitude (right) as functions of n B with m = 1.0, w = 2 and N = 48.

50 Figure 18 .
Figure 18.Quark distribution function n(p) for m = 1.0 with N = 48, w = 2.0, and n B = 0.25, 0.33, 0.42, 0.50.Colored dots represent the intersections of n(p) and the Fermi momentum of the free theory p F = πn B .The solid black line shows n(p) = 0.5.

(A. 46 )ψiγ 5
Here, we used eqs.(A.18) and (A.19) for the vacuum part.w is the vacuum contribution that can be eliminated through the renormalization of the charge.The remaining terms are the oscillation terms due to the open boundary conditions, which will vanish in the large volume limit L → ∞.On the other hand, the spatial component of the current vanishesψγ 1 ψ(n) = S < S,L (n, n) + S < L,S (n, n) = 0. (A.47)Similarly, ψψ(n) and ψiγ 5 ψ(n) can be evaluated asψψ(n) = S < L,L (n, n) − S < S,S (n, n) ψ(n) = iS < S,L (n, n) − iS < L,S (n, n) p) sin(p(4n − 1)) ,(A.49)respectively.Here, we decompose them into static and oscillating parts.Let us finally provide an expression for the two-point function that appeared in the distribution function: S < L,L (n, n ′ ) + S < S,S (n, n ′ cos(p) cos p 2n + 2n ′ − 1 + wδ n,n ′ .(A.50) 1.