{'fraction_non_alphanumeric': 0.033933862491078576, 'fraction_numerical': 0.01016501935679218, 'mean_word_length': 4.468070009460738, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Baryons and Low-Density Baryonic Matter in 1+1 Dimensional Large N_c QCD with Heavy Quarks

This paper studies baryons and baryonic matter in the combined large N_c and heavy quark mass limits of QCD in 1+1 dimension. In this non-relativistic limit, baryons are composed of N_c quarks that interact, at leading order in N_c, through a color Coulomb potential. Using variational techniques, very accurate calculations of single baryon masses and interaction energies of low-density baryon crystal are performed. These results are used to cross-check a general numerical approach applicable for arbitrary quark masses and baryon densities recently proposed by Bringoltz, which is based on a lattice in a finite box with periodic boundary conditions. The Bringoltz method differs from a previous approach of Salcedo, et al. in its treatment of a finite box effect - namely gauge configurations that wind around the box. One might expect these effects to be small for large enough boxes, in which the baryon density approaches zero to high accuracy at the edges. However, the effects of these windings appear to be quite large even in such boxes. The large mass infinite volume calculations performed here are consistent with the results of numerical calculations using the Bringoltz method. The calculation of the baryon crystal interaction energy requires the assumption that at low-densities the ground state is composed of individual baryons, each in a color-singlet state and orthogonal to each other. This assumption is plausible but ad hoc in that one can construct configurations in which the entire state is color-singlet but cannot be broken into individual color-singlet baryons. The interaction energy of low-density baryon crystals calculated with the assumption is consistent with numerical results based on Bringoltz's approach suggesting that the assumption is justified. This further supports a similar assumption that was made in 3+1 dimensions, where no alternative means of calculation exist.


I. INTRODUCTION
The problem of understanding cold nuclear matter from QCD first principles has remained quite challenging due to the intrinsic non-perturbative nature of the problem [1,2]. While many non-perturbative problems such as the finite temperature problem can be studied using lattice QCD techniques, the problem of finite density QCD is marred by the presence of the fermion sign problem [1,[3][4][5]; only low densities can be studied and that too only at large temperatures [6] as all available methods are only valid for small values of µ T . As such it is useful to find certain limiting cases for the study of cold nuclear matter, in which the problem of finite density QCD becomes somewhat more tractable.
Large N c QCD is a limit where the physics simplifies considerably [7]. Only planar diagrams contribute at leading order in N c while quark loops effects are suppressed. Despite these simplifications, the problem of large N c QCD has remained unsolved. However, the problem in 1+1 dimension is more tractable [8,9]; here the physics simplifies considerably because gluons, due to the absence of magnetic fields in one spatial dimension, do not propagate. Additionally, with appropriate choice of gauge (either axial or Coulomb), gluon-gluon interactions are absent and hence the theory looks very much like QED, in which color charges (i.e. quarks) interact via a Coulomb potential. The self-energy corrections to the quark propagator come from an infinite number of ladder diagrams, which remarkably can be summed up exactly in spite of the need to solve a bootstrap equation. Furthermore, the meson spectrum can be solved -albeit only numerically via a Bethe-Salpeter equation; a linear "Regge trajectory" (with no continuum) is observed.
Unlike the meson problem in 1+1 large N c QCD, the baryon problem is highly non-trivial. Many attempts have been made to study the problem: Ref. [10] studies the problem in the chiral limit, where baryons assume the form of a Sine-Gordon soliton. Ref. [11] invokes Eguchi-Kawai reduction [12] to argue controversially [13] that 1+1 large N c QCD is independent of baryon chemical potential. Ref. [14] formulates the large N c baryon problem for arbitrary quark masses at infinite volume. Numerical solutions of the problem require a lattice treatment in a finite-sized box. The most exhaustive work was done recently by Bringoltz in Ref. [13]. It is based on the axial gauge representation of QCD [15] and the coherent states formalism [16]. Ref. [13] solves both the single baryon and baryon crystal problems in a finite sized box for arbitrary quark masses. The principal difference between the approaches of refs. [13] and [14] is the treatment of the effects of finite box sizes. In Ref. [14] the infinite volume Hamiltonian is derived and simply used on the finite box by imposing boundary conditions. In contrast, the idea underlying Ref. [13] is that the quarks making up a baryon in a finite box interact via a color-Coulomb potential of the form |x − y + QL|, where |x − y| is the physical separation of the color charges, Q is an integer and L is the size of the box. These additional windings are due to the inability to gauge away the gluon degrees of freedom in the problem [13,15]. In effect, the treatment of Ref. [14] corresponds to including only the Q = 0 sector.
The study of these windings in Ref. [13] was motivated largely for theoretical reasons associated with QCD volume reduction. Clearly in small boxes the two treatments might be expected to be quite different. A priori, however, one might expect that the difference between calculations based on Refs. [13] and [14] for a baryon mass would be quite small in cases where the box sizes are large enough to comfortably contain a baryon since the difference between the two is a finite volume effect. However, as will be discussed in section III, the differences can be quite large, even in such circumstances. Indeed, the differences are large enough that one might worry the differences suggest a possible problem with the formulation of Ref. [13]. Unfortunately, It is difficult to verify numerically that this is the case by going to exceptionally large box sizes.
Given this situation, it is quite useful to have a regime in which one can independently compute the mass of a baryon to verify that the winding effects are as large as they appear to be for practical-sized boxes. The regime of heavy quark masses and large N c serves this purpose extremely well. In this paper, we perform a highly accurate variational calculation of single baryons using Witten's mean-field argument [17], which is similar to an analogous 3+1 dimensional calculation that was done recently in Ref. [18]. It was argued in Ref. [17] that baryons are composed of N c quarks in a color singlet state, where each quark moves in an average potential due to the remaining N c − 1 quarks. In the non-relativistic limit with quark masses taken to infinity, the average potential has a simple form because the only contributions to the interaction come from single gluon exchange. Note that the principal difference between the heavy quark system using Witten's approach and the general result of Ref. [13] is that the vacuum is trivial in the former case.
We find that the single baryon results using Witten's variational formulation agree with the results based on Ref. [13] that incle the nontrivial effects of the winding with remarkable accuracy. These results support the approach taken in Ref. [13] that include a large effect of the windings, albeit only in the restricted regime of large quark masses. Our result is also quite useful in providing a simple and physically intuitive picture for at least one regime of QCD in 1+1 dimension and, as such, nicely complements the Bringoltz method, which obtains results largely through numerical means.
The study of the heavy quark regime in 1+1 dimension is useful for another purpose: the study of baryonic matter at low density. Studies of low-density baryonic matter in the heavy quark regime of large N c QCD were previously done in 3+1 dimensions [18]. The principal attraction of this regime is that it appears to be tractable due to Witten's mean-field arguments and the non-relativistic nature of heavy quark dynamics. However, as discussed in Ref. [18] for the case of baryonic matter (as opposed to a single baryon) a plausible but ad hoc assumption was needed: that at low densities the ground state is composed of individual baryons each in a color-singlet state with the wave function of quarks for each baryon orthogonal to the wave function of quarks in any other baryon. This structure appears quite plausible, but it is ad hoc in that one can construct mean-field configurations in which the entire state is color-singlet but cannot be broken into individual color singlet baryons. In 3+1 dimensions, it was not easy to verify whether this assumption was justified. However, in 1+1 dimension, one can use the methods of Ref. [13] to numerically compute the interaction energy and compare it with the result in the mean-field approach using the assumption of no hidden-color states. We find good agreement between them, which suggests that this assumption is indeed valid. This agreement in 1+1 dimension greatly increases our confidence that the assumption made for the 3+1 dimensional case is valid.
The paper is organized as follows. In the following section, we calculate the single baryon energy and wave function in the limit of large quark masses using Witten's mean-field argument. We compare this variational solution to the numerical results based on the recent Hamiltonian approach proposed of Ref. [13] and with a previous Hamiltonian in Ref. [14] that did not include the windings. In the next section, we calculate the interaction energy per unit cell of a low density baryon crystal, again using Witten's mean-field argument. We compare the interaction energies with numerical results for the many baryon sector based on Ref. [13]. The paper concludes with a brief discussion of the results.

II. BARYONS IN THE COMBINED LARGE Nc
AND LARGE mq LIMITS

A. The Energy Functional in the Mean Field Approximation
The nonrelativistic mean-field energy functional for baryons in 3+1 dimensional large N c QCD for heavy quarks has been calculated in Ref. [18] based on Witten's classic analysis in Ref. [17]. QCD in 1+1 dimension at large N c is generally simpler than in 3+1 dimensions. In 1+1 dimension, quarks are spinless, nonrelativistic fermions and gluon-gluon couplings can be eliminated by appropriate choice of gauge (e.g. the axial gauge as in refs. [8] and [14]). Thus the leading order contribution to the energy comes from a coherent one-gluon exchange between pairs of N c quarks. Interactions involving multiple gluons are suppressed at least at a relative order of 1 Nc . Furthermore, quarks loops are suppressed at relative order of at least 1 Nc as a result of the large N c limit. However, for the case of heavy quarks, these simplifications do not help. For heavy quarks, the interactions are pairwise color-Coulomb in any event so the lack of gluongluon interactions adds no new simplifications. Similarly, spin effects are suppressed for heavy quarks. Thus, the leading-order energy functional in the dual expansion, is derived in an essentially identical manner to the 3+1 dimensional case [18]. However, the color-Coulomb potential in the context of one spatial dimension is linear in distance as opposed to dropping off as 1 r . The strength of the coupling is chosen to be g 2 4 , adopting the conventions in refs. [8], [13] and [14].
The energy of the baryon per color at leading order in where m q is the quark mass, ψ(x) is the normalized baryon wave function and λ ≡ g 2 N c . Minimizing the energy functional with respect to the wave function while imposing the normalization leads to a Schrődinger equation for the ground state baryon.
where ǫ is the Lagrange multiplier that ensures the normalization of the wave function.

B. The Variational Calculation
We choose to calculate the ground state baryon wave function ψ min (x) using a variational approach. From here on, ψ(x) will be taken to be the ground state wave function ψ min (x). To proceed, we define a scaling variable R and a scaled wave function f (y) where y is the dimensionless variable defined as y ≡ x R . Then, the baryon wave function is: The normalization of the wave function ψ(x) is independent of the scaling variable R as long as f (y) is normalized, i.e. dy|f (y)| 2 = 1.
Next we define functionals for the kinetic and potential energies, T and V: Rewriting the energy functional in terms of the new variables R, T [f ] and V[f ], we get: For simplicity, we choose a normalized functional form for f (y) to be a polynomial times a Gaussian: where n indicates the number of terms kept. Only even powers of y appear in the wave function because the ground state is expected to be symmetric. We found that truncation at n = 9 was sufficient to obtain the energy with a very high accuracy. We now have an energy functional that is a function of R and a i . The next step is to minimize the functional with respect to each of the variables, i.e. ∂ h ψ|H|ψ = 0 for h ǫ {R, a i }. Minimizing with respect to R gives which yields: To proceed we minimize ψ|H|ψ with respect to each a i for i = 1 to 9 with a 0 chosen to be 1. Including R, there are a total of 10 parameters in the variational calculation.
The binding energy of the baryon is the second term in Eq. (10) and is characterized by the quantity T which fixes the wave function.
It is sensible to present the results in units of λ i.e. g 2 N c , which is the only scale present in large N c QCD. Hence masses will be written in the units of √ λ and lengths in the units of 1 √ λ . The ground state energy per color in units of √ λ is: The variational calculation converged well and the interaction energy is estimated to be accurate to five significant figures. The parameter R is given by: the dimensionless Lagrange multiplier is C. Asymptotic Behavior of the wave function The nature of variational calculations is such that they determine energies more accurately than the wave functions. The Gaussian form of the trial wave function we chose was very convenient and provided for an extremely accurate calculation for the energy. Moreover, with a sufficient number of terms, it provides a rather accurate form for the wave functions over the bulk of the region that contributes substantially to the energy. However, it does not capture the correct behavior of the wave function at asymptotically large distances. At large x, the drop off of the true solution of Eq. (2) is slower than a Gaussian [18]. Since one principal goal of this paper is to compute the interaction energy of low density baryonic matter -which is highly sensitive to the long range part of the wave function -it is important to determine the correct behavior at asymptotic distances.
To proceed, we study Eq. (2) in the limit where |x| ≫ |x ′ |. In doing so, we obtain a Schrődinger equation that captures the asymptotic behavior of the baryon wave function: The general solution is a linear combination of the Airy functions. Taking only the normalizable solution, where a is a constant. The asymptotic form is thus fixed entirely by the determination of the constant a. This can be done numerically provided that there is a region for which x is large enough so that the asymptotic form of the wave function of Eq. (14) is accurate while simultaneously being small enough that the variational wave function is accurate. If so, one can match to the two forms in the region of simultaneous validity of the two forms. It is easy to verify that such a region exists: a plot of the ratio of the two forms of the wave function must have a plateau with a value near unity over an extended region. Note that the size of such a region of stability depends on the accuracy of the variational wave function. It was largely for this reason, that we used as many terms as we did in the variational wave function. Such a stable region was found in the vicinity ofx = 3 and the asymptotic wave function takes the form: where the asymptotic form of the Airy function has been used.

III. COMPARISON OF SINGLE BARYON WITH MORE GENERAL APPROACHES
In this section, we compare our single baryon energy results with numerical results based on the Hamiltonian approach of Refs. [13] and [14]. These approaches are more general than ours as they are valid for arbitrary quark masses. Our approach complements these; while the regime is far more limited, the approach is more transparent and far less reliant on numerical work. Moreover, there is real value in having the ability to independently check the analysis of Ref. [13] in a regime where the results of infinite box size are clearly under control.
Ref. [13] analytically finds the most general Hamiltonian for large N c finite density QCD in the baryon sector for arbitrary quark masses in a finite box. The Hamiltonian for a box of size L is given by [13]: where the the trace is only over Dirac indices. The Hamiltonian is given in terms of density matrices ρ j (x, y) but can also be represented in terms of the conjugate space density matrix P Q (x, y) as The index Q represents integers and P Q (x, y) is defined as: U is a path-ordered (denoted P) spatial gauge field operator and |C[A] represents appropriately normalized coherent states Note that |C[A] has implicit dependence on gauge fields (denoted A) via g(x, y), the gauge-invariant generators of the fermionic coherent group [13] g(x, y) = ψ † (x)U (x, y)ψ(y) .
C(x, y) are the corresponding weights to g(x, y) and |0 is a color singlet reference state that satisfies the condition a n |0 =0, where a n is a member of a set of single particle annihilation operators chosen to ensure that the reference state has the baryon number of interest in the box.
To recap, the diagonal projection onto coherent states defined in Eq. (18) is for a particular gauge configuration; however, P Q reflects averaging over all possible gauge configurations. In this averaging, the only non-trivial structure that survives are the windings Q around the periodic box.
Note that ρ j (y, x) has the following properties: Here ρ j vac (x, x) is the density matrix associated with the baryon Dirac sea and gives a formally divergent contribution to the integral in the absence of regularization. B is the baryon number in the box, which can be written in terms of b, the baryon density as B ≡ bL.
The baryon Hamiltonian of Eq. (16) can only be solved numerically. To make the problem numerically tractable, we truncated the theory such that M in Eq. (16) is kept finite. This truncation introduces corrections to observables (including the baryon mass) that scale as 1 M . To obtain the observables in the full theory, the observables are calculated for different values of M and the result in the limit M → ∞ is obtained by extrapolation.
To formulate the problem numerically, the Hamiltonian is also discretized in space and the continuum results are obtained by extrapolation to the continuum limit. A latticized Hamiltonian using staggered fermions is discussed in detail in Ref. [13]. Note that the act of latticizing the theory introduces an ultraviolet cutoff thereby regularizing the theory and rendering the vacuum density matrix, ρ j vac (x, x), finite. Here, we only present the form of the density matrix ρ j (x, y) latticized using staggered fermions. These matrices satisfy the properties in Eq. (23) can be written in terms of a finite dimensional basis of orthogonal wave functions φ j n (x) for each j with the indices n, x = 1, 2...L s , where L s is the number of lattice sites. The discretized density matrix is with the following orthonormality constraints on the wave functions: The sum from n = 1 to Ls 2 in Eq. (24) is the lattice version of the Dirac trace of vacuum density matrix ρ p vac (x, x) in Eq. (23). The integral becomes regularized on the lattice and gives a finite contribution. The Hamiltonian of Eq. (16) was derived in the axial gauge representation of QCD [15]. It is well-known that in the Hamiltonian approach to either QCD or QED, Gauss' law appears as a constraint and not as an equation of motion. Ref. [15] explicitly constructed a QCD Hamiltonian, with the Gauss' law constraint built in a manner that only the physical states of the Hilbert space play a role.
Ref. [13], uses the constrained Hamiltonian of Ref. [15] and the coherent states approach of Ref. [16] to solve the Hamiltonian in the baryonic sector, which is the leading order physics in a 1 Nc expansion. This follows from the fact that the coherent states form an overcomplete set so that any physical state can be expressed as an integral over them and that the overlap between coherent states becomes sharply peaked as N c → ∞.
The Hamiltonian has a number of interesting features. First, for translationally invariant states, the classical Hamiltonian for the baryons exhibits volume independence and partial volume independence [13,19]. Second, color charges interact via a linear Coulomb potential of the form |x − y + QL|, where L is the size of the box and Q refers to the number of windings around the box. The sign of Q determines the direction of the windings. The presence of these windings is a result of the inability to gauge away all the spatial gauge degrees of freedom in a finite box. Specifically, there are N c − 1 spatial gluon windings that remain in the axial gauge after the imposition of the Gauss' law constraint.
Ref. [14], on the other hand, formulates the baryon Hamiltonian in infinite volume, where the windings represented by Q are absent; the form of the Hamiltonian is exactly the same as that of Eq. (16) except that Q = 0. However, in constructing the lattice regularized version of the resulting Hartree-Fock equations, Ref. [14] assumes that the correct approach is to use a box large enough such that finite volume effects on the baryon structure is small and with lattice spacing small enough that the baryon density can be easily resolved. For this purpose, Ref. [14] constructs a lattice potential, which has the form |x − y| in the continuum and infinite volume limits, while continuing to ignore the windings. Thus, the resulting latticized theory of Ref. [14] is equivalent to the M = 1 truncation of the Hamiltonian in Eq. (16) and, for small boxes, may have significant errors. One expects these errors to go to zero as the box size goes to infinity.

A. Comparison of Single Baryon Mass and Density Profile
We find that in the limit of large quark masses, there is excellent agreement between the variational calculation of section II and a numerical calculation based on the latticized baryon Hamiltonian of Ref. [13]. In Figs  It is useful to see the extent to which including the windings (i.e. M = 1) is crucial for getting interaction energies accurately for system with relatively large-sized boxes. From the plots in Figs. 3 and 4 it is clear that the baryons fit well within the boxes in the sense that the baryon density has dropped extraordinarily close to zero well before the edge of the box. Thus, one might be tempted to conclude that these boxes are big enough so that the effects of winding should be quite small. However, this is not the case. Comparing the interaction energy of the M = 1 baryon mass (which we label by E h ) with the interaction energy after M → ∞ extrapolation (labelled E H ) using the quantity EH −E h Thus, we find that the relative difference in interaction energy between the M = 1 result is approximately 7 per cent, which is quite significant given the degree to which the box fully contains the baryon.
The large difference between the M = 1 and M = ∞ results show the practical importance of the windings in realistic calculations even for very large boxes. Of course, the effect of the winding is expected to be large for small boxes and was introduced largely to explain the behavior of the system in small boxes including large N c volume independence [13]. To some extent, the large effect of the windings on the boxes computed here is rather surprising. Indeed, they are large enough that one might question whether there is something wrong with the formalism in [13]. However, the fact that the results extrapolated to M going to infinity agree so well with the results derived with the formalism of this paper, which are known to be valid for heavy quark masses, shows that this is not the case.

A. Introduction to the Multi-baryon Problem
Now that we have established that Witten's heavyquark formulation for single baryons is consistent with results based on the Hamiltonian approach of Ref. [13], we turn our attention to baryonic matter. It is believed that baryonic matter forms a crystal at large N c : the baryons are heavy and have a tendency lie in equilibrium at the potential wells created by the other baryons.
We consider infinite nuclear matter at low enough densities such that only interbaryon interactions between nearest neighbors are relevant.
Here, we assume that the many-baryon sector consists of distinct color-singlet clusters of N c quarks with identical wave functions, as in Witten's discussion of baryon scattering [17]. An analogous assumption was made in Ref. [18] for heavy quark, large N c QCD in 3+1 dimensions. This assumption is highly plausible but as noted in Ref. [18], it has not been rigorously established. A possible concern is that hidden color states could contribute at leading order. Thus, it is useful to test this assumption in a 1+1 dimensional setting. We can do that by compar-ing the results based on this assumption using our heavy quark formalism with numerical results using the methods of Ref. [13], which are known to be valid at large N c generally and do not involve any additional assumptions. Let us begin by deriving the interaction energy per baryon in a 1+1 dimensional crystal at low density using the heavy quark non-relativistic formalism. The energy per unit-cell of the low-density baryon crystal in terms of a reference baryon is: whereΨ i (x) is the wave function of the reference baryon, Ψ i±1 (x) are the wave functions of the nearest-neighbor baryons, the first term is the contribution to the intrabaryon kinetic energy, the second term is the contribution to the intrabaryon potential energy, and the last two terms are the contributions to the interbaryon energy. In this expression both λ and m q have been scaled out of the energy density using the following change of variables: By assumption, the density is sufficiently small so that the effect of interactions with next-to-nearest neighbors is negligible. The wave function for the baryon crystal is constrained by the Pauli exclusion principle, which precludes two fermions from simultaneously occupying the same position. This implies that the wave function of a reference baryon stateΨ(x) must be completely orthogonal to the the wave function of the neighboring baryon states Ψ(x +d) andΨ(x −d).
whereφ(x) is some normalized wave function. Here, we have assumed that densities are low enough that the reference baryon wave function is dominated by the single baryon wave functionψ(x) except for corrections centered aroundx = ±d. The choice of κ and α is determined by the orthogonality constraint of the neighboring baryons and by normalization of the baryon wave function. We find κ and α in terms of γ ≡ dxψ 0 (x)φ(x), which characterizes the overlap at the same spatial points of the neighboring baryon wave functions and the quantity A ′ ≡ dxψ 0 (x)φ(x +d), which characterizes the orthogonality of the neighboring wave functions: We have disregarded O(A 3 ) effects because A becomes exponentially small with decreasing density and in the low density regime, the higher order effects are parametrically small. Effectively, we are writing the interaction energy in the form A 2 f w d (where w is the characteristic width of a baryon and d is the separation of neighboring baryons) by developing an appropriate power expansion for f w d . To proceed, we first assume thatφ(x) =ψ 0 (x), which appears to be the most efficient way to orthogonalize neighboring crystal wave functions. We will show a posteriori that this assumption is correct in the sense that it will minimize the total energy up to corrections which are parametrically small at low density. The orthogonality constraint in doing so gives As can be seen from the explicit form of A(d), which we calculate in the following subsection, for parametrically low densities A(d) ∼ e −d 3 2 . Hence, it is plausible that in the low density regime, where only nearest-neighbor interactions are relevant, the choice ofφ(x) =ψ 0 (x) is indeed the correct one. In the section IV D, we introduce corrections such thatφ(x) = η(ψ 0 (x) + δ∆(x)), where∆ is a normalized wavefunction and argue that contributions to energy at O(δ) do not contribute at the order to which we work.

B. Calculation of overlap function A
Before proceeding with the calculation of energy, it is useful to calculate the overlap function A, which is defined in Eq. (29). Most of the contribution to A comes from the localized baryon density atx =d 2 . In the low density regime, the contributions to A, localized atx = 0 andx =d are exponentially suppressed. The contribution to the center comes from the tails of neighboring baryons, which have the form of Airy functions. Using the asymptotic form of Airy functions, we get where k ≈ 1.21 is the strength of the asymptotic wave function of Eq. (15). The function in the exponent has a sharp minima at x =d 2 and hence the integrand contributing to A has a maxima atx =d 2 . We use the standard steepest descent method to calculate the integral, which involves series expanding the exponent around the turning point, evaluating the denominator of A atx =d 2 and performing the integral. This gives with a leading order correction term that is proportional The correction is suppressed for parametrically large distances, which is the relevant regime for this paper.

C. Energy at low density
Here, we calculate the interaction energy per unit cell of a low density baryon crystal up to O(A 2 ) and leading order in O w d using units of Eq. (27). Note that w is the characteristic width of an isolated baryon and d is the interbaryon separation. We divide the energy into two parts: the first line of Eq. (26), which label the intrabaryon contributions, and the second line, which label the interbaryon contribution. The intrabaryon contributions contain a contribution from the kinetic energy and another from a potential energy. Note that the intrabaryon energy is not the same as in the free space baryon since the wave function has been modified due to the interactions with the other baryons.

Interbaryon Energy
The contribution to the interbaryon energy under the assumption that only nearest neighbor interactions are relevant is presented in Eq. (33).
The two terms in the equation correspond to the interbaryon energies in a unit cell associated with two nearest neighbors of a reference baryon stateΨ(x) in the crystal. Note that all terms in the interbaryon energy in Eq. (33) are of the form dx dȳ ρ(x)ρ(ȳ)|x −ȳ|, which is of the same structure as the Coloumb energy associated with a one-dimensional charge density ρ(x). Thus, for example, the first term in Eq. (33) has an effective charge densityΨ(x)Ψ(x +d). Moreover, note that in one spatial dimension, unlike in three spatial dimensions, the Coloumb energy of a charge distribution grows with the size of a system.
In Eq. (34), we present the effective charge density that determines the interbaryon energy in terms of the single baryon wavefunctionψ 0 (x), α and κ. Equation (29) completely determines the two parameters. Sincē φ =ψ 0 for a low-density baryon crystal, we have γ = 1 and A ′ = A. Therefore, α = A 2 and κ = 3 4 A 2 . To proceed, we note that the effective charge density relevant for the interbaryon contribution is localized. In  6) and (11). The same properties are true for the charge density centered atx =d. However, the effective charge density centered aroundx =d 2 behaves differently. For large interbaryon separationd, the width (in units of Eq. (27)) scales as (d − 2ǭ) 1 4 , which can be easily calculated from the standard steepest descent calculation in the integrand of Eq. (29). Also note that the total effective charge in the center peak is A for parametrically small baryon densities. Although the width of the central peak depends ond (i.e. the baryon density), it is much smaller than the interparticle spacing at low density. Thus in the low-density regime, these three regions of charge are well separated from each other with a spacing much larger then the width of each lump.
To calculate the leading contribution interbaryon energy per unit cell of a baryon crystal in the low-density regime, the effective charge densities can be approximated as delta functions with the appropriate baryon charges. The leading correction is the fact that the spatial extent of the charge density centered aroundx =d 2 grows withd. The size of this correction can be calculated using the approximation to the asymptotic form of the Airy function that was used to calculate the overlap integral A in Eq. (33). This integral can be calculated easily using the method of steepest descent. For parametrically low baryon densities, the total interbaryon energy becomes: where ∆Ē interbaryon Nc is the correction to the interbaryon energy associated with the fact that the center peak scales withd: The interbaryon energy above has the form -charge squared times a width -which is as expected. Also, note that the energy contributions due to the finite extent of the two regions of effective charge at the the ends is O(d 0 ) and hence is parametrically smaller at larged.

Intrabaryon Energy: Potential Contribution
In this subsection and the next, we calculate intrabaryon energy for the low-density baryon crystal at O w d . There are two contributions to this, the potential and the kinetic.
The potential contributions to the intrabaryon energy of the baryon crystal comes from Eq. (37). The effective charge density |Ψ(x)| 2 is presented in Eq. (38). Additionally, we present in Fig. 5 a schematic plot of |Ψ(x)| 2 − |ψ 0 (x)| 2 , the shift in the effective baryon charge density from the free case, |Ψ(x)| 2 , minus the isolated baryon density, |ψ 0 (x)| 2 , to emphasize the various charge enhancements. This quantity is even aroundx = 0 and hence we only plot the region 0 <x <d.
Furthermore, at O(A 2 ) the intrabaryon energy contributions for the baryon crystal come from a cross-term between the effective charge from the unperturbed baryon, ψ 0 (x), and the effective charge due to the perturbation induced by the near-neighbor baryons. To the extent that d is large, the original charge is well-localized at the origin, it creates an effective electric field which is constant for | x w | ≫ 1 (or a potential linear ind). With a linear field dependence, the potential energy is sensitive to the center of charge the effective charge due to the perturbation. Thus, we can approximate the charge densities in Eq. (38) as localized delta functions with the appropriate strengths, which are associated with the total baryon charges of each of the terms. Now for the case of the first term, which has a strength of 1 + 3 2 A 2 and is localized at the origin, the delta function approximation is problematic, since it is not in the region where the field from the unperturbed effective charge is linear. However, the potential energy from this term is of O(A 2w ), while the leading behavior scales as A 2d . Thus, this term can be neglected as being parametrically small. The contributions from the other two terms yield, As noted above, there are corrections to this of relative order w d from the perturbation in the effective charge due to κ. There are other corrections due to the fact that unperturbed effective charge has a spread and hence the effective field is not precisely constant. These corrections are exponentially small ind.

Intrabaryon Energy: Kinetic Contribution
The kinetic energy contribution to the intrabaryon energy at leading order in w d is presented in Eq. (40). Note that the first three terms are independent of the interbaryon spacingd, whereas the remaining terms are not; in the low-density regime, the contribution from these terms independent ofd are parametrically small. Furthermore, we ignore the O(1) term because it is the contribution to the single baryon kinetic energy and not the crystal interaction energy. We also ignore the last term because it is not nearest neighbor.
The remaining terms are trivially calculated using integration by parts, the asymptotic form of the Airy function in Eq. (14) and the steepest descent method to evaluate the integrals dx (∂xψ 0 (x))(∂xψ 0 (x∓d)), where the integrand is sharply peak aroundx = ±d 2 .
The kinetic energy contribution to the interaction energy then becomesĒ intrabaryon kinetic

The interaction energy
Finally, putting all the interaction energies together, the interaction energy per unit cell of a low-density baryon crystal in the heavy quark mass limit up to O(A 2 ) and to leading order ind is: Equation (42) contains as a factor (d − 2ǭ) 1/4 . We are interested in the limit of larged. So one could, in principle, expand (d − 2ǭ)  . This slower rate of divergence suggests that the expression will remain accurate for larger crystal densities.

D. Self consistency
Up until this point, we have assumed that the wave function for a baryon in unit cell was given by Eq. (28) withφ =ψ 0 . Such a choice seemed natural in that it orthogonalizes the wave functions associated with different unit cells in a very efficient manner. It turns out that this ansätz can be shown to be correct: the corrections to it yield changes in the interaction energy that are parametrically smaller than the result given in Eq. (42).
To show this, we introduce corrections toφ(x) of the formφ where δ is taken to be a small parameter,∆(x) is an arbitrary wave function and η is a normalization constant ensuring thatφ(x) is normalized. Without loss of generality, one can always take take∆(x) to be normalized and orthogonal to bothψ(x) andψ(x ±d). Any component in∆(x), which are not orthogonal toψ(x) and ψ(x ±d) simply serves to renormalize η and κ without altering the wave function.
Our ansätz thatφ =ψ 0 is self-consistently correct if it minimizes the energy up to corrections, which are parametrically small. This is true provided that for any choice of∆, since the leading term in the energy is parametrically O(A 2d 1 4 ). To proceed, one simply inserts the form in Eq. (43) into Eqs. (28) and (26) and expands to first order in δ. Our self consistency condition will be satisfied if the terms of O(δ) are parametrically of O(A 2d0 ). There are a total of 31 terms in the interaction energy at order δ -eight coming from the intrabaryon potential contribution, 21 from the interbaryon energy and two from the interbaryon kinetic contribution. Of these, four terms -two from the intrabaryon potential contribution and the two from the intrabaryon kinetic contribution -are each of O(A 1 ), the other 27 terms are all of O(A 2 ). It is straightforward to show that while each of the four terms of O(A 1 ) is nonzero, their sum is exactly zero.
Finally, there are O(A 2 ) terms from the interbaryon and intrabaryon energies. It can be shown on a termby-term basis that these terms all contribute at O(A 2d0 ) or less. In doing so, the orthogonality of∆(x)) with bothψ 0 (x) andψ 0 (x ±d) play a critical role since this implies that the integrated effective charge density due to cross-terms between∆(|barx) andψ 0 (x) yield zero. In demonstrating that each term contributes at O(A 2d0 ) or less, it is often necessary to separately consider the cases where∆(x) is localized on a scale ofw (that is,d 0 ) or in which it is spread over a parametrically larger fraction of the unit cell since the reason that the contribution to the energy is O(A 2d0 ) or less depends on the scale over which∆(x) is spread. However, while the reason that a given contribution is parametrically small may depend on the scale of∆(x), the fact that it is small does not.
The fact that all contributions at O(δ) either cancel or are parametrically of O(A 2d0 ) or less, shows that the our ansätz was self-cosistent to the order at which we are working. Therefore, the energy per unit cell up to parametrically small correction is whereĒ 0 Nc is the leading order contribution, which is stated in Eq. (42) andĒ 2 Nc is the energy at O(δ 2 ).

V. COMPARISON WITH DIRECT NUMERICAL SOLUTIONS
In this section, we compare the result of the lowdensity baryon crystal energy to the numerical results based on the Hamiltonian approach of Ref. [13], which as noted before is valid for large N c at any quark mass. Below, we plot the interaction energies for three different quark masses: mq √ λ = 10, 40, 100. We plot interaction energy (per color) versus the baryon density b = B L √ λ , where B is the baryon number and L is the physical size of the box. The points represent the results from numerical finite density calculations based on Ref. [13]; the lines represent the results based on our heavy-quark-lowdensity formalism.
The calculation was done for M = 10 and L s = 130 at a fixed volume for each of the quark masses. The mq √ λ = 10 calculation was done in a box of size L √ λ = 4 √ 2π and the mq √ λ = 40, 100 calculations were performed in a box of size L √ λ = 2 √ 2π. Obviously, there are numerical corrections to the interaction energy associated with discretization in space and the fact that M is finite. For the large quark masses involved the corrections associated with choosing a finite M = 10 and a non-infinite number of lattice sites L s = 130 are relatively small [13]. The discretization correction depends on the number of lattice points that define the wave function of a single baryon. Since, we probed different densities in a box of the same size (for each quark mass), the error associated with discretization is comparable for different baryon densities. We find that the agreement between the heavy-quark, low-density calculation of Eq. (42) and the results based on Ref. [13] is remarkabably good for small b. Indeed, it is surprising just how well it works. In deriving Eq. (42), rather strong assumptions were made -including the fact that contributions of O(A 2d0 ) were neglected as compared to terms of O(A 2d 1 4 ). One might expect these to become very important as b increases. Note from Fig. 7 that the approximation breaks down dramatically at certain fixed values, which depends on the quark mass. As noted above, this break down occurs atd = 2ǭ. It is remarkable, however, how accurate the large-mass-lowdensity expression remains even quite close to the break down value.

VI. DISCUSSION
In this paper, the nonzero baryon number sector of large N c QCD in 1 + 1 dimension in the limit of large quark masses was explored. The energy and density profile of individual baryons for different quark masses and the interaction energies of low density baryon crystals were computed. These observables were calculated directly using Witten's mean-field approximation, which simplifies greatly for infinite quark masses [17] since the result depends entirely on interactions involving single gluon exchanges. This approach, although valid in a limited regime of quark masses has the virtue of providing a simple intuitive picture for nonzero baryon density which complements the numerical approach Ref. [13]. The computed mean-field energy and baryon density profiles compare very well with ones calculated numerically based on the Hamiltonian of Ref. [13]. The agreement in the large quark mass regime is excellent, and demonstrate the significance of the effect of all the gluon windings even in cases where the baryon density is exceptionally small at the boxes' edge.
The analysis of low density baryonic matter developed in this paper relied on an ad hoc assumption that the hidden-color states -in which the wave function cannot be factorized into the product of color singlet baryons, play a subleading role at large N c . This assumption is highly plausible; however, unlike for the single baryon case, has not been explicitly demonstrated to be correct [18]. The fact that it agrees so well with the numerical solution, which was derived to be valid at large N c without making such an ad hoc assumption, is strong evidence that the assumption is in fact valid. Analogous reasoning was used for the case of large N c QCD with heavy quarks in 3+1 dimensions [18] in the low density regime; the success of the approach in 1+1 lends credence to the 3+1 calculation for which there is no numerical cross-check.
Finally, the success of our low density approximation with its ad hoc assumption in matching the numerical results based on Ref. [13], which was designed to work at large N c without such ad hoc assumptions, suggests a way forward in formally establishing the validity of the approximation used here and quite conceivably for the 3+1 dimensional case. Note that the approach of Ref. [13] was based on properties of a class of coherent states. It is highly plausible that a formal treatment based on the coherent states approach of Ref. [13] but taken for the heavy quark limit will yield the formalism developed here and thus replace our ad hoc assumption with something more solid on theoretical grounds. Moreover, it is plausible that some extension of such a method may work for the 3+1 dimensional case as well. While the general 3+1 dimensional case has a very complicated structure even at large N c due to the gluodynamics, in the heavy quark limit the gluons only contribute to produce static color Coulomb interactions and one might imagine an appro-priate coherent state formalism applying for the quark fields. These possibilities will be explored in future work.