Ab initio and phenomenological studies of the static response of neutron matter
aa r X i v : . [ nu c l - t h ] A p r Ab initio and phenomenological studies of the static response of neutron matter
Mateusz Buraczynski and Alexandros Gezerlis
Department of Physics, University of Guelph, Guelph, ON N1G 2W1, Canada
We investigate the problem of periodically modulated strongly interacting neutron matter. Wecarry out ab initio non-perturbative auxiliary-field diffusion Monte Carlo calculations using anexternal sinusoidal potential in addition to phenomenological two- and three-nucleon interactions.Several choices for the wave function ansatz are explored and special care is taken to extrapolatefinite-sized results to the thermodynamic limit. We perform calculations at various densities aswell as at different strengths and periodicities of the one-body potential. Our microscopic resultsare then used to constrain the isovector term from energy-density functional theories of nuclei atmany different densities, while making sure to separate isovector contributions from bulk properties.Lastly, we use our results to extract the static density-density linear response function of neutronmatter at different densities. Our findings provide insights into inhomogeneous neutron matter andare related to the physics of neutron-star crusts and neutron-rich nuclei.
I. INTRODUCTION
Neutron matter is integral to the study of neutron starsand neutron-rich nuclei [1]. The equation of state (EOS)of neutron matter has been studied extensively using abinitio approaches [2–7]. While neutron matter does notobtain in nature in pure form, its EOS is closely con-nected to that of physical systems. More specifically,it is direct input into the Einstein’s field equations (typi-cally cast as the Tolman-Oppenheimer-Volkov equations)that lead to basic observable properties of a neutron star.On the other hand, the neutron-matter EOS is also con-nected to nuclei via the use of nuclear energy densityfunctionals (EDFs). EDFs take on a variety of forms [8–11] and are typically fit to empirical data such as nuclearmasses and radii. Other constraints may include the EOSof neutron matter [9, 12–17], the neutron-matter pairinggap [18], the energy of a neutron impurity [19, 20], orthe properties of a neutron drop [21–25]. The EDF ap-proach is especially useful as it allows for predictions tobe made across the nuclear chart. Investigating neutronmatter is an excellent opportunity for benchmarking bothphenomenological [26–30] and chiral [31–40] nuclear in-teractions and many-body methods.While homogeneous matter is an intriguing system, itis not fully representative of either finite systems or as-trophysical settings, since nuclei and neutron-star matterare inhomogeneous systems. The matter in a neutron-star crust is rich with unbound neutrons and also con-tains a lattice of nuclei. Focussing only on the unboundneutrons, these experience the interaction with the latticeas a periodic modulation. Similarly, nuclei are finite sys-tems: their density eventually drops off as one goes far-ther away from the center of the nucleus. Thus, as a firstapproximation, one can study the effects of a one-bodyexternal periodic potential on pure infinite neutron mat-ter. This system, periodically modulated neutron mat-ter, directly mimics the situation in a neutron star; also,if used as an input constraint to EDFs, it can inform usabout the physics of neutron-rich nuclei. This problemof an external periodic modulation is known as the static response of neutron matter: it amounts to a comparisonbetween externally modulated and unmodulated infiniteneutron systems. This problem has been tackled usinga variety of approximations in the literature [41–48], seeRef. [49] for a recent review. The static response prob-lem has also received a lot of attention in other areas ofphysics [50]. This includes pioneering quantum MonteCarlo calculations for strongly correlated systems suchas liquid He at zero temperature and pressure [51] andthe three-dimensional electron gas [52].In a recent Letter [53] we reported on the first ab ini-tio calculation of the static response of neutron mat-ter. This included a calculation of the linear density-density static response function of neutron matter fol-lowing from a Quantum Monte Carlo (QMC) T = 0 ap-proach. Specifically, this involved the Variational MonteCarlo (VMC) and auxiliary-field diffusion Monte Carlo(AFDMC) methods: these are non-perturbative and ac-curate methods for computing properties of many-bodynuclear systems. Reference [53] also mentioned the roleof finite-size effects, the wave-function ansatz, and con-clusions that can be drawn on the isovector gradient co-efficient of EDFs.In the present article, we expand on Ref. [53], by dis-cussing finite-size effects and the wave-function ansatzin more detail. In addition to this, we have carriedout new QMC calculations for the periodically mod-ulated system at many different densities, extractingfrom there EDF parameters for several different Skyrmeparametrizations. We have also carried out several newcalculations of the modulated system for several differentstrengths and periodicities of the external one-body po-tential: this includes a study of the lower-density regime,as well as separate simulations (at both low and inter-mediate densities) with and without three-nucleon in-teractions. These results provide a more detailed un-derstanding of different effects and, since they are non-perturbative and accurate, they provide “synthetic data”that can be used in (or compared to) other approaches.We make some first steps in this direction by examiningthe response function coming from selected EDFs.We begin with some background on the methods em-ployed (section II), mainly establishing our notation inwhat follows. We also provide a more extensive discus-sion of the theory of static response (section II C), ex-plaining what the density-density response is, and how itcan be extracted from our calculations. We then proceedto study finite-size effects in the non-interacting problemin detail (section III), before turning to the QMC andEDF results for the interacting and periodically modu-lated problem (section IV). II. METHODSA. Auxiliary-field diffusion Monte Carlo
The Hamiltonian of the full interacting problem westudy is made up of a non-relativistic kinetic term, atwo-nucleon (NN) interaction, a three-nucleon (NNN) in-teraction, and a one body potential:ˆ H = − ~ m X i ∇ i + X i Another method for computing the energy of this sys-tem comes from density functional theory (DFT). In thiscase, the nuclear force takes a phenomenological form ofthe Skyrme type. The many-body wave function takesthe form of a Slater determinant of single-particle orbitals ψ i ( r ). The ground-state energy is given by: E = Z H ( r ) d r, (5)where H is the energy density energy functional (EDF)[8]: H = ~ m τ + 2 v q cos( q · r ) n + E Sk (6)and n ( r ) = X i [ ψ i ( r )] τ ( r ) = X i [ ∇ ψ i ( r )] (7)are the nucleon and kinetic energy densities. The Skyrmeinteraction terms of the EDF in the isospin representa-tion are: E Sk = X T =0 , h ( C n,aT + C n,bT n σ ) n T + C ∆ nT ( ∇ n T ) + C τT n T τ T i (8)For pure neutron matter n = n = n and the same istrue for τ . We have done calculations for the SLy4, SLy7,and SkM* parametrizations. C ∆ n is known as the isovec-tor gradient term. C ∆ n = − , − , and − 17 MeV fm inSLy4, SLy7, and SkM* respectively [9]. In what follows,we will adjust this parameter based on our AFDMC re-sults. Eq. (5) provides us with a method to approximatethe energy using the density functions. This is calledthe local-density approximation. Our method of approx-imating the energy is similar to the VMC optimization ofthe trial wave function. The ψ i ’s are the orbitals of thenon-interacting Hamiltonian with our external potential(Mathieu functions). We minimize the right hand side ofEq. (5) to get the local-density approximation energies. C. Static-response theory An objective of this work is to compute the lineardensity-density static response function of neutron mat-ter. This gives a quantitative (up to first order) descrip-tion of the effect of an external perturbation on the phys-ical properties of a homogeneous neutron gas. Let ˆ H denote the unperturbed Hamiltonian. This is Eq. (1)without the v ext = P i v ( r i ) term. The 0 subscript isour notation for the unperturbed system. (Note thatthis is different from the, standard, n notation used inEq. (8) in the isospin representation of the EDF.) Theground-state density of the system is a functional of theexternal potential: n v ( r ) = n ( r , [ v ]) The density-densityresponse functions are defined as the functional deriva-tives of density with respect to v [56]: n v ( r ) = n + ∞ X k =1 k ! Z d r . . . d r k χ ( k ) ( r − r , . . . , r k − r ) v ( r ) . . . v ( r k ) , (9)where the χ ( k ) ’s are the response functions. χ (1) ( r ) is thelinear density-density response function. Likewise, theground state energy can also be expressed as a functionalof v : E v = E ([ v ]). The energy can be expressed as: E v = E + Z dλ Z d rn ( r , [ λv ]) v ( r ) (10)This follows from first-order perturbation theory: δE v /δv ( r ) = n v ( r ). This is easy to see if the interac-tion term is cast as R d r ˆ n ( r ) v ( r ), where the one-bodydensity operator is ˆ n = P i δ ( r − r i ).The energy and density can be expressed in termsof the Fourier components of the potential v ( r ) = P q v q exp( i q · r ) and the Fourier transforms of the re-sponse functions with respect to their spatial arguments χ ( k ) ( q , . . . , q k ): n v ( r ) = n + ∞ X k =1 k ! X q ,..., q k χ ( k ) ( q , . . . , q k ) v q . . . v q k × exp[ i ( q + . . . + q k ) · r ] E v N = E N + v + 1 n ×× ∞ X k =1 k + 1)! X q + q + ... + q k =0 χ ( k ) ( q , . . . , q k ) v q v q . . . v q k (11)For the one-body external potential v ( r ) = 2 v q cos ( q · r )the density is given by: n v ( r ) = n + 2 n q cos( q · r ) n q = χ ( q ) v q + χ ( q , q , − q )2 v q + . . . (12)The change in the density n v ( r ) − n depends only onodd powers of v q . The change in energy is: E v N = E N + χ ( q ) n v q + χ ( q , q , − q )4 n v q + . . . (13)The energy change only depends on even powers of v q .If the energy per particle is known at several differentvalues of v q then Eq. (13) gives a method to calculatelower-order response functions by fitting to a polynomialof even powers. The coefficient of the quadratic termgives the linear density-density response function. Thecoefficient of the quartic term is very small in our calcu-lations, on the order of 10 − MeV − or smaller. Higher-order fits, with more points, are required to reliably ex-tract the third-order response function.The response of a non-interacting Fermi gas can becomputed analytically. It is given by the Lindhard func-tion [57]: χ L = − mq F π ~ " q F q − (cid:18) q q F (cid:19) ! ln (cid:12)(cid:12)(cid:12)(cid:12) q + 2 q F q − q F (cid:12)(cid:12)(cid:12)(cid:12) (14)We compare our results in later sections with this re-sponse. III. NON-INTERACTING PROBLEM Many of the concepts needed for an understanding ofinteracting particles require an understanding of a muchsimpler problem: the 3D non-interacting free-Fermi gas.For a finite system of N particles it is standard to restrictpositions to a cubic box of volume V = L and imposeperiodic boundary conditions on the wave function, whenan extended system is the end goal of the study. At T = 0particles occupy states corresponding to the lowest avail-able energy levels. A state is identified by its momentum 10 100 1000 N | E F G ( N ) / E F G - | 20 40 60 80 100 N E F G ( N ) [ E F G ] FIG. 1. FS dependence of the non-interacting free-Fermi gas.FS effects go to zero in the thermodynamic limit. The mini-mum at N = 67 motivates a study at the closed shell N = 66.The inset plots the same relationship on a linear scale. wave-vector k = (2 π/L )( n x , n y , n z ) where n x , n y , n z areintegers. For a spin-1/2 system a maximum of two par-ticles with different spin-projection can occupy the samewave vector, due to the Pauli exclusion principle. A par-ticle placed in the state with wave-vector k has energy E = ~ k / m and occupies the single particle orbital ψ k ( r ) = e i k · r / √ V . Closed shell configurations where theenergy level populations are filled to capacity occur at N = { , , , , , ... } . It is preferable to workwith closed shells to avoid any ambiguity in determiningwhich states are occupied.Since a neutron star is a macroscopic system we areparticularly interested in the thermodynamic limit (TL)where N → ∞ , V → ∞ and n = N/V is constant. Inthe TL the number density n is related to k F , the max-imal wave vector magnitude by k F = (3 π n ) / . Theenergy per particle is given by E F G = (3 / ~ k F / m =(3 / E F . Differences in properties between the ther-modynamic limit and a finite number of particles arecalled finite-size (FS) effects. FS effects go to zero inthe thermodynamic limit. They are also largest at small N . This can be seen in Fig. 1 which plots FS effectsin the energy per particle versus particle number for thenon-interacting free-Fermi gas. This result is density-independent. As mentioned earlier, we are primarily in-terested in shell closures. These appear at the cusps inthe inset of Fig. 1. There is a minimum in FS effects oc-curring at 67 particles. This leads us to study the closedshell at 66 particles.We now extend the above problem to include our ex-ternal potential. The Hamiltonian is:ˆ H = − ~ m X i ∇ i + v ext , (15)where v ext = P i v ( r i ) and v ( r i ) = 2 v q cos ( q · r i ), as be-fore. The orbitals are given by Mathieu functions and the N E N I ( N ) [ E F G ] FIG. 2. FS dependence of the non-interacting Fermi gas inthe presence of a one-body potential of fixed strength 2 v q =0 . E F and periodicity q = 4 π/L . n is fixed at 0 . − . energies by the corresponding characteristic values. Sim-ilarly to the free gas, intensive properties converge in theTL. We calculate the energy per particle of the perturbedgas in the same manner as the free gas. Doing so demon-strates convergence in the TL. This is shown in Fig. 2where the line gives energy per particle versus particlenumber. The density is fixed at 0 . − and q is set to4 π/L so that 2 periods of the potential fit inside the box.The amplitude of the potential is 0 . E F = 21 . 363 MeV.The curve in Fig. 2 shows a convergence near 0.8: thisis not 1, as could be expected from the presence of theexternal potential.We use the energy per particle of the non-interactingsystem to handle FS effects in the interacting problem(see also Ref. [58]). Our goal is to describe an extendedneutron system. To accomplish this, we keep the po- N | E N I ( N ) / E N I ( ∞ )- | N E N I ( N ) [ E N I ( ∞ )] FIG. 3. FS dependence of the non-interacting Fermi gas inthe presence of a one-body potential of fixed strength 2 v q =0 . E F and periodicity q = 1 . − . n is fixed at 0 . − . q/q F - χ ( q ) / n [ M e V - ] Lindhard function66 particles66000 particles FIG. 4. Static-response function of the non-interacting free-Fermi gas at a density of 0 . − . The squares and circlesare for 66 and 66000 particles respectively. The line is theLindhard function describing the response in the TL. tential fixed and compute the energy per particle versusparticle number. Keeping in mind translational invari-ance, we require an integer number of periods in the box.Thus we choose to perform calculations for a discrete setof particle numbers. This is displayed in Fig. 3 wherethe density and amplitude are 0 . − and 21 . 363 MeVrespectively. q is fixed so that two periods fit in the boxat N = 66. Of course, the energy per particle convergesin the TL. FS effects in the energy per particle E I ( N ) ofthe interacting system will be handled by extrapolatingto the TL: E I ( ∞ ) = E I ( N ) − E NI ( N ) + E NI ( ∞ ) (16)We tested Eq. (16) by applying it to energy calculationsof homogeneous neutron matter. This was done for theSLy4 energies of 66 neutrons at various densities. Theresults are shown in Table I. They agree with SLy4 en-ergies of homogeneous neutron matter to within 0.5%.This boosts our confidence in Eq. (16).We further highlight the importance of FS correctionsby comparing calculations of the response function tothe analytically known response in the TL. The responsefunction for 66000 particles at a density of 0 . − (cir-cles in Fig. 4) matches the Lindhard function (solid line).This makes sense since 66000 particles is practically inthe TL, as per Fig. 3. The response for 66 particles atthis density (squares) does not match the Lindhard func-tion except at large q . This stresses that 66 particles isnot in the TL so it is important that FS effects be han-dled in order to study infinite neutron matter.(Note thatthe FS handling in Ref. [53] suffered from a numericalerror in the, near-trivial, calculation for E NI (66); a sim-ilar error was present in E NI (66000) but was immaterialthere. This error is corrected here.) We note that thebehavior exhibited by the 66-particle results in Fig. 4 isobserved at other densities also: there is a dip as the q is n (fm -3 ) E / N [ M e V ] no V ext q = 0.5 E F (PW)2v q = 0.5 E F (opt.) FIG. 5. AFDMC neutron-matter energy per particle for 66particles as a function of density using only NN interactions.Squares denote the case without a one-body potential; dia-monds a one-body potential of fixed strength 2 v q = 0 . E F ,periodicity q = 4 π/L , and plane-wave single-particle orbitals;circles a one-body potential of fixed strength 2 v q = 0 . E F ,periodicity q = 4 π/L , and optimized single-particle orbitals(opt.). lowered, before the response goes back up for the lowest- q point. IV. INTERACTING PROBLEMA. Equation of state We computed the equation of state of 66 neutrons bothwith and without an external potential. We first presentresults that do not include NNN interactions. The NNinteractions are given by the Argonne v8’ potential [59].Calculations were performed for densities in the rangeof 0 . 02 fm − to 0 . 12 fm − since these are similar to thedensities found in the crust and outer-core of a neutronstar. The nodal structure of the trial wave function is im-portant, so the single-particle orbitals used in the Slaterdeterminant must be chosen carefully. We use the so- TABLE I. SLy4 energies of a free-Fermi gas of neutrons com-puted using the local density approximation. Eq. (16) wasused to extrapolate to the TL and compare to SLy4 energiesof homogeneous neutron matter. n (fm − ) E I (66) (MeV) E I ( ∞ ) (MeV) TL (MeV)0 . 04 7 . 15 7 . 22 7 . . 06 8 . 62 8 . 71 8 . . 08 9 . 99 10 . 10 10 . . 10 11 . 40 11 . 53 11 . 10 12 14 16 18 20 22 24 26 28 30 β (MeV) V M C E / N ( M e V ) FIG. 6. Neutron-matter VMC energy per particle asa function of the variational parameter β using NNinteractions. n = 0 . − , 2 v q = 0 . E F = 21 . 363 MeV, and q = 4 π/L . lutions of the one-body problem with the same externalpotential and no interactions. For the case of no externalpotential these are the plane-waves of the non-interactingfree-Fermi gas. The AFDMC results for these are thesquares in Fig. 5. The energy increases with increasingdensity. The energies agree with known values [60]. Wecomputed energies using plane-wave orbitals for a onebody potential of strength 2 v q = 0 . E F and two periodsof the potential in the box (diamonds). This yields lowerenergies than the unperturbed problem. However, this isnot good enough since we have not optimized the trialwave-function. The energies for optimized single-particleorbitals (circles) are up to 1 MeV different from theplane-wave results. Overall, the AFDMC energies withan external potential are several MeV smaller than thosewithout an external potential. This reflects the parti-cles’ tendency to stay away from the repulsive regionsof the potential and collect in the wells of the potential.Note that the amplitude of the external potential usedis density dependent. Also, the period of the potentialdecreases with increasing density. This is consistent withwhat happens in a neutron-star crust where the latticespacing decreases with increasing density.We now describe the optimization procedure. Firstly,Mathieu functions yield lower VMC energies than plane-waves in the Slater determinant. Since this is a varia-tional optimization the goal is to minimize the VMC en-ergy. We do this using a variational parameter β wherethe solutions to the one-body non-interacting problemwith external potential v ( r ) = β cos( q · r ) [51] are usedin the Slater determinant. For the case β = 2 v q theseare the orbitals of the one-body non-interacting prob-lem with the same external potential as the system un-der study. Consider a density of 0 . − for the casewhere 2 v q = 0 . E F and two periods of the potential fit inthe box. The minimum VMC energy that we simulated q /E F E / N [ M e V ] no FS fixwith FS fix FIG. 7. Neutron-matter energy per particle as a function ofthe one-body potential strength using NN+NNN interactionsand AFDMC. n = 0 . − and q = 4 π/L . Circles corre-spond to energies prior to handling FS effects and squares toenergies extrapolated to the TL. occurs at β = 18 MeV (Fig. 6). The AFDMC energyfor β = 18 MeV is 8.19 MeV. The AFDMC energy for β = 2 v q is 8.15 MeV. Thus, the difference in energy dueto the β optimization procedure is much smaller than thedifference due to using plane-waves vs Mathieu functions.Therefore, it is safe to set β = 2 v q for practical purposes.We are primarily interested in the static response ofneutron matter. This requires us to perform calculations q /E F E / N [ M e V ] AFDMCSLy4SLy4 (mod.) FIG. 8. Neutron-matter energy per particle for 66 particles asa function of the one-body potential strength at a density of n = 0 . 10 fm − , using NN+NNN interactions and a one-bodypotential periodicity q = 4 π/L . Circles denote AFDMC re-sults, and the solid line follows from the SLy4 energy-densityfunctional. The dashed line corresponds to SLy4 results witha modified isovector gradient term. n [fm -3 ] E / N [ M e V ] , v q = AFDMCSLy4 FIG. 9. TL AFDMC and TL SLy4 results in the absence ofa one-body potential. across a set of strengths of the external potential. It isimportant that the strengths are large enough so that theenergy is statistically different from homogeneous neu-tron matter. At the same time the filling of the singleparticle orbitals changes at larger v q in comparison to thefree non-interacting orbital filling. We performed simu-lations for various orbital fillings and found no signifi-cant change in the energy. We did this for 2 v q = 0 . E F at n = 0 . 08 fm − for two different fillings and foundAFDMC energies of 7 . 072 and 7 . 062 MeV. The samewas done for 2 v q = 0 . E F at n = 0 . 04 fm − yielding1 . 818 and 1 . 822 MeV. We decided to calculate energiesfor 2 v q = 0 . , . , . , . , . E F . NNN interactionswere included in these simulations. Specifically, the Ur-bana IX potential was used, as is appropriate for neutron-rich matter [54]. Turning up the potential strength re-sults in a decrease in energy. This is seen in the circlesin Fig. 7 for 66 particles at a density of 0 . − and twoperiods of the potential in the box. We have also appliedEq. (16) to extrapolate to the TL (squares). We extractthe response function from such results. We present thoseresults in the next section. B. Constraining the isovector gradient coefficient We have also used the response of neutron matterto constrain energy density functionals. The energy inEq. (5) is minimized with respect to the variational pa-rameter β and different orbital fillings. This is done fora range of potential strengths. Fig. 8 displays the resultsof this procedure for SLy4 (solid line) for a 66 particlesystem at a density of 0 . 10 fm − and two periods of thepotential in the box. The 66 particle AFDMC NN+NNNenergies (circles) are more repulsive than the SLy4 ener-gies. The separation between the AFDMC and SLy4 en-ergies is largest at v q = 0. The separation decreases as v q q /E F E / N [ M e V ] AFDMCSLy4SLy4 (mod.) FIG. 10. Neutron-matter energy per particle for 66 particlesas a function of the one-body potential strength at a density of n = 0 . 04 fm − , using NN+NNN interactions and a one-bodypotential periodicity q = 4 π/L . Circles denote AFDMC re-sults, and the solid line follows from the SLy4 energy-densityfunctional. The dashed line corresponds to SLy4 results witha modified isovector gradient term. increases. ∇ n = 0 at homogeneity so the isovector gradi-ent term does not contribute to the v q = 0 energy. Thusthe difference between the AFDMC and SLy4 energiesat homogeneity is due to the bulk energy of the system.This homogeneous mismatch in energy must be respectedin fitting the EDF to AFDMC results. The isovector gra-dient term has the effect of bringing the SLy4 energiescloser to the AFDMC energies at larger v q . Thus, our fit-ting of the isovector gradient term maintains the v q = 0difference between the EDF and QMC energies. We fitto low strengths 2 v q = 0 . 25 and 0 . F , in order to ensurethat the density perturbation magnitude is not compa-rable to the unperturbed density. We found a modifiedisovector gradient term of C ∆ n = − 29 MeV fm at thisdensity (dashed line in Fig. 8).The homogeneous energy difference between AFDMCand SLy4 impacts how the isovector gradient term shouldbe modified. Above n = 0 . 06 fm − the homogeneousAFDMC energy (circles in Fig. 9) is more repulsive thanthe homogeneous SLy4 energy (solid line in Fig. 9). Be-low this density the AFDMC is more attractive thanthe SLy4 energy. We see this at n = 0 . 04 fm − wherethe AFDMC NN+NNN energies (circles in Fig. 10) aresmaller than SLy4 (solid line in Fig. 10) for small enough v q . The separation between AFDMC and SLy4 decreaseswith increasing v q at both densities larger and smallerthan 0 . 06 fm − . This means that the fitted isovectorgradient term is density dependent. This result is alsofound using the density-matrix expansion [61, 62]. For n > . 06 fm − the isovector gradient fit must make theSLy4 energies more attractive if they are to be equidis-tant from the AFDMC energies. This requires a decrease n [fm -3 ] -40-200204060 i s ov ec t o r g r a d i e n t c o e ff . [ M e v f m ] SLy4 modifiedSLy7 modifiedSkM* modified FIG. 11. Isovector gradient coefficients for the modified SLy4,SLy7, and SkM* Skyrme potentials. All calculations weredone with two periods of the potential in the box. Coefficientswere modified to make Skyrme responses match the QMCresponses. in the isovector gradient term. For n < . 06 fm − an in-crease in the isovector gradient term is required. At adensity of 0 . 04 fm − we find a modified isovector term of C ∆ n = 9 MeV fm (dashed line in Fig. 10). Since we fitthe isovector term to low strengths, the modified SLy4still approaches the AFDMC results at some large v q .We have carried out calculations such as those in Fig. 8and Fig. 10 for several other densities. They exhibit qual-itatively the same trends as discussed above. We thenused these AFDMC results to constrain the isovector co-efficient for several functionals. Specifically, Fig. 11 liststhe isovector gradient term of modified SLy4, SLy7, andSkM*. All of these were done using two periods of the po-tential in the box. The errors were determined by fittingto different strengths and examining the spread of themodified isovector gradient term. Note that the densitydependent energy versus v q behaviour described abovedoes not hold for SkM*. Nevertheless, our fitting pre-scription yields the same density dependence in the mod-ification of the isovector term for all three parameteriza-tions. We see a decrease is required at large densities. Atlow densities there is an increase in the isovector term.The isovector term is least modified at n = 0 . 06 fm − ,where in the homogeneous case the AFDMC and SLy4results agree reasonably well, as per Fig. 9.Note that the isovector coefficient fits discussed abovewere all carried out by focusing on the EDF and QMCresults for L = 2 d , i.e., using two periods of the potentialin the box. As we will see in the following subsection,our attempt to make the QMC and EDF results equidis-tant essentially amounts to trying to match (not QMCenergies to EDF energies, but) the EDF response func-tion to the QMC response function. In other words, weare attempting to modify the SLy4 response from Fig. 12 q/q F - χ ( q ) / n [ M e V - ] Lindhard functionSLy4 TL SLy4SLy4 modified FIG. 12. Static-response functions of neutron-matter at adensity of 0 . 10 fm − . The circles are the SLy4 response ex-trapolated to the TL limit. The diamonds are for the modifiedSLy4 with C ∆ n = − 29 MeV fm extrapolated to the TL. Theresponse was extracted by fitting to 2 v q = 0 . , . , . , . , and 0 . E F . The dashed curve is the SLy4 response pro-duced in the TL [63]. The solid line is the Lindhard functiondescribing the response of a non-interacting Fermi gas. below to match that in Fig. 13 below. All this, for thespecific case of L = 2 d . As the difference between SLy4and modified SLy4 in Fig. 12 shows, for L = 2 d we wouldneed a more attractive modification, whereas for L = d we would need a more repulsive one (and for L = 3 d wewould need a modification that is more attractive thanthat for L = 2 d ). Thus, the optimal thing to do is to tryto optimize results for as many periodicities as possiblesimultaneously. We have done this at the two densitiesof n = 0 . 10 fm − and n = 0 . 04 fm − , where we have pro-duced AFDMC results for many different periodicities,as discussed below. C. Response functions Using calculations like those discussed above, we ex-tracted linear density-density response functions at both0 . 04 fm − and 0 . 10 fm − for AFDMC, SLy4 and modifiedSLy4 results. (To do this, we fit to even-powered poly-nomials up to fourth order in Eq. (13).) Since we arestudying neutron matter we use Eq. (16) to extrapolateenergies to the TL. It was previously mentioned that weonly consider q such that an integer number of periods ofthe potential fit in the box. We have performed simula-tions for q = 2, 4, 6, 8, 12, 16, 20 times π/L correspondingto 1, 2, 3, 4, 6, 8, and 10 periods of the potential insidethe box.We have also taken advantage of the compressibilitysum rule: this provides us with a way to calculate χ (0)starting from the energy per particle as a function ofdensity of the unperturbed system:1 χ (0) = − ∂ ( n E/N ) ∂n (17)We used Eq. (17) to compute χ (0) for the various re-sponse functions that we extracted and checked for con-sistency with our modulated results.We first present results at n = 0 . 10 fm − . We havefound that the SLy4 response (circles in Fig. 12) doesnot match the Lindhard function (solid line in Fig. 12)although there are similarities. Both responses have afinite χ (0) and go to 0 at large q . We compare to theSLy4 response function in the TL [63] (dashed line inFig. 12). Our response agrees with it pretty well for allexcept the smallest q values. The compressibility sumrule gives a value of − χ (0) /n = 0 . 057 MeV − for SLy4which matches the TL SLy4. We are also interested inthe modified SLy4 response (diamonds in Fig. 12). Thisresponse is similar in shape but larger in magnitude thanthe SLy4 response we extracted. This makes sense sincethe modified SLy4 is more attractive than SLy4. Thecompressibility sum rule gives the same χ (0) for modifiedSLy4 as SLy4 since the unperturbed system is indepen-dent of the gradient term.In Fig. 13 we show updated AFDMC NN+NNN results(circles): note that these include the corrected FS han-dling and therefore differ (in the lowest- q response value)from the circles in Fig. 3 of Ref. [53]. The diamondsin Fig. 13 show the AFDMC NN response function inthe TL. Our results in the TL are similar in shape tothe Lindhard function (line in Fig. 13). At larger q the q/q F - χ ( q ) / n [ M e V - ] Lindhard functionAV8+UIX AV8 Compressibility sum rule (AV8+UIX) FIG. 13. Static-response functions of neutron-matter at adensity of 0 . 10 fm − . Produced using AFDMC results. Cir-cles are with NN+NNN interactions extrapolated to the TLlimit. Diamonds are for NN interactions extrapolated tothe TL. The response was extracted by fitting to: 2 v q =0 . , . , . , . , and 0 . E F . The square is the responseat q = 0 predicted by the compressibility sum rule for theNN+NNN case. The curve is the Lindhard function describ-ing the response of a non-interacting Fermi gas. q/q F - χ ( q ) / n [ M e V - ] Lindhard functionAV8+UIXAV8 FIG. 14. Static-response function of neutron-matter at a den-sity of 0 . 04 fm − . The circles are with NN+NNN interactionsextrapolated to the TL. Diamonds are for NN interactions ex-trapolated to the TL. The AFDMC responses were extractedby fitting to: 2 v q = 0 . , . , . , and 0 . E F . The solid lineis the Lindhard function describing the response of a non-interacting Fermi gas. response goes to zero and matches the Lindhard func-tion. At smaller q the NN+NNN response is smallerthan the Lindhard function. The NN response is largerthan the Lindhard function at the lowest q value, butother than that it is very similar to the NN+NNN re-sponse. The compressibility sum rule gives a value of − χ (0) /n = 0 . 043 MeV − for neutron matter with ourNN+NNN interactions (square in Fig. 13) and a value of0 . 089 MeV − for NN interactions only. These are largerthan the corresponding value of 0 . 035 MeV − for theLindhard function. Note that Fermi liquid theory yields − χ ( q ) /n ≈ . 035 MeV − at n = 0 . 10 fm − [41, 64].We now examine some of the responses at a den-sity of 0 . 04 fm − . We have found that the AFDMC re-sults do not follow the Lindhard function (solid line inFig. 14) as well as the 0 . 10 fm − AFDMC results do.Both the AFDMC NN response (diamonds in Fig. 14)and the AFDMC NN+NNN response (circles in Fig. 14)are larger than the Lindhard function at small q . Thecompressibility sum rule gives a value of − χ (0) /n =0 . 14 MeV − for neutron matter with our NN+NNN inter-actions and 0 . 19 MeV − for NN interactions only. Theseare larger than the − χ (0) /n = 0 . 065 MeV − of the Lind-hard function. Fermi liquid theory yields − χ ( q ) /n ≈ . 083 MeV − at n = 0 . 04 fm − [41, 64].For all results it was found that the response functiongoes to zero as q goes to infinity and χ (0) is finite. Inaddition, the response functions extracted from AFDMCand the Lindhard functions show more similarity to oneanother than to the SLy4 responses. Overall, the SLy4responses are narrower and steeper than the other re-sponses. It is also interesting to contrast these to theresponse functions of He [51] and the 3D electron gas0 n [fm -3 ] - χ ( q ≈ q F ) / n [ M e V - ] Lindhard functionAFDMC FIG. 15. Static-response function of neutron-matter acrossseveral densities. The circles correspond to AFDMC resultswith NN+NNN interactions (and are extrapolated to the TL).The solid line corresponds to the Lindhard function describingthe response of a non-interacting Fermi gas. [52] both of which have χ (0) = 0.Similarly to what was shown in Figs. 13 & 14, one canextract the response function of neutron matter for otherdensities. We have carried out precisely such an extrac-tion and show the results in Fig. 15. These correspondto AFDMC calculations using NN+NNN interactions forthe case where two periods fit inside the box. They arecompared to the free-gas result, which follows from theLindhard function. Overall, we see that the microscopicresults are roughly similar to the free-gas results regard-less of the density. At a more fine-grained level, we ob-serve that the answer to whether or not the microscopicresponse is higher or lower than the free-gas one dependson the density. One could say that this behavior is sim-ilar to what is seen in Fig. 9, but such a comparison ismisleading for two reasons: a) there we were comparingAFDMC results to EDF results, not to free-gas values,and b) our new results in Fig. 15 show the answer forthe response, i.e. at finite one-body potential strength.Thus, these responses cannot be simply extracted fromthe homogeneous-gas answers and constitute microscopic benchmarks. V. SUMMARY & CONCLUSION To summarize, in this article we have investigated theproperties of periodically modulated neutron matter, us-ing a combination of large-scale simulations and qualita-tive insights. We started from the non-interacting prob-lem, examining in detail finite-size effects: since 66 isthe number of particles commonly used for homogeneousneutron matter, we studied the adjustments that need tobe carried out in order to use that particle number forthe inhomogeneous problem. We then reported on ourauxiliary-field diffusion Monte Carlo simulations, under-lining the importance of optimizing the trial wave func-tion by minimizing the VMC energy. This depended ona detailed understanding of the single-particle orbitals.AFDMC allowed us to compute the ground-state en-ergy of neutron matter at various densities, potentialstrengths, and periodicities of the potential. In particu-lar we studied the inhomogeneous problem by increasingthe strength of the potential starting from homogeneity.We then examined several consequences of our ab ini-tio results. We first saw the impact that they have onenergy density functionals. We used the response of neu-tron matter to constrain the isovector term while care-fully disentangling the contributions of bulk and gradi-ent terms. We found a density-dependent isovector termand provided our estimate for its magnitude at each den-sity. Next, we extracted the linear density-density staticresponse function of neutron matter from AFDMC andEDF results at two different densities. This required aset of ab initio results for each of the periodicities thatwe studied. We then compared and contrasted the re-sponse function of neutron matter to that of other sys-tems. More than a proof-of-principle, this work providesdetailed benchmarks that other ab initio calculations cancompare to, or EDF approaches can use as input.The authors are grateful to A. Bulgac, J. Carlson, S.Gandolfi, J. W. Holt, C. J. Pethick, S. Reddy, and A. Riosfor many valuable discussions. They would also like tothank D. Lacroix and A. Pastore for sharing their Skyrmeresults as well as several insights. This work was sup-ported by the Natural Sciences and Engineering ResearchCouncil (NSERC) of Canada and the Canada Founda-tion for Innovation (CFI). Computational resources wereprovided by SHARCNET and NERSC. [1] S. Gandolfi, A. Gezerlis, J. Carlson, Ann. Rev. Nucl.Part. Sci. , 303 (2015).[2] B. Friedman and V. R. Pandharipande, Nucl. Phys. A361 , 502 (1981).[3] A. Akmal, V. R. Pandharipande, and D. G. Ravenhall,Phys. Rev. C , 1804 (1998).[4] A. Schwenk and C. J. Pethick, Phys. Rev. Lett. , 160401 (2005).[5] A. Gezerlis and J. Carlson, Phys. Rev. C , 032801(R)(2008).[6] E. Epelbaum, H. Krebs, D. Lee, and U. -G. Meißner, Eur.Phys. J. A , 199 (2009).[7] N. Kaiser, Eur. Phys. J. A , 148 (2012).[8] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. , 121 (2003).[9] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R.Schaeffer, Nucl. Phys. A , 231 (1998).[10] R. Sellahewa and A. Rios, Phys. Rev. C , 054327(2014).[11] C.J. Yang, M. Grasso, D. Lacroix, Phys. Rev. C ,031301(R) (2016).[12] S. A. Fayans, JETP Lett. , 169 (1998).[13] B. A. Brown, Phys. Rev. Lett. , 420 (2008).[15] F. J. Fattoyev, C. J. Horowitz, J. Piekarewicz, and G.Shen, Phys. Rev. C , 055803 (2010); F. J. Fattoyev,W. G. Newton, J. Xu, B.-A. Li, Phys. Rev. C , 025804(2012).[16] B. A. Brown and A. Schwenk, Phys. Rev. C ,065801 (2016).[18] N. Chamel, S. Goriely, and J. M. Pearson, Nucl. Phys. A , 72 (2008).[19] M. M. Forbes, A. Gezerlis, K. Hebeler, T. Lesinski, A.Schwenk, Phys. Rev. C , 2416 (1996).[22] F. Pederiva, A. Sarsa, K. E. Schmidt, S. Fantoni, Nucl.Phys. A , 255 (2004).[23] S. Gandolfi, J. Carlson, and S. C. Pieper, Phys. Rev.Lett. , 012501 (2011).[24] H. D. Potter, S. Fischer, P. Maris, J.P. Vary, S. Binder,A. Calci, J. Langhammer, R. Roth, Phys. Lett. B ,445 (2014).[25] P. W. Zhao and S. Gandolfi, Phys. Rev. C , 041302(R)(2016).[26] J. Carlson, J. Morales, Jr., V. R. Pandharipande, and D.G. Ravenhall, Phys. Rev. C , 025802 (2003).[27] S. Gandolfi, A. Yu. Illarionov, K. E. Schmidt, F. Ped-eriva, and S. Fantoni, Phys. Rev. C , 054005 (2009).[28] A. Gezerlis and J. Carlson, Phys. Rev. C , 025803(2010).[29] S. Gandolfi, J. Carlson, and S. Reddy, Phys. Rev. C ,032801 (2012).[30] M. Baldo, A. Polls, A. Rios, H.-J. Schulze, and I. Vida˜na,Phys. Rev. C , 064001 (2012).[31] K. Hebeler and A. Schwenk, Phys. Rev. C , 014314(2010).[32] I. Tews, T. Kr¨uger, K. Hebeler, and A. Schwenk, Phys.Rev. Lett. , 032504 (2013).[33] A. Gezerlis, I. Tews, E. Epelbaum, S. Gandolfi, K.Hebeler, A. Nogga, and A. Schwenk, Phys. Rev. Lett. , 032501 (2013).[34] L. Coraggio, J. W. Holt, N. Itaco, R. Machleidt, and F.Sammarruca, Phys. Rev. C , 014322 (2013).[35] G. Hagen, T. Papenbrock, A. Ekstr¨om, K. A. Wendt,G. Baardsen, S. Gandolfi, M. Hjorth-Jensen, and C. J.Horowitz, Phys. Rev. C , 014319 (2014). [36] A. Gezerlis, I. Tews, E. Epelbaum, M. Freunek, S. Gan-dolfi, K. Hebeler, A. Nogga, and A. Schwenk, Phys. Rev.C , 054323 (2014).[37] A. Carbone, A. Rios, A. Polls, Phys. Rev. C , 054322(2014).[38] A. Roggero, A. Mukherjee, F. Pederiva, Phys. Rev. Lett. , 221103 (2014).[39] G. Wlaz lowski, J. W. Holt, S. Moroz, A. Bulgac, and K.J. Roche, Phys. Rev. Lett. , 182503 (2014).[40] I. Tews, S. Gandolfi, A. Gezerlis, and A. Schwenk, Phys.Rev. C , 024305 (2016).[41] N. Iwamoto and C. J. Pethick, Phys. Rev. D , 313(1982).[42] E. Olsson, P. Haensel, and C. J. Pethick, Phys. Rev. C , 025804 (2004).[43] N. Chamel, Phys. Rev. C , 035801 (2012).[44] N. Chamel, D. Page, and S. Reddy, Phys. Rev. C ,035803 (2013).[45] D. Kobyakov and C. J. Pethick, Phys. Rev. C , 055803(2013).[46] D. Davesne, A. Pastore, J. Navarro, Phys. Rev. C ,044302 (2014).[47] A. Pastore, M. Martini, D. Davesne, J. Navarro, S.Goriely, and N. Chamel, Phys. Rev. C , 025804 (2014).[48] D. Davesne, J. W. Holt, A. Pastore, and J. Navarro,Phys. Rev. C , 014323 (2015).[49] A. Pastore, D. Davesne, and J. Navarro, Phys. Rep. ,1 (2015).[50] D. Pines and P. Nozi`eres, The Theory of Quantum Liq-uids, Vol. I , (Benjamin: Reading, 1966).[51] S. Moroni, D. M. Ceperley, G. Senatore, Phys. Rev. Lett. , 1837 (1992).[52] S. Moroni, D. M. Ceperley, G. Senatore, Phys. Rev. Lett. , 689 (1995).[53] M. Buraczynski and A. Gezerlis, Phys. Rev. Lett. ,152501 (2016).[54] B. S. Pudliner, V. R. Pandharipande, J. Carlson, S. C.Pieper, and R. B. Wiringa, Phys. Rev. C , 1720 (1997).[55] K. E. Schmidt and S. Fantoni, Phys. Lett. B , 99(1999).[56] G. Senatore, S. Moroni, D. M. Ceperley, Quantum MonteCarlo Methods in Physics and Chemistry edited by M. P.Nightingale and C. J. Umrigar, p. 183 (Kluwer, 1999).[57] C. Kittel, Solid State Physics , edited by F. Seitz, D.Turnball, H. Ehrenreich , 1, (Academic Press: NewYork, 1968).[58] H. Kwee, S. Zhang, and H. Krakauer, Phys. Rev. Lett. , 126404 (2008).[59] R. B. Wiringa and S. C. Pieper, Phys. Rev. Lett. ,182501 (2002).[60] P. Maris, J. P. Vary, S. Gandolfi, J. Carlson, and S. C.Pieper, Phys. Rev.C , 054318 (2013).[61] J. W. Holt, N. Kaiser, W. Weise, Eur. Phys. J. A , 128(2011).[62] N. Kaiser, Eur. Phys. J. A , 36 (2012).[63] D. Lacroix, private communication (2016).[64] S.-O. Backman, C.-G Kallman, and O. Sjoberg, Phys.Lett.43B