Phase separation in lowdensity neutron matter
Abstract
Lowdensity neutron matter has been studied extensively for many decades, with a view to better understanding the properties of neutronstar crusts and neutronrich nuclei. Neutron matter is beyond experimental control, but in the past decade it has become possible to create systems of fermionic ultracold atomic gases in a regime close to lowdensity neutron matter. In both these contexts pairing is significant, making simple perturbative approaches impossible to apply and necessitating ab initio microscopic simulations. Atomic experiments have also probed polarized matter. In this work, we study populationimbalanced neutron matter (possibly relevant to magnetars and to density functionals of nuclei) arriving at the lowestenergy configuration to date. For small to intermediate relative fractions the system turns out to be fully normal, while beyond a critical polarization we find phase coexistence between a partially polarized normal neutron gas and a balanced superfluid gas. As in cold atoms, a homogeneous polarized superfluid is close to stability but not stable with respect to phase separation. We also study the dependence of the critical polarization on the density.
pacs:
21.65.f, 03.75.Ss, 05.30.Fk, 26.60.cI Introduction
The phase structure of polarized (or populationimbalanced) Fermi gases is very rich. The energetically favored phase depends on the masses and relative populations of the species, the interaction strength, and whether the interactions are long ranged. Polarized gases can exist in a variety of physical systems, including ultracold trapped atomic Fermi gases (1), dense hadronic phases near the cores of neutron stars (2); (3); (4), and possibly quark matter (5).
In this paper we focus on the possibility that neutrons in the inner crust of neutron stars are spinpolarized (6). Theoretical calculations suggest that this region, which is about a kilometer thick, features unbound neutrons while the protons cluster in neutron rich nuclei. While interactions with nuclei affect the properties of the unbound neutrons (7); (8), in this work we consider an isolated neutron fluid and compute the conditions for the coexistence of an unpolarized phase and a polarized phase at densities roughly times smaller than the nuclear saturation density. A polarized phase in the neutron star crust will have a larger specific heat and therefore its presence may affect the time scales associated with observed crustal transient phenomena (9). Neutron superfluidity also affects thermal transport properties (10); (11). If the polarization is large enough all effects of superfluidity are expected to disappear.
At the low densities we consider, the dominant term in the interaction between spinup and spindown neutrons is an wave potential with a scattering length fm and range fm. In a balanced system this attraction drives the formation of Cooper pairs of spinup and spindown neutrons of opposite momenta resulting in a neutron superfluid whose properties can be understood qualitatively using BCS theory. All the neutrons are paired and the number of up and down species is equal. Even at densities much smaller than the nuclear saturation density, is much larger than the interparticle spacing and obtaining quantitatively reliable results requires nonperturbative techniques. Detailed Quantum Monte Carlo (QMC) simulations have calculated the energy density up to densities equal to or larger than the nuclear saturation density (12); (13). At low densities, the wave pairing gap obtained peaks at nearly MeV.
Magnetic fields tend to split the up and down Fermi surfaces thus disrupting the pairing and creating a polarization. For magnetic fields () smaller than Gauss, pairing wins (assuming no effect coming from the lattice of nuclei), while larger fields can create a polarized phase. In BCS theory, the critical field, or the resulting chemical potential split ( where is the neutron magnetic moment), is called the ClogstonChandrasekhar point. At this point, the BCS superfluid can coexist in equilibrium with a partially polarized normal phase which is essentially a gas of free neutrons. Numerically where is the gap in the BCS phase. A weakly polarized gas of neutrons confined in a volume forms a mixed phase (14) consisting of the BCS and normal phases with volume fractions depending on the polarization. Hartree corrections substantially modify the ClogstonChandrasekhar point even for modest couplings , where is the Fermi wave vector (15).
At strong coupling (mathematically given by ) the energies of both the superfluid and the normal phase need to be computed nonperturbatively. In this paper, we perform these calculations for the normal phase and find the critical polarization at which it can coexist with the superfluid phase. We also compare the energy of a phase separated state with the previously evaluated (6) homogeneous polarized superfluid phase. In this paper we only consider competition between homogeneous and isotropic phases and ignore the possibility of promising but still unconfirmed scenarios (like wave (16) and inhomogeneous phases (17)).
As mentioned above, our microscopic ab initio calculations may impact the observed behavior of neutron stars with extremely large internal magnetic fields. It is interesting to concurrently examine possible connections of such simulations with terrestrial experiments. Our results may be useful for constraining Density Functional Theories (DFTs) that are employed to study properties of heavy nuclei. Equation of state results at densities close to the nuclear saturation density have been used for some time to constrain Skyrme and other other density functional approaches to heavy nuclei (18). The potential significance of such calculations has led to a series of publications on the equation of state of lowdensity neutron matter over the last few decades (19); (20); (12); (21); (22); (23); (24); (25); (26); (27); (28); (29). The densitydependence of the gap in lowdensity neutron matter has also been used to constrain SkyrmeHartreeFockBogoliubov treatments in their description of neutronrich nuclei (30). Similarly, an ongoing project is studying the limit of extreme polarization in neutron matter (31). This limit is known in condensedmatter physics as the “polaron” problem: one impurity embedded in a sea of fermions. Quantities of interest here include the polaron binding energy and its effective mass, which can be used to constrain Skyrme and other functionals. In this line of study, DFTs can be constrained by the equation of state of polarized unpaired neutron matter that we report here.
Another connection of our calculations with experiment may be possible through the field of atomic ultracold Fermi gases. There, the scattering length of the interaction between two hyperfine states near a Feshbach resonance can be tuned by changing the magnetic field. Most experiments have been performed in the region where the scattering length is infinite (unitarity limit) and the effective range is negligible (32); (33); (34); (35); (36). QMC calculations performed with the same techniques that we use here in the same limit give results that agree well with the experiments for both unpolarized (37); (38); (13) and polarized (39); (40); (41); (42) Fermi gases and this gives us confidence in applying the same techniques to the problem with finite and . More directly, by tuning the scattering length and the effective range (43) it may be possible to obtain a system with similar values of and as the inner crust and test our results.
Ii Construction of the mixed phase
The conditions for a polarized phase () to be in phase equilibrium with an unpolarized superfluid () phase are as follows:

The pressures of the two phases should be equal.

The average chemical potential of the two species in the phases, , should be equal.

The chemical potential splitting, , should be less than the gap in .
We construct the mixed phase between the superfluid and the partially polarized normal phase (). The pressure in a polarized phase depends on the net density of fermions () and their relative fraction where (or ) is the majority species.
Motivated by the unitary Fermi gas (44); (45), we write down the energy density as:
(1) 
For convenience we denote the energy per particle by and the corresponding energy for a free gas by . In terms of , .
For the unitary gas, the discussion of mixed phases is simplified greatly since the function is only a function of . In neutron matter, the interaction potential has inherent scales (these can be thought of simplistically as and ). To satisfy the conditions of phase equilibrium, we need to calculate the chemical potentials and the pressure, which include the functional dependence of on both and . The chemical potential is:
(2) 
The chemical potential splitting is:
(3) 
The pressure is simply , where .
For the unpolarized superfluid (), the function depends only on the density () and we have:
(4) 
To organize the discussion it is convenient to start from the completely polarized () gas. Here the interactions between the only species present are weak and the ground state is well described by a normal phase. (wave pairing and interactions will modify the energy of the fully polarized neutron gas from a free Fermi gas. As we discuss in Sec. III, we include the effect of wave interactions exactly, but ignore wave pairing.)
Now consider the following thought experiment. Keeping the density constant, change a few fermions to fermions so that changes. For small we expect that the phase obtained will be a partially polarized normal phase. As we increase the fraction of fermions, it is possible that some polarized superfluid phase () becomes the favored state of matter. Another possibility is that at some fractional density , it becomes favorable to form a mixed phase between the partially polarized normal phase () and the unpolarized superfluid phase (). In the present section, we will focus on the mixed phase between and . In Sec. III we will calculate the energy for the polarized normal gas () phase and in Sec. IV we will look at the competition between the homogeneous polarized superfluid () of Ref. (6) and the mixed phase.
Equating the chemical potentials and the pressures gives two equations in two variables, and . is the total superfluid density and is the total density of the normal phase. For one can find the coexistence curve in space. Suppose that a volume fraction is occupied by the phase and the rest () by the superfluid phase. Then:
(5) 
(6) 
where .
The parametric dependence on can be eliminated to find along the coexistence curve. The coexistence curve (for example, see dashed curve in Fig. 2) should not be seen as the continuation of the curve at constant density, since the density changes along it. It is simply a projection of the coexistence curve in space, on the plane.
To proceed with the calculation of the coexistence curve, we need to calculate for the phase and for the phase. We use Quantum Monte Carlo techniques to calculate the energies at a few values of and interpolate to determine the functional dependence. This is discussed next.
Iii Green’s Function Monte Carlo simulations
The Hamiltonian for lowdensity neutron matter is:
(7) 
where is the total number of particles. The neutronneutron interaction is not purely wave but still somewhat simple if one considers the AV4’ formulation: (46)
(8) 
In the case of =0 (singlet) pairs this gives:
(9) 
However, it also implies an interaction for =1 (triplet) pairs:
(10) 
Ref. (47) explicitly included such wave interactions in the samespin pairs (the contribution of which was small even at the highest density considered), and perturbatively corrected the pairs to the correct wave interaction.
In these calculations it is customary to first employ a standard Variational Monte Carlo simulation, which minimizes the expectation value of the Hamiltonian given a variational wave function . At a second stage, the output of the Variational Monte Carlo calculation is used as input in a fixednode Green’s Function Monte Carlo simulation, which projects out the lowestenergy eigenstate from the trial (variational) wave function . This is accomplished by treating the Schrödinger equation as a diffusion equation in imaginary time and evolving the variational wave function up to large . The ground state is evaluated from:
(11)  
evaluated as a branching random walk. The fixednode approximation gives a wave function that is the lowestenergy state with the sames nodes (surface where = 0) as the trial state . The resulting energy is an upper bound to the true groundstate energy.
The groundstate energy can be obtained from:
(12) 
The variational wave function is taken to be of the following form:
(13) 
In our earlier works we studied superfluid neutron matter and therefore used a trial wave function of the JastrowBCS form with fixed particle number. (13); (47); (6) In the present work we are attempting to determine the stability of different phases. In order to do this, we have completely mapped out the energy of a normal (i.e. nonsuperfluid) neutron gas as a function of the polarization at different densities. Thus, the we used in Eq. (13) describes the particles as being in a free Fermi gas (i.e., lets all the correlations lie within the Jastrow functions). In other words, the wave function is composed of two Slater determinants (one for spinup particles and one for spindown ones):
(14) 
where
(15) 
and
(16) 
The primed (unprimed) indices correspond to spinup (spindown) neutrons and . is the antisymmetrizer and .
The Jastrow part is usually taken from a lowestorderconstrainedvariational method calculation described by a Schrödingerlike equation:
(17) 
for the oppositespin and by a corresponding equation for the samespin . Since the and are nodeless, they do not affect the final result apart from reducing the statistical error.
We calculate groundstate energies at different total number densities [], more specifically at , , and fm. To put these densities into perspective we can compare them to nuclear matter saturation density: they are , and percent, respectively, of fm. These simulations were performed at values of the relative fractions chosen specifically in order to ensure a full coverage of the axis. More specifically, simulations were carried out at relative fractions of , corresponding to particle numbers of , respectively. These systems are quite large (and therefore computationally demanding) so as to ensure a minimization of finitesize effects (see also Refs. (13); (47); (48)).
The results are shown in Fig. 1 as points (circles, squares, and diamonds, respectively), along with cubic fits to the microscopic data. The latter will be useful to us in the following section, when we try to check the relative stability of different phases. To facilitate the use of these results in connection with Eq. (1), we have divided the groundstate energies with the energies of corresponding free spinup Fermi gases, and raised the result to the 3/5 power. When we see that the values are very close to 1, though still above it: this is due to the samespin wave interactions. Similarly, as we increase the density the energy is decreased, a fact that, as we see in Fig. 1, holds at any relative fraction . The statistical errors of the Quantum Monte Carlo results are smaller than the symbols shown in the figure. Also shown to the right of the figure are the QMC results for an upolarized superfluid at the three densities of interest (hollow circles, squares, and diamonds, respectively). The results for and were taken from Ref. (47); the value at was reoptimized for the purposes of the present work.
Iv Phase competition
Having arrived at the functional dependence of the energy of a normal neutron gas on the relative fraction (thereby getting information on the function defined in Eq. (1) and appearing in the rest of Sec. II), we can now use the simple cubic fits for the three densities to determine the competition between different phases. More specifically, Fig. 1 shows the behavior of the NP (normal polarized) equation of state. From the hollow points on the same figure we also have access to the energy of the SF (unpolarized superfluid) phase. Using these two dependences, we can find out the critical values of the relative fraction (, respectively) at the three densities we are studying. Above a critical relative fraction the system phase separates into a mixture of SF and NP, whereas below that value the neutron gas is normal, and therefore does not exhibit any of the widely studied features of a superfluid in the neutron star crust.
To determine we need the functional dependence of on both and , which is given by the interpolations for , and . Let us first focus on . In Eq. (2) we also need the derivative with respect to . We estimate this derivative by finite differences: and . We estimate the derivatives for used in Eq. 4 with a linear form in and . In practice, is small and doesn’t affect the value of and significantly. We can then solve the conditions for coexistence and find and . The numerical values show slight dependence on the interpolation method used. A simple cubic fit as a function of gives and fm. A three parameter fit in motivated by (1); (40) gives and fm. This gives some idea about the errors made in the interpolation: we estimate our error bar on to be for the three densities. This detailed procedure gives results very similar to a simpler procedure, namely constructing a tangent to the constant density curve for passing through the value of at the same density, which gives . The close agreement between the two approaches stems from the fact that is small — for , the coexistence curve is given exactly by the tangent construction (45). Therefore, we will use the tangent construction to find for and below.
This analysis requires that for the normal phase is thermodynamically stable. That is, the eigenvalues of the matrix are positive. For we can check and have checked explicitly that this is the case. Finally, we find at is given by where . Comparing with obtained from Ref. (47), we see that satisfying the third condition for coexistence.
In Fig. 2 we have plotted as a black solid line the results from Fig. 1; at this density the perturbative correction in the Hamiltonian mentioned in Sec. III is the smallest, leading to the highest degree of confidence in the accuracy of our microscopic results. We focus on large values on the axis for reasons that will soon become clear. We also show as red square points the unpolarized superfluid result from Ref. (47) (at ) and the homogeneous polarized superfluid results of Ref. (6) (at smaller ). In addition, we show a coexistence curve which is arrived at by a tangent construction from the unpolarized superfluid to the normal polarized curve when the axis is chosen as we have (45). This coexistence line meets the NP curve at (). As can be clearly seen in this Figure, the coexistence line lies below the homogeneous polarized superfluid Quantum Monte Carlo results, implying that these are not energetically favored (foreseeing this possibility, Ref. (6) included a critical assumption explicitly excluding the possibility of phase separation). This situation follows the behavior of ultracold fermionic gases at unitarity: (40); (41); (1) there the critical value of the relative fraction was closer to . The value we find is different for a variety of reasons: first, corresponds to , which is somewhat smaller than . More importantly, the neutronneutron interaction is characterized by a finite effective range and the presence of higher partial waves. The latter point explains why when we repeat this exercise for () and () the corresponding results (from tangent constructions) for the critical fraction are () and (), respectively. The trend appears to be that at larger density (stronger coupling) the critical fraction is moving toward 1, implying that farther down toward the core of the neutron star polarized neutron matter quickly becomes normal for the vast majority of possible polarizations.
The value of is relevant for neutron star phenomenology since whether is greater or less than determines the phase. This is seen more simply in the grand canonical ensemble where we can compare the pressure as a function of the chemical potentials. In Fig. 3 we plot the pressure as a function of for fixed MeV. For this value of the density of the superfluid phase at the critical is . For the pressure of the phase is larger and the system is unpolarised. Even if the magnetic field is not large enough to polarize the system, certain thermal observables, for example the specific heat, are sensitive to the difference between and .
Determining the critical relative fractions is also relevant: if a magnetic field is strong enough to polarize neutron matter sufficiently (), the gas ceases to be superfluid even at zero temperature. At such large magnetic fields, processes that rely on the presence of a bulk superfluid in the inner core of neutron stars, for example a recently proposed heatconduction mechanism which requires superfluid phonons (10), would no longer be operational. On a different note, a precise determination of these fractions (along with the rest of the equation of state) is also significant for Skyrme functional practitioners: for the region where the NP phase is the true ground state of the system, the results shown in Fig. 1 are a microscopic constraint to nuclear energydensity functional theory, which can help improve its predictions on neutronrich heavy nuclei.
V Summary & Conclusions
In summary, we have studied spinpolarized lowdensity neutron matter at many values of the polarization. Our Quantum Monte Carlo approach provides tight upper bounds which are expected to be quite close to the true ground stateenergy of the system. At verylow density the Hamiltonian becomes quite simple and we include its dominant wellknown terms. We have calculated the energy as a function of polarization using the AV4 interaction at three different densities.
One of our main results is a determination of the critical relative fraction values above which the system phase separates into an unpolarized superfluid and a polarized normal neutron gas. Above these values, a homogeneous polarized superfluid is found not to be energetically favored. Below these values, the neutron gas is completely normal. Both these facts can lead to astronomically important consequences. Importantly, the trend of our results seems to imply that at slightly larger densities would be very close to one: this would mean that even at zero temperature a small polarization would be enough to close the superfluid pairing gap, and thus easily lead to larger polarization. In a nonastrophysical context, it is conceivable that our results could be tested directly by using ultracold fermionic atom gases with unequal spin populations and a finite effective range. In cold atoms, Quantum Monte Carlo simulations of spinpolarized matter have in recent times been repeatedly used as input to computationally less demanding densityfunctional theory approaches. Similarly, we expect that the results presented in this work can also be used as input to selfconsistent meanfield models of nuclei.
A possible direction for future work, that would be interesting to both nuclear astrophysicists and Skyrme practitioners, would be a study of the behavior of polarized and unpolarized superfluid and normal neutron matter in an external periodic potential. Such a study of the static response of neutron matter could also in principle be guided (or guide) analogous research related to optical lattice experiments with cold atoms.
In the neutron star context, it will be useful to analyze the finite temperature phase diagram to understand how the specific heat is modified as a function of the chemical potential splitting, and how this modification affects transient phenomena in the neutron star crust.
Acknowledgements.
A.G. thanks George Bertsch, Kai Hebeler, and Achim Schwenk for stimulating discussions. R.S. acknowledges discussions with Jeremy Heyl. The authors thank Mark Alford and Sanjay Reddy for valuable comments and suggestions. The authors acknowledge the hospitality of the Perimeter Institute and the organizers of the Micra 2011 workshop where this work was initiated. This work was supported by DOE Grant No. DEFG0297ER41014 and the National Sciences and Engineering Research Council of Canada (NSERC). Computations were performed at the US DoE’s National Energy Research Scientific Computing Center (NERSC).References
 S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 80, 1215 (2008).
 S. Fantoni, A. Sarsa, and K. E. Schmidt, Phys. Rev. Lett. 87, 181101 (2001).
 I. Vidaña, A. Polls, and A. Ramos, Phys. Rev. C 65, 035804 (2002).
 A. Rios, A. Polls, and I.Vidaña, Phys. Rev. C 71, 055802 (2005).
 M. G. Alford, J. Berges, K. Rajagopal, Nucl. Phys. B558, 219242 (1999).
 A. Gezerlis, Phys. Rev. C 83, 065801 (2011).
 C. J. Pethick, N. Chamel, S. Reddy, Prog. Theor. Phys. Suppl. 186, 916 (2010).
 V. Cirigliano, S. Reddy, and R. Sharma, Phys. Rev. C 84, 045809 (2011).
 E. F. Brown and A. Cumming, Astrophys. J. 698, 1020 (2009).
 D. N. Aguilera, V. Cirigliano, J. A. Pons, S. Reddy, and R. Sharma, Phys. Rev. Lett. 102, 091101 (2009).
 M. G. Alford, S. Reddy, and K. Schwenzer, arXiv:1110.6213 (2011).
 J. Carlson, J. Morales, Jr., V. R. Pandharipande, and D. G. Ravenhall, Phys. Rev. C 68, 025802 (2003).
 A. Gezerlis and J. Carlson, Phys. Rev. C 77, 032801(R) (2008).
 P. F. Bedaque, H. Caldas, and G. Rupak, Phys. Rev. Lett. 91, 247002 (2003).
 R. Sharma and S. Reddy, Phys. Rev. A78, 063609 (2008).
 A. Bulgac, M. M. Forbes, and A. Schwenk, Phys. Rev. Lett. 97, 020402 (2006).
 A. Bulgac and M. M. Forbes, Phys. Rev. Lett. 101, 215301 (2008).
 M. Bender, P.H. Heenen, and P.G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).
 B. Friedman and V. R. Pandharipande, Nucl. Phys. A361, 502 (1981).
 A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Phys. Rev. C 58, 1804 (1998).
 A. Schwenk and C. J. Pethick, Phys. Rev. Lett. 95, 160401 (2005).
 J. Margueron, E. van Dalen and C. Fuchs, Phys. Rev. C 76, 034309 (2007).
 M. Baldo and C. Maieron, Phys. Rev. C 77, 015801 (2008).
 E. Epelbaum, H. Krebs, D. Lee, and U. G. Meissner, Eur. Phys. J. A 40, 199 (2009).
 E. Epelbaum, H.W. Hammer, and U. G. Meissner, Rev. Mod. Phys. 81, 1773 (2009).
 T. Abe and R. Seki, Phys. Rev. C 79, 054002 (2009).
 S. Gandolfi, A. Yu. Illarionov, F. Pederiva, K. E. Schmidt, and S. Fantoni, Phys. Rev. C 80, 045802 (2009).
 A. Rios, A. Polls, and I. Vidaña, Phys. Rev. C 79, 025802 (2009).
 K. Hebeler and A. Schwenk, Phys. Rev. C 82, 014314 (2010).
 N. Chamel, S. Goriely, and J.M. Pearson, Nucl. Phys. A812, 72 (2008).
 M. M. Forbes, A. Gezerlis, K. Hebeler, T. Lesinski, and A. Schwenk, under preparation (2011).
 L. Luo and J. E. Thomas, J. Low Temp. Phys. 154, 1 (2009).
 N. Navon, S. Nascimbène, F. Chevy, and C. Salomon, Science 328, 729 (2010).
 J. Carlson and S. Reddy, Phys. Rev. Lett. 100, 150403 (2008).
 A. Schirotzek, Y. I. Shin, C. H. Schunck, and W. Ketterle, Phys. Rev. Lett. 101, 140403 (2008).
 Y. Shin, C. H. Schunck, A. Schirotzek, and W. Ketterle, Nature 451, 689 (2008).
 J. Carlson, S. Y. Chang, V. R. Pandharipande, and K. E. Schmidt, Phys. Rev. Lett. 91, 050401 (2003).
 G. E. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini, Phys. Rev. Lett. 93, 200404 (2004).
 J. Carlson and S. Reddy, Phys. Rev. Lett. 95, 060401 (2005).
 C. Lobo, A. Recati, S. Giorgini, and S. Stringari, Phys. Rev. Lett. 97, 200403 (2006).
 S. Pilati and S. Giorgini, Phys. Rev. Lett. 100, 030401 (2008).
 A. Gezerlis, S. Gandolfi, K. E. Schmidt, and J. Carlson, Phys. Rev. Lett. 103, 060403 (2009).
 B. Marcelis, B. Verhaar, and S. Kokkelmans, Phys. Rev. Lett. 100, 153201 (2008).
 T. D. Cohen, Phys. Rev. Lett. 95, 120403 (2005).
 A. Bulgac and M. M. Forbes, Phys. Rev. A 75, 031605 (2007).
 R. B. Wiringa and S. C. Pieper, Phys. Rev. Lett. 89, 182501 (2002).
 A. Gezerlis and J. Carlson, Phys. Rev. C 81, 025803 (2010).
 M. M. Forbes, S. Gandolfi, and A. Gezerlis, Phys. Rev. Lett. 106, 235303 (2011).