Spin Susceptibility in Neutron Matter from Quantum Monte Carlo Calculations

The spin susceptibility in pure neutron matter is computed from auxiliary field diffusion Monte Carlo calculations over a wide range of densities. Calculations are performed for different spin asymmetries, using twist-averaged boundary conditions to reduce finite-size effects. The employed nuclear interactions include both the phenomenological Argonne AV8$^{\prime}$+UIX potential and local interactions derived from chiral effective field theory up to next-to-next-to-leading order.


Introduction
Spin susceptibility in neutron matter is known to be quite small [1][2][3]. Estimated values are of the order of 10 −3 MeV −1 fm −3 . This suggests that the magnetic fields that are commonly found in neutron stars are very unlikely to produce any significant deviation from the prediction of a model assuming strict zero spin polarization of the constituent matter. However, recent observations have revealed the existence of isolated neutron stars, the so called magnetars, where surface magnetic fields can reach intensities of the order of ∼ 10 14 − 10 15 G [4]. Moreover, in violent phenomena like supernova explosions or neutron star mergers, the fluctuations of the magnetic field magnitude can reach very high peaks [5]. General relativity simulations of such events might therefore be sensitive to the details of the assumed values for the susceptibility [6][7][8]. For this reasons, an accurate benchmark of the existing results is needed.
In this work we present new results for the spin susceptibility in pure neutron matter (PNM) obtained from quantum Monte Carlo (QMC) calculations. We make use of the fact that it is possible to strongly reduce finite-size effects thanks to a modification of the system boundary conditions, overcoming some of the limitations of previous attempts to estimate the susceptibility with QMC methods [1]. In particular, we fix the external magnetic field and use the auxiliary field diffusion Monte Carlo (AFDMC) method [9,10] to study the energy per particle of the system as a function of the spin polarization, keeping the total number of particles fixed. To reduce finite-size effects we use twist-averaged boundary conditions (TABC) [11].
Calculations are carried out for two different nuclear interactions that give realistic mass-radius relations for neutron stars. The first includes a phenomenological two-plus three-neutron potential, namely AV8 +UIX. This interaction has been widely used to study light nuclei and neutron matter properties (see Refs. [12][13][14] and references therein). The second employs local potentials derived from chiral effective field theory (EFT) up to next-to-next-to-leading order (N 2 LO) [15][16][17][18][19]. Among the different available formulations we consider the interaction with coordinate-space cutoff R 0 = 1.0 fm and the E1 parametrization of the three-body contact term V E (see Ref. [19] for details). Such a potential has been successfully used for nuclear structure and nuclear dynamics studies in nuclei [19][20][21][22][23][24][25], and it has been recently employed to derive the equation of state of nuclear matter and the symmetry energy, with good comparison with constraints from terrestrial experiments and multi-messenger astronomy [26].
The manuscript is organized as follows. In Section 2, we provide a short overview of the periodic boundary conditions used in QMC calculations, and the technical details on the implementation of the TABC. In Section 3 we describe the procedures employed to estimate the spin susceptibility starting from the evaluation of the energy as a function of the strength of an external, static magnetic field. Section 4 contains the results of our AFDMC calculations, and Section 5 is devoted to conclusions.

Boundary Conditions in QMC
In QMC calculations, neutron matter is usually modelled as a system of N particles in a box of fixed size L, so that the particle density is ρ = N/L 3 . In order to avoid surface contributions, periodic boundary conditions (PBC) [10] are applied at the borders. This implies that the system is actually made up by an infinite number of boxes, periodically repeated, and all identical. From a dynamical point of view, this means that a particle exiting the box in a given direction re-enters on the opposite side coming from a neighbouring box. Wave functions are chosen coherently to this symmetry, and must be such that: The ansatz for the PNM wave function used in QMC calculations consists of a symmetric Jastrow factor, encoding all the information from the short-range correlations among particles, multiplied by a Slater determinant of plane waves with wave vectors compatible with the periodicity of the system [10]. The definition of the PBC is not univocal. In general, one can allow particles to pick up a phase θ when they wrap around the boundaries. This fact can be expressed by a more general formulation of the boundary conditions (here for a wrap along the x direction): The extra phase θ x is called a "twist" [11]. When the twist angle θ x = 0, Eq. (2) yields back the standard PBC. Allowing for a generic twist angle θ x = 0 defines the TABC. Of course it is possible to introduce separate twists for each system dimension, defining therefore a twist vector.
One of the shortcomings of the standard PBC applied to a homogeneous system of Fermions lays in the fact that the isotropy of the wave function in momentum space is lost. The average over randomly chosen twists helps to restore such a symmetry, thereby reducing shell effects and kinetic energy finite-size errors [11].
Under TABC, physical observables need to be periodic in the twist angle, F(θ i + 2π) = F(θ i ). Due to this symmetry, the values on each component of the twist angle θ i are restricted to the interval: In order to respect the twisted boundary conditions, the wave vectors k n of the plane waves must satisfy: where n is an integer vector.
The energy of each single particle state is given by E n = (h 2 /2m)k 2 n , and the ground state of the system is obtained by filling the lowest energy states. One should notice that the ordering of the single particle energies depends on the value of the twist [11]. As expected from the notion of thermodynamical limit, the contributions of the twist angles to the energy decrease rapidly when increasing the number of particles.
There are two aspects of using TABC that are extremely relevant for the evaluation of the spin susceptibility. First, it significantly reduces the number of particles necessary to describe the infinite system, at the cost of extra computation needed to sample the twist angles (much cheaper). Second, an arbitrary number of spin-up and spin-down particles can be used (in contrast to the PBC, where a number of particles corresponding to closed-shells has to be used), and thus a system with an arbitrary spin polarization can be described.
In Fig. 1 we show results for the energy of PNM as a function of the number of neutrons N using PBC and TABC at saturation density ρ 0 = 0.16 fm −3 (solid symbols). Alongside, we report the energies computed for the free Fermi gas (FG) (dotted-dashed and dashed curves). The results for the free FG are rescaled to match the PNM values asẼ(N) = E(N) − E ∞ + E(66), where E(N) is the energy per particle of the free system, E ∞ is the correct result for the free FG at closed-shell configurations, and E(66) is the energy of PNM with 66 particle. Brown triangles are the PBC results corrected for the kinetic energy of the FG as done in Ref. [27]. The cusps in the energy curves computed using PBC correspond to closed-shell configurations. The use of TABC strongly reduces shell effects: the convergence to the thermodynamic limit is much smoother, and a very reasonable approximation is already obtained for N ≥ 30. To check the robustness of the calculations against a specific choice of the set of twist vectors for selected cases we performed more calculations, sampling new twist vectors for each run. For clarity, the results for a given particle number but different sets of twist vectors are plotted slightly shifted with respect to the correct particle number. Calculations were performed for N = 14, 20, 26, 30, 34, 38, 42, 46, 66. Since the free FG for N = 66 yields an energy particularly close to the thermodynamic limit, QMC calculations of PNM using PBC are usually carried out using that particle number. Presently, AFDMC calculations with realistic nuclear interactions are limited to N ∼ 100.

Computation of the spin susceptibility
The spin susceptibility quantifies the response of a system to the presence of a (static) external magnetic field. The Hamiltonian describing N neutrons in the presence of a magnetic field can be written as: where b = µB, with µ = 6.030774 × 10 −18 MeV/G. H 0 is the non-relativistic Hamiltonian of the interacting system with no external magnetic field, i.e., the sum of a non-relativistic kinetic term and an interaction term that include two-and three-nucleon interactions [13].
The spin susceptibility is defined as: where ρ is the density, and E 0 (b) is the ground-state energy per particle as a function of the magnetic field b. We use two different approaches to estimate the spin susceptibility: (i) Following the approximations of Ref. [1], the derivative in Eq. (6) is estimated from an expansion of the energy as a function of the magnetization, computing the ground-state energy for a few selected values of field strength and spin polarization. (ii) Using directly the definition of the spin polarization ξ where the spin susceptibility can also be written as: The derivative in the above equation is evaluated directly from the value of ξ as a function of the magnetic field strength.

Use of the Pauli expansion
Following Ref. [1], we can use the Pauli expansion of the energy as a function of the polarization ξ. From the expression of the energy per particle E(b) = H /N, and remembering that assuming that b = bk, we obtain: By minimizing E(ξ) with respect to ξ, one can easily rewrite the definition of the spin susceptibility as: The quantity E (0) can be estimated by means of a few ground-state calculations for some fixed spin polarization J z = N ↑ − N ↓ , where N ↑ and N ↓ correspond to Fermion closed-shells in a periodic box (N ↑,↓ = 1, 7, 19, 33, . . .), that can be easily realized within the AFDMC scheme. For any given J z we will have a different value of the energy as a function of the magnetic field E(J z , b), and consequently a different value of the spin polarization ξ(J z ). By using the chain rule, it is possible to express E (0) as: By considering the ground-state energy (∂E 0 /∂J z = 0), the above equation simplifies to: Assuming that: Fig. 3); 3. the spin polarization is linear in J z , the two derivatives in Eq. (14) can be rewritten as: and: where J z0 is the spin asymmetry of the ground state for a given external magnetic field b 0 . These assumptions become exact in the limit of an infinite system with J z and b small. The energies entering Eqs. (15) and (16) can be directly computed using the AFDMC method [10,26]. As previously mentioned, the main difference from the work of Ref. [1] is that with TABC we can access all spin polarizations with a fixed number of particles, while PBC limit the calculations to closed-shell configurations with roughly the same number of particles and only few spin polarizations can be analyzed. Moreover, the spin polarization of the ground state was assumed to be the one of the Fermi gas, while here we compute it directly for the interacting system.

Use of the spin polarization
Another possibility, that also makes explicit use of TABC, is to try to directly determine from AFDMC calculations the ground-state polarization of the neutrons for a given value of the external magnetic field b. By interpolating the values end extracting the b → 0 limit, it is possible to directly use the definition of Eq. (9).

Results
We first verified that the assumptions yielding the approximate expressions of Eqs. (15) and (16) are reliable. In Fig. 2 we report the results for the energy at saturation density for different values of J z /N and b = 0, using N = 38 neutrons. The employed Hamiltonian includes the AV8 +UIX potential. The values are well described by a quadratic fit, at least for values |J z /N| < 0.2. In Fig. 3, the energy as a function of b and for two different fixed values of J z is reported. Here, one expects the energy to be linear as a function of the strength of the magnetic field, and this behaviour is well confirmed. We then performed calculations of the energy per particle for seven different spin asymmetries for each external magnetic field. In order to implement TABC, five different twist vectors have been randomly generated, and independent AFDMC calculations have been performed. The energy for each spin asymmetry is obtained by averaging over the different runs.
The value of J z /N at which the energy is minimized can be assumed as an estimate of the spin polarization of the system in a magnetic field b. Therefore, we performed a quadratic fit on the AFDMC results and determined the energy minimum. Such energy value is then used in Eqs. (15) and (16) to calculate the ground-state spin polarization ξ(b). Uncertainties on the parameters, combined with the statistical errors from the Monte Carlo evaluation, give the total uncertainty on the susceptibility.
In Fig. 4 we report, as an example, the AFDMC results for PNM at saturation density with external field of 20 MeV and neutrons interacting via the AV8 +UIX potential. Similar calculations have been performed at different densities, different external magnetic fields, and also for the local chiral EFT interaction. It is evident that the ground-state spin asymmetry of the interacting system is in general not coincident with the one predicted for a non interacting FG at the same density.  The results of our QMC calculations are reported in Table 1. We obtain significantly smaller spin susceptibilities compared to those from Fantoni et al. [1] for the same Hamiltonian, especially at low densities. This is mainly a consequence of performing our calculation at the correct ground-state spin polarization ξ 0 of the interacting system, rather than at the polarization predicted by the free FG. This is confirmed by the fact that if we compute χ at the the same values of ξ as in Ref. [1], we obtain very similar predictions (see Table 1). Table 1. Spin susceptibility ratio χ/χ F in PNM at different densities (in fm −3 ). Results from Ref. [1] using the two-body potentials AV6 and AV8 supported by the three-body force UIX are reported in the second and third columns. The last two columns show the results of this work for the AV8 +UIX potential. AV8 +UIX ( ) indicates the results obtained at the spin polarization of the ground state ξ 0 predicted using a free FG, to be directly compared with the values reported in Ref. [1].
Ref. [1] This work As shown in Fig. 5, the ground-state spin polarizations computed with and without interaction are significantly different at all densities ρ and at all external magnetic fields b. It is interesting to notice that at densities below saturation, the predicted spin susceptibility is significantly larger, while we obtain the opposite behavior at saturation and higher densities. The results for PNM with the AV8 +UIX potential (solid symbols) are compared to those for the free FG (empty symbols) at three different densities. Analytic results for the FG are shown with solid lines, while linear fits to some mock data (generated from the analytic solution of the free FG at five different external magnetic fields b) are shown with dotted-dashed lines. The results for the interacting system obtained from the quadratic fits of approach (i) are shown with orange-red-brown (from low to high density) dashed lines.
The second approach mentioned in Section 3 directly relies on the definitions of the spin susceptibility and of the spin polarization. As an example, we report in Fig. 5 the results of the ground-state spin polarization ξ 0 as a function of the external magnetic field b for PNM with AV8 +UIX at three different densities. Similar results can also obtained for local chiral EFT potentials, and at different densities.
In Fig. 6 we report the estimates of the spin susceptibility in PNM with neutrons interacting with both the phenomenological AV8 +UIX potential and the local chiral EFT interaction. A comparison between the two approaches discussed in the previous section for the AV8 +UIX case is also reported in the same figure. The direct estimation of the susceptibility gives a smoother dependence on the density. This fact might be related to the absence of residual systematic errors that are still likely to plague the indirect estimate described above. For as concerns the results from the EFT Hamiltonian, they do not significantly differ from those obtained from a a phenomenological interaction over all the density range considered. In order to estimate the uncertainty due to the assumption of a linear dependence of the energy on the magnetic field b, we can once again resort to the FG analysis. In fact, in the case of the non-interacting system we can compute the analytic solution and analyze the discrepancies with the results obtained under the assumption of a linear dependence. As shown in Fig. 7, the linear approximation yields to slightly underestimated values for the susceptibility, in particular at larger densities. In principle, there is no reason to assume a different behaviour for an interacting system, although the correct result is unknown.
In Fig. 8 we compare our results for the spin susceptibility to those available from Brueckner-Hartree Fock (BHF) calculations. In black, we report the calculations by Bombaci et al. [3] for the two-body AV18 potential. The results of Vidaña et al. [2], obtained using the two-body NSC97e interaction in a parametrization for high-spin polarizations, are shown with the orange band. The results look overall qualitatively consistent, even though a fair comparison is difficult due to the different employed nuclear potentials, many-body methods, and strategy used to estimate the spin susceptibility. For phenomenological potentials, the effect of the three-body force appears to be non negligible, as shown by the ≈ 10 − 15% difference between BHF AV18 and AFDMC AV8 +UIX results. Realistic interaction models including both two-and three-body forces (AV8 +UIX and local chiral EFT) appear to provide consistent susceptibilities, regardless the scheme of the interaction. Finally, we analyze an alternative procedure of computing the spin susceptibility that only relies on the equation of state (EOS) of spin unpolarized PNM and fully spin polarized PNM. We first define an energy density functional ε of the form: where ε 0 (ρ) and ε 1 (ρ) are two parametrization with the same functional form for the spin unpolarized and fully spin polarized system, respectively. The quadratic dependence of the energy as a function of the spin polarization ξ is a common assumption deriving from the expansion of the energy of the partially polarized FG. Inserting the energy function ε in the definition of Eq. (6), we obtain that the spin susceptibility can be calculated as: where we can define the spin symmetry energy as E σ sym (ρ) = ε 1 (ρ) − ε 0 (ρ). This holds for any expression of ε i (ρ). Using AFDMC calculations for spin unpolarized PNM (from Ref. [27] for the AV8 +UIX potential, and from Ref. [18] for the local chiral EFT interaction) and fully spin polarized PNM (from Ref. [28]), we directly compute the spin symmetry energy, as shown in Fig. 9.
Using Eq. (18), we can directly compute the spin susceptibility. The results are shown in Fig. 10 and are compared to the values obtained using the approach (ii). The spin susceptibility for the phenomenological potential provides a better agreement between the two approaches, while the comparison is worse for the local chiral EFT interaction.
This approach, although really simple, provides a reasonable estimate for the spin susceptibility in PNM. In contrast, the free FG solution predicts a spin susceptibility larger by a factor ≈ 3 than the one obtained for the interacting case.

Conclusions
We performed substantially improved AFDMC calculation of the spin susceptibility in PNM. The main result is that the predicted ground-state polarization of the interacting system is lower than the  Figure 8. Spin susceptibility ratio χ/χ F in PNM as a function of the density. Green squares (red circles) are results for the AV8 +UIX (local chiral EFT) potential from approach (ii). The black line shows the BHF results of Ref. [3] for the two-body interaction AV18. The orange band reports the BHF results of Ref. [2] for the two-body NSC97e interaction in a parametrization for high-spin polarizations.
one predicted by the free FG. Although previous calculations did not take such an effect into account, this does not lead to significant changes in the magnitude of the spin susceptibility. Implementing TABC allowed us to check the reliability of the assumptions made in previous works [1]. Moreover, the spin susceptibility is now calculated more consistently: on the one hand we can perform calculations with an arbitrary number of particles, therefore reducing significantly finite-size effects; on the other hand we can perform calculations with arbitrary spin polarization, keeping the total number of particles fixed. The same technique can also be directly extended to isospin-asymmetric matter. We performed and compareed calculations using both the phenomenological potential AV8 +UIX and a local chiral EFT interaction up to N 2 LO. We presented a method, namely approach (ii), to directly compute the spin susceptibility from the energy of a system subject to an external magnetic field, that relies on the possibility of performing AFDMC calculations for arbitrary number of particles in a periodic box, minimizing the finite size effects by means of TABC. This approach has an analogous counterpart for the free FG, which provides a reasonable way to estimate uncertainties.    Figure 10. Same as Fig. 6 but results are shown for PNM from approach (ii) and from E σ sym .