{'fraction_non_alphanumeric': 0.03403523207540113, 'fraction_numerical': 0.01107080076298762, 'mean_word_length': 4.412550607287449, 'pattern_counts': {'":': 0, '<': 0, '<?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': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Saturation Properties of Nuclear Matter with a Nonlocal Confining Solitons

We examine saturation properties of a quark-based picture of nuclear matter. Soliton matter consisting of nonlocal confining solitons is used to model nuclear matter. Each composite nucleon is described by a non-topological soliton as given by the Global Color Model. We apply techniques and concepts from the discription of crystal lattices. In particular, the Wigner-Seitz approximation is used to calculate the properties of the soliton lattice at the mean-field level. We focus on infinite nuclear matter at around standard nuclear matter density with the simplest one-parameter assumption for the gluon propagator. the saturation density and incompressibility are calculated as functions of the single parameter of the model


I. INTRODUCTION
The Global Color Model (GCM) 1 has been developed to model quantum chromodynamics (QCD) at energies low on the scale of high-energy particle physics [1]. The successes of the GCM in the hadronic sector [2] include the reproduction of chiral perturbation theory results [1,3,4], meson form factors [5] and spectra [6], and both, the soliton [7] and the Faddeev [8] picture of the nucleon. These successes motivate the application of this particular soliton model when a description of nuclear matter in terms of quarks and gluons is attempted.
Vladimir Naumovich Gribov efficiently used condensed-matter techniques in studies of quark confinement [9]. This is but one example of the recent convergence in language and methods between nuclear/particle physics and condensed-matter physics. In this paper we use the theory of crystal lattices in a description of nuclear matter in terms of quarks and gluons, with the GCM providing the physics of the primitive cell, which corresponds (for large values of the lattice constant) to an isolated composite nucleon.
To model nuclear matter in terms of a many-soliton problem where each individual soliton is accountable to the underlying quark and gluon degrees of freedom is an ambitious undertaking [10][11][12]. At low excitation the kinetic energy of the composite solitons can be neglected as a simplification. Due to the short-range repulsion between the solitons, a soliton crystal emerges as the lowest-energy solution at the mean-field level [13]. This crystal, however, is complicated relative to the structures familiar from condensed matter physics in that the lattice sites are themselves populated by extended objects (solitons), localized to some extent for large lattice constants. As the density increases, we expect a mean-field solution which no longer has localized centers of density.
An early success of condensed-matter physics is the Kronig-Penney model [14], which treats a crystal as a one-dimensional periodic set of square wells. This simple description allows analytic solutions to be found for the Bloch wave functions and for the energy bands usually associated with crystal lattices. In the case of nuclear matter, deeply bound light quarks call for a relativistic approach. Glendenning and Banerjee have solved the relativistic Kronig-Penney model and obtained analytic solutions for the eigen-value spectrum which retains the band structure [15].
In more complete geometries analytic solutions of the Dirac equation for the Bloch wave functions and eigen-energies are not available, so some kind of approximation is needed.
Here we apply the Wigner-Seitz approximation [16] to nuclear matter in the ground state. This entails the assumption of a spherically symmetric mean field, which takes into account the effect of the surrounding matter (lattice sites) on the primitive cell in a directionally averaged manner. This approximation may quite naturally take into account that nuclear matter is more like a fluid than a crystal [17].
As discussed above, the primitive cell will be given by the GCM, which admits soliton solutions with an intrinsically generated, extendedqq meson field. In addition, the individual GCM solitons are confining, as can be seen by the lack of poles in the quark propagator outside the region where the meson field is nonzero [7,18]. We have recently applied this kind of picture to draw attention to the intersection of the ground state band with the next unoccupied band at high densities [18,19]. The intersection represents a color insulator-toconductor transition in the model, and thus signals the onset of quark deconfinement. In this sense these ideas of condensed matter physics are applied in the spirit of Gribov [9], albeit to a model situation in place of full QCD.
Our goal in the present paper is to look at the saturation properties of this description more closely. The paper is organized as follows. In Section 2, we describe the model used for the primitive cell. The application of the GCM to nuclear matter is reviewed in Section 3. In Section 4 we present new results around saturation density, including nuclear-matter incompressibility. Section 5 contains the discussion and the summary.

II. GCM: THE PRIMITIVE CELL
As quarks and gluons are confined within hadrons, the description of low-energy nuclear matter has traditionally been accomplished in terms of effective hadronic degrees of freedom.
The effective degrees of freedom should be derived from the QCD Lagrangian. However, until an appropriate solution of QCD at low energies and large length-scales becomes available, it is necessary to use a model in place of full QCD. Clearly, we attempt to keep as much of the essential physics in the model as possible. An approach to systematize the needed approximations and to connect QCD to an effective hadronic field theory can be formulated in terms of functional integral methods. The strategy is to transform the integration variables from quark and gluon fields to hadron fields.
One particular implementation of the above ideas is in the framework of the GCM, which starts with a truncation of QCD [1], leading to the Euclidean action § Here j a ν (x) =q(x) λa 2 γ ν q(x) is a local quark current, with Euclidean Dirac matrices γ ν and Gell-Mann matrices λ a . The two-point gluon function, D µν , can be considered the phenomenological input point for the model. Using a Feynman-like gauge, D µν = δ µν D(x − y), the gluon propagator is particularly simple, and provides a single parameter function for the GCM. In (1), m is a current quark mass, which will be taken to be zero in the following. The GCM has the global symmetries of QCD, but lacks local gauge invariance. Our choice of the gluon propagator is dictated by the requirements of simplicity and quark confinement. We take a delta function in momentum space, The strength parameter α is the only input of the model. This form of the gluon propagator has been shown to produce quark confinement and reproduce nucleon properties [20].
To exhibit nonlocal quark-antiquark structures in the action, a Fierz reordering may be performed [21], which transforms the current-current term in (1) as Here, J θ (x, y) =q(x)Λ θ q(y) can be looked upon as a quark-antiquark bilocal current with quantum numbers specified by θ. The quantity Λ θ is a direct product of Dirac, flavor, and color matrices, and contains both, color-singlet and color-octet terms. We focus on the color-singlet sector in this work, ignoring correlations that correspond to diquark degrees of freedom.
To cast the partition function in terms of Bose fields, auxiliary nonlocal fields, B θ (x, y), are introduced, and the partition function is multiplied by where N is a normalization constant. After the transformation B θ (x, y) → B θ (x, y)+g 2 D(x− y)J θ (y, x), the action is bilinear in terms of the quark fields and the Grassman integration can be performed in the usual way. This yields the action in terms of bilocal Bose fields.
The replacement of the quark fields with Bose fields is, in principle, an exact functional change of variables. Observables calculated from the action (1) are not affected by the variable transformation, but are now expressed in terms of the Bose degrees of freedom, provided the entire sum over θ is kept. This is impossible in practice, and the truncation scheme used can be developed into a systematic method of approximation. To retain the chiral content of the QCD action, at least two Bose fields need to be kept.
The classical vacuum configuration B θ 0 is identified by δS/δB θ = 0. This produces a quark self-energy, Σ(x − y) = Λ θ B θ 0 (x − y) satisfying a Schwinger-Dyson equation. In momentum space This can be considered an integral equation for the self-energy amplitudes A(p 2 ) and B(p 2 ). Numerical solutions for these amplitudes are now available at varying degrees of sophistication [22]. Our choice for this study is governed by simplicity within the context of the requirement of confinement. The delta-function gluon propagator (2) simplifies the integral equation (5) resulting in an algebraic equation. Analytical solutions for A and B can be found as [23] A(p 2 ) = 2 where the input parameter α controls the coordinate space width of A and B. The lack of solutions to the equation p 2 + M 2 (p 2 ) = 0, (where M = B/A is the dynamic quark mass) indicates that (6) produces quark confinement as there is no on-mass-shell point, thus the propagation of a quark in the normal vacuum is prohibited [18]. It is important to note that the amplitude B(p 2 ) plays a dual role in the model: it also acts as the distributed vertex for coupling the quarks to theqq Goldstone modes [7,24]. The fluctuations B θ (x, y) = B θ (x − y) − B θ 0 are identified as the propagating Bose fields. If the color-singlet scalar-isoscalar and pseudoscalar-isovector fluctuations are retained, the formalism can be adopted to the requirement of chiral invariance by the variable transformation where r = x − y and R = (x + y)/2, are relative and cm-like coordinates, respectively, f π is the pion decay constant, and it has been assumed that the on-shell form factor B can also be used off-shell. As a further simplification, the φ = 0 point on the chiral circle can be fixed. In this case the radial fluctuations away from the chiral circle coincide with the scalar-isoscalar field variable prior to the transformation. In the numerical work that follows the single fluctuation field χ corresponding to this situation will be used. Letting m → 0, the action (up to a constant) can be written as a sum of fermionic and bosonic terms: The chemical potential (µ) dependence of the fermion term in equation (8) ensures that a meson source from the valence quarks will be generated [7]. The U( χ 2 ) term is the effective meson self-interaction [1]. For µ = 0 the inverse quark propagator takes the form and the saddle-point configuration is at χ = f π [20].
Since G −1 (x, y) is time independent, stationary eigenstates of the form u j (x) can be obtained from a self-consistent Dirac equation, which in coordinate space takes the form The new χ = χ − f π , and as p 4 = q 4 = iǫ j , where ǫ j is the energy eigenvalue, the meson vertex B has an energy dependence. It can also be seen that a wave-function renormalization appears with the renormalization constant Z j given by [7] The meson field equation δE δχ = 0 may be summarized as with the meson source provided by the valence quarks according to Equations (10) and (12) form a system of coupled differential equations for the quark wave functions and the meson field, which need to be solved selfconsistently, with the appropriate boundary conditions.

III. INFINITE NUCLEAR MATTER
The use of soliton matter to model high-density nuclear matter originated in the mideighties [10][11][12]. The new feature of the present work is that the GCM soliton represents the primitive cell. Some of the above authors also apply the Wigner-Seitz approximation [16], to which we turn next.

A. Wigner-Seitz approximation
As a means of describing nuclear matter, we consider an infinite collection of GCM solitons. At the lowest energies the solitons are expected to arrange themselves in a crystal lattice [15,13]. Accordingly, the single-quark eigen-energies will develop into energy bands. For simplicity, we assume a simple cubic crystal (sc). For a periodic lattice, the Dirac wave function must be invariant to a lattice translation, so the solutions must have the Bloch form [25] u lat m (r) = u m (r) e im · r , where m = (m x , m y , m z ) is the lattice momentum and u m (r) is a Dirac spinor which has the periodicity of the lattice. To solve for the Bloch functions we employ the Wigner-Seitz approximation [16]. This amounts to taking the primitive cell as given by the geometry of the lattice (a cell bounded by perpendicular bisectors of the lattice vectors as shown in Figure 1) and replacing it by a spherical cell of radius R and solving for m = 0 in (14). The full solution u lat m (r) is then approximated by u 0 (r) e im · r . Changing the density will be implemented in our calculation by varying the cell radius R.
The Wigner-Seitz approximation places specific boundary conditions on the Dirac spinors. These conditions express the requirement that the upper component of the wave function must be periodic and anti-periodic for the bottom and the top of the band, respectively. We focus attention on the lowest-energy state of the band, for which the above, together with the r = 0 boundary conditions, implies where g(r) and f (r) represent the radial parts of the upper and lower components of the Dirac wave function, respectively. In addition, the meson field of equation (12), which also appears in the source term of the Dirac equation (10) must now be periodic in r, so that where χ(r) is the radial part of the meson field.

B. GCM on the Lattice
It is convenient to solve the Dirac equation in momentum space, while the nonlinear Klein-Gordon equation is easier to handle in coordinate space. We seek solutions of the Dirac equation with the boundary conditions (15). These can be incorporated using a three dimensional Fourier expansion. Expanding A, B, u, and the meson-field source, we integrate over y and use orthonormality to get an equation for the Fourier components of the quark wave function where g n = g(k n ), f n = f (k n ), and To get this final form we have written the Dirac wave functions as where k n = nπ R and we used the spherical symmetry of the Wigner-Seitz cell. If m = 1, 2, . . . , M, then equation (17) is a 2M-by-2M eigenvalue problem for the energy eigenvalue ǫ. The quantity B/A plays the role of a dynamic mass and the scalar part of the self-energy B also acts to couple the quarks to the meson field via (18). The self-energy terms have an ǫ dependence which makes this a highly non-linear problem. The V l term (l = 1, 2) represents the Legendre coefficient for the meson field in the presence of the distributed coupling B.
One needs to solve the Dirac equation (17) and the Klein-Gordon equation for the meson field (12) selfconsistently.
To solve for a soliton lattice, we first pick a starting meson field and search for the lowest energy eigenvalue of Eq. (17). This means finding the energy ǫ which makes the determinant of equation (17) vanish. We start at ǫ = 0 and work upwards in energy until the determinant changes sign. We then use the bisection method to find the root. Care must be taken so that the initial steps are sufficiently fine in ǫ not to miss the lowest root. With the root in hand, we can solve for the Fourier components of the Dirac wave functions. For this we first perform a lower-upper triangular decomposition and use inverse iteration [26]. The momentum-space meson-field source term is constructed from the Dirac wave functions and we transform the source to coordinate space for use in the nonlinear Klein-Gordon equation for the meson field (12). To solve this nonlinear equation, we treat equation (12) as a functional of χ and use Newton's method. Once the Klein-Gordon equation is solved for the new meson field, we start over with the Dirac equation in this modified meson field. We iterate until convergence of the quark wave functions is achieved, which takes between three to six iterations to reasonable accuracy.

C. Energy Bands
As known from condensed-matter physics, each soliton of the lattice contributes one level to each energy band [25]. In the Wigner-Seitz approximation we need to calculate only the energy for the case when the crystal momentum m in equation (14) is zero. The Dirac-Bloch wave function, in the Wigner-Seitz approximation, for an arbitrary crystal momentum is and we use (20) to calculate the expectation value for the square of the Dirac Hamiltonian to estimate the lattice momentum dependence of the energy levels in the band as where ǫ bot is the energy of the bottom band.
To obtain the possible values of m we refer back to the underlying lattice structure. For a simple cubic crystal of N solitons and sides of length L = 2RN, the allowed values of the component of the lattice momentum in the direction of any of the three axes are with the top of the energy band corresponding to m = Nπ L = π 2R . Thus for the top energy band we obtain We have performed the same estimate assuming body-centered and face-centered cubic lattices. This variation in assumed crystal structure introduces a less than 10 % uncertainty in our results.
The lowest band (l = 0, l ′ = 1 and j = 1 2 ) is labeled 1s 1/2 . The next lowest (at around standard nuclear density) band has nonzero orbital angular momentum in the large component of the Dirac wave function (l = 1, l ′ = 2, and j = 3 2 ) and is labeled 1p 3/2 . The next band is again an s-state, corresponding to a radial excitation. For very low density (R −→ ∞) the energy bands shrink to single levels and in the limit reproduce the energies of a single soliton. As the density increases, the bands spread out and approach each other. The intersection of the energy bands as the density is increased beyond twice standard nuclear density is discussed elsewhere [19]. However, properties of the model at around saturation density have not been examined in detail in Ref. [19]. This is the goal of the present study.

IV. RESULTS
Here we focus on the properties of the lowest energy band (ground state band) at around standard nuclear density, ρ ≈ ρ 0 = 0.17 fm −3 . With one soliton in each Wigner-Seitz cell, the energy of each primitive cell is given by the soliton energy. As mentioned above, the kinetic energy of the center-of-mass of the soliton is neglected in the present approximation. The total soliton energy has a quark contribution which contains the quark kinetic and potential energies. Assuming three-fold degeneracy of the energy levels (three colors), the quark energy thus equals three times the lowest quark eigen-energy (ǫ) in the self-consistent meson field. In addition, the kinetic and potential energy associated with the meson field contribute to the soliton energy. The meson field has a kinetic energy given by The potential energy of the meson field is given by where U(χ) is the meson self-interaction from eq. (8), expressed as a function of the fluctuation field χ = χ − f π . U(χ) has the usual "Mexican hat" shape [20]. The total energy takes the form This quantity is plotted against the density in terms of standard nuclear density in Figs. 2-4 for different values of the parameter α. We selected a range for the values of α earlier [27]. It was argued that while reproducing the pion decay constant f π closely is particularly important for studies at the hadron level, for the present work, which is concerned with the density of nuclear matter, the root mean square (rms) radius of the proton constitutes the most sensible benchmark. It is not realistic to expect that this simplified one-parameter model, which does not explicitly take the pion field into account, should simultaneously fit both of these quantities. Since the explicit pion cloud is neglected in the model and it is normally assumed that the pions will increase the value of the rms proton radius by about 20-30 %, we need to choose α to give an rms proton radius ≈ .6 or .7 fm. We are able to achieve this for α ≈ 1.35 ± .1 GeV. Thus, we display calculated results for the values α = 1.25, 1.35 and 1.45 GeV in the following.
In Fig. 2 we show the quark contribution and the total soliton energy for α = 1.25 GeV. On one of the calculated points we indicate a typical uncertainty we associate with our computation. The main source of this uncertainty, which is not more than a couple percent, is the freedom in prescribed tolerances at different stages of the calculation. We see that as the density increases from zero, the ground-state energy develops a minimum. The low-density attraction between the solitons is a consequence of the boundary conditions on the quark wave functions (15). In particular, the upper component of the quark wave function is forced to have less curvature than in the case of a single soliton, leading to a lower value of the quark kinetic energy. At higher densities, where the solitons and the quark wave functions begin to overlap, the resulting repulsion overcomes the attraction and the ground-state quark energy starts to increase. For this value of α there is a minimum in the total energy at ρ min /ρ 0 = 0.66. To calculate the minimum we performed a second-order fit to the calculated points around the minimum, up to ρ/ρ 0 = 1.5. The value of the total energy at the minimum is 1844 MeV. However, our scheme does not take into account the center-of-mass kinetic energy associated with the quarks. An estimate of the soliton mass (on the lattice) can be obtained by subtracting this motion; M s = [E 2 − P 2 ] 1/2 , where P 2 = 3 p 2 , with p 2 being the expectation value of the square of the quark momentum. A numerical estimate along these lines leads to a soliton mass ≈ 15 % below the soliton energy, yielding ≈ 1560 MeV for the soliton mass for α = 1.25 GeV. Scaling back to α = 1.05 GeV (used in the earlier singlesoliton work) would decrease this by another 20 % to ≈ 1, 200 GeV, not too far from the average mass of the nucleon-delta system. Explicitly including the pion will lower this value even further [20].
Figs. 3 and 4 contain the same information for α = 1.35 and 1.45 GeV, respectively. In Fig. 3 we observe a minimum of the total energy at ρ min /ρ 0 = 0.73. Fig. 4 yields a minimum at ρ min /ρ 0 = 0.98. We have included all calculated points below ρ/ρ 0 = 1.5 in all these fits. (We carried out calculations at higher densities for each value of α in connection with our studies of the QGP transition [19]. However, in the present work we focus on the secondorder behavior of the energy versus density curve near its equilibrium value (minimum). We estimate an uncertainty in the position of the minima on the order of 15%, as a consequence of combining uncertainties in crystal structure, tolerances and fits. Finally we calculate the incompressibility defined as around the observed minima. The resulting values are summarized in Table 1, together with the pion decay constant, the root mean square proton (single soliton) radius and the value of the density where the minimum occurs for each value of α. As α increases, the root mean square radius of a single soliton is moving away from the experimental value of the rms radius of the proton, r 2 1 2 ≈ 0.83 fm. However, as discussed above, the effect of the pion cloud is expected to increase the root mean square radius by about a third. The incompressibility appears to be in the generally accepted range for the larger values of α. The ≈ 20 % uncertainty in this quantity is again due to a combination of crystal-structure, tolerance and fit-induced uncertainties. As mentioned above, the experimental value of the pion decay constant, f π = 93 MeV is best reproduced in our simple single-parameter description for α = 1.05 GeV. The value of f π is important for calculations in the hadronic sector. It is anticipated that with a more complicated gluon propagator with several parameters, the quantities f π , ρ min and K can all have reasonable values simultaneously.

V. DISCUSSION
We have calculated saturation properties of nuclear matter with a nonlocal confining soliton model, the Global Color Model (GCM). The encouraging successes of the GCM on the hadronic level triggered the use of the model for the description of nuclear matter in terms of the underlying quark and gluon degrees of freedom. Since there is a large number of degrees of freedom to deal with in each individual soliton, the soliton matter picture of nuclear matter is necessarily restricted to the mean field level and the Wigner-Seitz approximation is utilized. This model provided robust conclusions on the QGP transition at high density [18,19]. We therefore felt it necessary to investigate the properties of the model at around standard nuclear matter density.
A reasonable reproduction of nuclear matter saturation properties has been found. Some caution is necessary though when interpreting the physical significance of this agreement, since the kinetic energy of the soliton as a whole has been neglected in our calculation, while the Fermi kinetic energy of the nucleons plays an important role in nuclear saturation. Another reason while the agreement does not signal a deep physical understanding of nuclear matter at the quark level is the fact that we used an extremely simplified, "bare-bones" oneparameter gluon propagator. Furthermore, explicit pions were not included in the model. We find it remarkable that such a simple effective description nevertheless gives values of the equilibrium density and of the incompressibility in the right ball park. We wish to inform the community about the existence of such a powerful effective description. It is of interest to calculate in-medium properties as a function of the density with the model. We made the first steps in this direction elsewhere [27] by calculating the axial-vector coupling constant and correlation function in medium. It should be kept in mind that more realistic gluon propagators are now available in the literature [2], and should improve the model and accommodate a simultaneous fit of the pion decay constant and the root mean square radius of the proton. A more complete model requires the treatment of explicit pion degrees of freedom for inclusion into the calculation of the root mean square proton radius. Developments along these lines will also facilitate the exposure of the chiral content of the model and direct comparisons to QCD-based calculations. With such improvements the model can become a practical tool in the description of high-density strongly interacting matter.