Small Bits of Cold Dense Matter

The behavior of QCD at high baryon density and low temperature is crucial to understanding the properties of neutron stars and gravitational waves emitted during their mergers. In this paper we study small systems of baryons in periodic boundary conditions to probe the properties of QCD at high baryon density. By comparing calculations based on nucleon degrees of freedom to simple quark models we show that specific features of the nuclear spectrum, including shell structure and nucleon pairing, emerge if nucleons are the primary degrees of freedom. Very small systems should also be amenable to studies in lattice QCD, unlike larger systems where the fermion sign problem is much more severe. Through comparisons of lattice QCD and nuclear calculations it should be possible to gain, at at least a semi-quantitative level, more understanding of the cold dense equation of state as probed in neutron stars.


I. INTRODUCTION
Understanding quantitatively the properties of QCD at low temperature and high baryon density remains one of the most challenging problems in nuclear physics. It is increasingly important as it governs the behavior of neutron stars including their mass-radius relations (see Refs. [1][2][3][4] and references therein), cooling, and the emission of gravitational waves in neutron star mergers [5]. At low to moderate densities it is natural to model neutron star matter via a system of interacting nucleons, and much progress has been made along these lines [6][7][8][9] At very high densities the problem is also tractable: At low temperatures, the ground state will be a color-flavorlocked superfluid phase of quark matter with large superfluid pairing gaps [10,11]. At intermediate densities (several times nuclear saturation density) the situation is much less clear. One expects the dominant degrees of freedom to transition from nucleons to quarks and gluons, but the density at which this occurs remains very difficult to determine even qualitatively.
Studying the neutron star mass-radius relationship and cooling of neutron stars has been an important goal to constraining the equation of state through astrophysical observations [12]. The recent observation of two solar mass neutron stars [13], for example, has severely constrained possible equations of state at high density. Nevertheless these provide only some constraints on the equation of state, and less about the relevant degrees of freedom at high density.
In this paper we study the behavior of small numbers of baryons at high density (or equivalently small volumes) in periodic boundary conditions. It may soon be possible to study these systems in lattice QCD, as the sign problem for small baryon number is less severe. The sign problem grows exponentially with imaginary time and is proportional to the number of nucleons times the mass of the nucleon minus three halves the mass of the pion [14]. We evaluate these systems in both nucleonic models and quark models, and identify specific features that arise in the spectra as the degrees of freedom change from nucleons to quarks. These features rely on the relatively high momenta and short distances that arise in small periodic volumes. The most important of these features are shell closures for small numbers of nucleons and pairing in open-shell systems. The latter is important even for very small systems, in particular for the N = 4 systems we study.
The small lattice length L needed to simulate high densities with a small number of baryons may be subject to significant corrections from long-range physics. However, we can expect that some of the longest-distance effects due to pions should be the same in nuclear and QCD simulations as long as nucleon-pion interactions are consistently included in the nuclear Hamiltonians. In any case the smallest box considered in this study (corresponding to N = 4 at ρ = 0.48 fm −3 ) has a box size of L ≈ 2.03 fm.
It may also be possible to perform lattice QCD and nuclear simulations for unphysical heavy up-and downquarks, resulting in heavy pions. This would reduce the sign problem and the corrections due to finite volume at the cost of unphysical pion masses. Studying the transition even for high pion masses may be instructive as pion degrees of freedom may not play a very important role in the high-density phase transitions. See Ref. [15] for some studies of N N phase shifts at high pion mass.
Small volumes will typically result in large excitation energies, reflecting the wider spacing in the singleparticle spectra, as discussed below. This will be less true for comparison of different pairing symmetries that are degenerate in the free-particle limit. However, in general, both the quantum Monte Carlo (QMC) and lattice QCD simulations will converge more quickly with imaginary time for small systems.
Although studies of small systems are ill-suited to capture critical behavior and cannot precisely identify possible phase transitions, this work is well motivated because presently we lack even a qualitative understanding of how quark degrees of freedom might emerge at high density. The signatures we identify below suggest that these small systems could be quite valuable in this regime. Eventually one can add protons and/or hyperons to look at the density dependence of the symmetry energy and/or the presence of hyperons in neutron stars. In these initial studies, we concentrate on pure neutron matter.

A. Nucleonic Models:
To study nucleonic matter we consider nonrelativistic nucleons interacting via two-nucleon (N N ) interactions: where the N N interaction is taken as either the Argonne AV18 [16], the AV8 [17], or the local next-to-next-toleading order (N 2 LO) chiral interactions of Ref. [18]. At low densities such interactions should be able to faithfully reproduce the properties of QCD. Of course the three-(and eventually four-and many-) nucleon forces will be present, and eventually play a significant role. The simple spectral features we identify below will remain, though, even in a more sophisticated picture. These features depend primarily upon the single-particle states available to a nucleon in periodic boundary condition, indeed some are present even for noninteracting neutrons.
In this paper we consider only cubic simulation volumes with periodic boundary conditions. Other geometries, such as elongated volumes or different types of boundary conditions, might allow one to identify additional spectral features [19], but since present lattice calculations typically use cubic symmetry, we will adopt it in this study. Below we show results for different densities with various numbers of neutrons in cubes of length L on each side, and we show results as a function of the density ρ = N/L 3 . The nuclear interaction models described above are defined in the continuum. To maintain periodic boundary conditions, we add the contributions from periodic images for each pair: where r ij is the minimum separation in the periodic box and M = 1 to 2 images in each direction is sufficient to obtain periodic solutions since the N N interaction is at most of pion range. We note that the longest-range corrections from pions "wrapping around" the box are also present in lattice QCD simulations [20]. For heavier pion masses these long-range periodic images will play much less of a role.
For the larger number of neutrons (N > 4), we restrict ourselves to solutions fully symmetric under cubic rotations. For the smallest systems, we also investigate states that would correspond most closely to p-wave (N = 3, 4) or d-wave (N = 4) solutions in the continuum. The couplings to nonzero angular momenta prove quite interesting in comparing neutron potential models to quark models.
Calculations are performed using QMC methods: Either the Green's function Monte Carlo or Auxiliary Field Diffusion Monte Carlo methods. More details are described in [21]. The simulations are fairly simple as there are only a modest number of neutrons and the small volumes raise the energies of the excited states allowing for a more rapid convergence.

B. Quark Model:
For comparison, we also consider very simple quark models of high-density QCD in periodic boundary conditions. These models are not intended to be realistic or predictive of the behavior of QCD in this regime, but they should illustrate possible alternative behaviors when deconfined quarks are the dominant degrees of freedom.
We consider both free and interacting quarks in periodic boundary conditions. In general the models can be written as: In the free case we consider for the kinetic energies T i of the quarks relativistic and nonrelativistic dispersion relations For the interacting quark model, we shall assume that chiral symmetry is not fully restored when quark degrees of freedom first manifest in the spectrum and use a relatively large value of m = 300 MeV. Under these conditions we treat the quarks as nonrelativistic even for the small volumes that we consider.
The pair potential is chosen in order to reproduce the pairing pattern expected from a Nambu-Jona-Lasiniolike model [22] where interactions are antisymmetric with respect to color. Furthermore, for 3 colors and 2 flavors, we look for a color superconducting ground-state called the 2SC phase where two species are paired in a color and flavor anti-symmetric channel while the third does not interact with either and is effectively decoupled [23]. In the following we will consider the blue quarks to be the decoupled species.
For a given neutron number N the occupation numbers in flavor-color-spin space are chosen in order to satisfy the following constraints baryon number: When quark degrees of freedom are manifest, it is appropriate to neglect confinement at the short-distance scales of relevance to our small volumes. We assume that the average confining interaction per baryon is a function of the density only independent of baryon number at fixed density. We then compare the evolution of the energy with baryon number at fixed density to the nucleonic models.
The shell structure can be influenced by interactions, though, especially those that lead to pairing. To include this we shall consider a simple local potential of the form where i, j are multi-indices containing color, flavor and spin projection. In order to obtain the desired pairing structure expected in high density QCD, we choose a simple short-range interaction The matrix Λ ij is chosen to be anti-diagonal with entries equal to 1, µ is the reduced mass, and where a s and r e are the (S-wave) scattering-length and effective-range respectively. In our calculations we choose a s = 10 fm and r e = 0.1 fm in order to be close to the unitary limit. It is useful to compare relevant energy scales for 4 and 14 neutrons in periodic boundary conditions. In Table  I we compare the lowest finite energy (k = 2π/L) freeparticle modes in the box. Since the quarks are light, it takes substantial energy to raise them to higher momentum states, as indicated in the table. There are more degrees of freedom available, however, meaning that quarks will have a substantially lower Fermi energy for the case of weak interactions at high density.
In very simple models the strength of the quark interaction is adjusted to reproduce the N − ∆ mass splitting. The total splitting is 320 MeV, which can be compared to the single-particle energy splittings above. For example, for four neutrons, two nucleons have a momentum |k| = 1, while with quarks one can accommodate all quarks with |k| = 0 at the cost of twice the N − ∆ mass splitting. At low densities (large volumes) the fourneutron system would be preferred, but at high densities (small volumes) having all the quarks at |k| = 0 would be preferable. At present it is difficult to compute many-neutron states in lattice QCD because of the rapid growth in the number of correlators required and because the signal to noise ratio for small pion mass grows exponentially with baryon number. Very small systems of 4 baryons may be easiest to simulate in lattice QCD. For these systems we calculate states with different quantum numbers as an additional probe of hadronic versus quark degrees of freedom. We find distinctively different behavior even for 4 neutrons when comparing the nuclear and quark models for different quantum states.
Two neutrons in finite volume have been studied in lattice QCD (see for example Refs. [24][25][26]) and as two nucleons using QMC methods in Ref. [27]. These results are directly tied to the phase shifts of the neutron-neutron interaction via the Lüscher formula. Systems of three and four neutrons in external wells at low density have been studied in Ref. [28]; here we are interested in the behavior at high densities in periodic boundary conditions to mimic lattice simulations.
For three neutrons we study low-lying states with the quantum numbers of two neutrons with spin and total momentum zero, and with the extra neutron in a |k| = 1 state. The spin of this unpaired neutron can be oriented along or anti-aligned to the lattice equivalent of the angular momentum giving something similar to P 3/2 or P 1/2 states. As expected the former are slightly lower in energy due to the spin-orbit splitting in the neutronneutron interaction. The total (including center-of-mass kinetic energy) ground-state energies of three neutrons at different densities with the AV8 and AV18 interactions are compared with free neutrons with the same boundary conditions in Fig. 1.
For four neutrons we study states with two neutrons paired to total momentum and spin zero, and then either s-, p-, or d-wave pairing of the remaining two dominantly We find the initially surprising result (see Fig. 2) that neutron-neutron interactions favor pairing in the d-wave state for small box sizes. For free neutrons the three different (s-, p-and d-wave) states would be degenerate. The interaction between the two neutrons in |k| = 1 states coupled to zero total momentum dominates the spectrum. The relative momentum for this pair is quite large in these small volumes, in the region of the repulsive s-wave neutron-neutron interaction. The relevant phase shifts are shown in Fig. 3. At saturation density, two neutrons with momenta +1 and −1 are at a total energy of nearly 200 MeV, while at twice saturation den-sity the center-of-mass energy is 300 MeV. The strong s-wave repulsion disfavors the s-wave state, and the periodic boundary conditions favor the L = 2 state, as the d-wave pairing is symmetric across the periodic boundaries. That is, a pair orbiting with L = 2 feels an attractive interaction, while for L = 1, the periodic images interfere as the the relative coordinates r and L − r and consequently L · S are oriented in opposite directions. Results for the four-neutron calculations at different volumes are shown in Fig. 2 for the AV18 and N 2 LO N N interactions. The different states all have very similar energies at half saturation density, while the d-wave state is favored in all these models at ρ = 0.16 fm −3 and above. At twice saturation density the s-wave state is roughly 100 MeV higher than the d-wave state.
For free or paired quarks it is always advantageous to keep all the quarks dominantly in k = |0| states. The pairing energy that can be gained by promoting some quarks to higher momentum is small compared to the energy cost of promoting two or more quarks to |k| = 1 states. These conclusions are unaltered by the addition of a gluon-exchange spin-interaction (which historically was invoked to explain the N − ∆ mass difference) of the form [29] where S (2) ij is the tensor operator. The energy difference between the s-wave and d-wave states favors the s-wave as the ground state, since the s-wave is lower by ∆E per baryon ≈ 27 MeV and 45 MeV at 2 and 3 times nuclear saturation density, respectively. This ordering is opposite to that seen in the neutron calculations.

III. RESULTS FROM N = 4 TO 14 NEUTRONS
We also consider nuclear ground states from N = 4 to 14 with even numbers of neutrons paired to spin zero with full cubic symmetry. In the continuum, this would correspond to an s-wave superfluid, which is expected to be the ground state at low densities [30,31]. The special cases N = 2 and 14 correspond to filled singleparticle nuclear shells and hence are expected to have lower energy per particle than the remaining (open-shell) systems. Similar behavior has been observed in calculations of neutrons in external fields with either Woods-Saxon or harmonic wells [32,33]. This is confirmed by our numerical results as shown in Fig. 4. The rest mass of the nucleons are not included in this figure. Note that the differences between different particle numbers are quite significant at high density, of order 10 MeV per particle between adjacent N and of order 50 MeV per nucleon lower for N = 14 compared to N = 6. Such strong shell dependence, absent for free quarks, is also observed in the limit of strongly paired quark matter.
The nuclear model dependence is fairly small at modest densities, but increases substantially at the highest densities considered. This is illustrated in Fig. 5, where we plot the ratio of energies to Fermi gas energies for different particle numbers and densities for both the AV8 and AV18 interactions. The AV8 model only fits the lower partial waves, and is therefore less reliable at moderate to high densities. The difference between AV18 and AV8 is treated perturbatively in these calculations. The pattern of E/N versus N is the same for both interactions. In Fig. 6 we compare results of the nuclear model with free neutrons and also free and paired quarks at twice saturation density. The confinement energy in the quark model, which is assumed to be constant with density, is adjusted to match the nuclear result for 14 interacting neutrons. Note the dramatically different behavior versus baryon number in the nuclear and quark models, particularly for small N . The same behavior is observed also at three times saturation density in Fig. 7 For both quarks and neutrons we find that the addition of interactions does not change the qualitative behavior apart from the cases N = 16 and N = 6 at high density. The sudden increase in energy for N = 16 compared to N = 14, which is caused by the filling of the |k| = 2 momentum shell in the neutron case, has a similar nature for the quark model: At N = 16, it is preferable for the (down) blue quarks to fill the |k| = 2 shell instead of populating the (up) blue states as this will require breaking interactive pairs. The large energy at N = 6 for the quark model compared to N = 4 is due to the fact that for the latter all of the interacting quarks can reside in the |k| = 0 shell while at N = 6 two pairs of quarks have to fill the |k| = 1 state. This large shell effect of ≈ 310 MeV per baryon at 3 times saturation density will necessarily indicate the dominance of quarks degrees of freedom. The feature is less pronounced at twice saturation density where the, still large, gap is ≈ 120 MeV. This effect seems rather robust since even in the limit of noninteracting quarks the gap between N B = 6 and N B = 4 is still rather large, ≈ 170 MeV for ρ = 0.48 fm −3 and ≈ 120 MeV for ρ = 0.32 fm −3 .

IV. CONCLUSIONS
In this work we have studied the behavior of the energy per baryon in high density neutron matter for small volumes. These types of simulations can be valuable in understanding the relevant degrees of freedom at high baryon chemical potentials, in particular to begin to identify the regime where nucleon degrees of freedom dominate and where full QCD simulations of quarks and glu-ons are required.
While these simulation volumes are small, particularly for N = 4, and may suffer from important finite-size effects beyond the periodic pion potentials discussed here, further studies with larger quark masses and the N N potentials derived from them are warranted. Eventually lattice QCD studies with larger N at physical pion masses will be important.
Even for very small systems, we find a particular pairing pattern in neutron simulations. For four neutrons the d-wave state is favored at high density for the two neutrons in the |k| = 1 state, due to a combination of the sand d-wave phase shifts and the periodic boundary conditions. In contrast, for quark models the s-wave pairing is always favored to avoid filling the |k| = 1 states.
A clear signature of the dominance of strongly paired quarks as degrees of freedom can also be seen in the large energy gap between N B = 4 and N B = 6, in particular for paired quarks compared to nucleons at 2 (or 3) times saturation density. For the spatially symmetric states, we expect an increase of ≈ 120 MeV per baryon (≈ 310 MeV per baryon) to be compared with a much smaller gap of ≈ 4 MeV per baryon (≈ 15 MeV per baryon) for the interacting nuclear case.
Simulations of three-and four-nucleons (A = 3 and 4 nuclei) have already been completed with higher quark masses, so it is reasonable to expect that similar comparisons of neutron simulations to full lattice QCD could be completed in the near future. A direct comparison would require reasonable extractions of the N N scattering phase shifts for these quark masses. Such a direct comparison could have dramatic implications for the behavior of QCD at low temperatures and high baryon density and hence the mass-radius relation of neutron stars and the gravitational waves emitted during their mergers.