Vortices in low-density neutron matter and cold Fermi gases

TL;DR

This study uses diffusion Monte Carlo simulations to analyze vortices in low-density neutron matter and cold Fermi gases, revealing similarities and differences in vortex properties and density profiles, with implications for neutron star physics.

Contribution

It provides the first detailed comparison of vortex properties in low-density neutron matter and cold atomic gases using quantum Monte Carlo methods.

Findings

Agreement between cold gases and neutron matter at very low densities.

Density depletion at vortex cores is about 25% in neutron matter.

Vortex excitation energies vary with interaction strength.

Abstract

Cold gas experiments can be tuned to achieve strongly-interacting regimes such as that of low-density neutron matter found in neutron-stars' crusts. We report $T$=0 diffusion Monte Carlo results (i) for the ground state of both spin-1/2 fermions with short-range interactions and low-density neutron matter in a cylindrical container, and (ii) properties of these systems with a vortex line excitation. We calculate the equation of state for cold atoms and low-density neutron matter in the bulk systems, and we contrast it to our results in the cylindrical container. We compute the vortex line excitation energy for different interaction strengths, and we find agreement between cold gases and neutron matter for very low densities. We also calculate density profiles, which allow us to determine the density depletion at the vortex core, which depends strongly on the short-ranged interaction in cold atomic gases, but it is of $\approx$ 25% for neutron matter in the density regimes studied in this work. Our results can be used to constrain neutron matter properties by using measurements from cold Fermi gases experiments.

Figures

Tables

Table 1. Table 1: Bulk energies per particle and the parameters E 0 cyl superscript subscript 𝐸 0 cyl E_{0}^{\rm cyl} and λ S subscript 𝜆 𝑆 \lambda_{S} fitted to the functional form of Eq. ( 1 ). The bulk energies and E 0 cyl superscript subscript 𝐸 0 cyl E_{0}^{\rm cyl} are reported in units of the free Fermi gas energy, E F ​ G subscript 𝐸 𝐹 𝐺 E_{FG} (see Eq. ( 2 )), while λ S subscript 𝜆 𝑆 \lambda_{S} is reported in units of E F ​ G ​ k F − 2 subscript 𝐸 𝐹 𝐺 superscript subscript 𝑘 𝐹 2 E_{FG}k_{F}^{-2} .

	Cold gases			Neutron matter
				$s$-wave AV18			modified Poschl-Teller
$-k_{F}a$	bulk	$E_0^{\mathrm{cyl}}$	${\lambda }_S$	bulk	$E_0^{\mathrm{cyl}}$	${\lambda }_S$	bulk	$E_0^{\mathrm{cyl}}$	${\lambda }_S$
0.5	0.8636(1)	0.90(3)	95(15)
1.0	0.7864(2)	0.79(1)	99(7)	0.814(5)	0.76(2)	-78(9)	0.821(4)	0.98(9)	204(54)
2.0	0.6806(2)	0.72(4)	70(19)	0.748(2)	0.75(4)	-242(10)	0.749(3)	0.70(7)	117(30)
3.3	0.5979(2)	0.66(1)	58(4)	0.667(2)	0.67(2)	-365(70)	0.681(2)	0.68(5)	85(24)
5.0	0.5407(2)	0.60(1)	64(5)	0.598(2)	0.60(5)	-445(80)	0.608(1)	0.61(3)	85(14)

Equations

$$E^{\mathrm{cyl}}(\mathcal{R},\mathcal{L})=E_0^{\mathrm{cyl}}+\frac{{\lambda }_S}{2\pi \mathcal{R}\mathcal{L}},$$

$$E_{FG}=\frac{3}{10}\frac{{\hslash }^2}{m}{\left(3{\pi }^2\frac{N}{V}\right)}^{2/3}.$$

$$k\cot \delta (k)=-\frac{1}{a}+\frac{r_e{k}^2}{2}+\mathcal{O}(k^4),$$

$$V_{\mathrm{mPT}}(r)=-v_0\frac{{\hslash }^2}{m_r}\frac{{\mu }^2}{{\cosh }^2(\mu r)},$$

$$a\mu =\frac{\pi }{2}\cot \left(\frac{\pi \lambda }{2}\right)+\gamma +\Psi (\lambda ),$$

$${\psi }_{\mathrm{BCS}}(\text{R},S)=\mathcal{A}[\varphi ({\text{r}}_1,s_1,{\text{r}}_{1^{\prime }},s_{1^{\prime }})$$

$$\varphi ({\text{r}}_2,s_2,{\text{r}}_{2^{\prime }},s_{2^{\prime }})\dots \varphi ({\text{r}}_{N/2},s_{N/2},{\text{r}}_{N/2^{\prime }},s_{N/2^{\prime }})],$$

$$\varphi (\text{r},s,{\text{r}}^{\prime },s^{\prime })=\tilde{\varphi }(\text{r},{\text{r}}^{\prime })\left[\frac{⟨s{s}^{\prime }|\uparrow \downarrow ⟩-⟨s{s}^{\prime }|\downarrow \uparrow ⟩}{\sqrt{2}}\right],$$

$$\left[-\frac{{\hslash }^2}{m}{\nabla }^2+V(r)\right]f(r<d)=\lambda f(r<d),$$

$${\psi }_{\mathrm{T}}(\text{R},S)=\prod _{i,j^{\prime }}f(r_{i{j}^{\prime }}){\psi }_{\mathrm{BCS}}(\text{R},S).$$

$${\tilde{\varphi }}_{\mathrm{bulk}}(\text{r},{\text{r}}^{\prime })=\sum _{n=1}^{n_c}{\alpha }_n{e}^{i{\text{k}}_{\text{n}}\cdot (\text{r}-{\text{r}}^{\prime })}+\tilde{\beta }(|\text{r}-{\text{r}}^{\prime }|),$$

$$\tilde{\beta }(r)=\{\begin{array}{ll}\beta (r)+\beta (L-r)-2\beta (L/2)& \text{for }r⩽L/2\\ 0& \text{for }r>L/2\end{array}$$

$$\beta (r)=[1+cbr][1-e^{-dbr}]\frac{e^{-br}}{dbr},$$

$${\Phi }_{n\nu p}(\rho ,\phi ,z)={\mathcal{N}}_{\nu p}{J}_{\nu }(k_{\nu p}\rho )\exp \left[i(k_{z}z+\nu \phi )\right],$$

$${\tilde{\varphi }}_{\mathrm{cyl}}(\text{r},{\text{r}}^{\prime })=\sum _{q=1}^{q_c}{\tilde{\alpha }}_q{\mathcal{N}}_{\nu p}^2{J}_{\nu }\left(\frac{j_{\nu p}}{\mathcal{R}}\rho \right)J_{\nu }\left(\frac{j_{\nu p}}{\mathcal{R}}{\rho }^{\prime }\right)$$

$$\times e^{i\nu (\phi -{\phi }^{\prime })}{e}^{i{k}_z(z-z^{\prime })}+\bar{\beta }(\text{r},{\text{r}}^{\prime }),$$

$$\bar{\beta }(\text{r},{\text{r}}^{\prime })=\{\begin{array}{ll}{\mathcal{N}}_{01}^2{J}_0\left(\frac{j_{01}\rho }{\mathcal{R}}\right)J_0\left(\frac{j_{01}{\rho }^{\prime }}{\mathcal{R}}\right)\times & \\ \left[\beta (r)+\beta (2\mathcal{R}-r)-2\beta (\mathcal{R})\right]& \text{for }r⩽\mathcal{R}\\ 0& \text{for }r>\mathcal{R}\end{array}$$

$${\tilde{\varphi }}_{\mathrm{vortex}}(\text{r},{\text{r}}^{\prime })=\sum _{q=1}^{q_c}{\bar{\alpha }}_q{\mathcal{N}}_{\nu p}{\mathcal{N}}_{\nu -1;p}\times$$

$$\{J_{\nu }\left(\frac{j_{\nu p}}{\mathcal{R}}\rho \right)J_{\nu -1}\left(\frac{j_{\nu -1;p}}{\mathcal{R}}{\rho }^{\prime }\right)e^{i(\nu \phi -(\nu -1){\phi }^{\prime })}{e}^{i{k}_z(z-z^{\prime })}$$

$$+J_{\nu }\left(\frac{j_{\nu p}}{\mathcal{R}}{\rho }^{\prime }\right)J_{\nu -1}\left(\frac{j_{\nu -1;p}}{\mathcal{R}}\rho \right)e^{i(\nu {\phi }^{\prime }-(\nu -1)\phi )}{e}^{i{k}_z(z^{\prime }-z)}\},$$

$$H=-\frac{{\hslash }^2}{2m}\left[\sum _{i=1}^{N_{\uparrow }}{\nabla }_i^2+\sum _{i=j^{\prime }}^{N_{\downarrow }}{\nabla }_{j^{\prime }}^2\right]+\sum _{i,j^{\prime }}V(r_{i{j}^{\prime }}),$$

$$\psi (\tau )=e^{-(H-E_T)\tau }{\psi }_T,$$

$$\underset{\tau \to \infty }{\lim }\psi (\tau )={\Phi }_0.$$

$$\psi (\text{R},\tau )=\int d{\text{R}}^{\prime }G(\text{R},{\text{R}}^{\prime },\tau ){\psi }_T({\text{R}}^{\prime }),$$

$$⟨O(\text{R}){⟩}_m=\frac{⟨{\Psi }_T(\text{R})|O(\text{R})|{\Phi }_0(\text{R})⟩}{⟨{\Psi }_T(\text{R})|{\Phi }_0(\text{R})⟩},$$

$$⟨O(\text{R}){⟩}_e=\frac{⟨{\Phi }_0(\text{R})|O(\text{R})|{\Phi }_0(\text{R})⟩}{⟨{\Phi }_0(\text{R})|{\Phi }_0(\text{R})⟩}.$$

$${\int }_V{d}^{3}r\mathcal{D}(\rho )=1,$$

$$\eta (R)=N{\int }_0^{\mathcal{L}}dz{\int }_0^{2\pi }d\phi {\int }_0^{R}d\rho \rho \mathcal{D}(\rho ),$$

$$\underset{q\to 0}{\lim }\frac{{\delta }_0}{2q}=\frac{1}{\lambda }-\frac{\pi }{2}\cot \left(\frac{\pi \lambda }{2}\right)+\sum _{n=1}^{\infty }\left(\frac{1}{\lambda +n}-\frac{1}{n}\right)=$$

$$\frac{1}{\lambda }-\frac{\pi }{2}\cot \left(\frac{\pi \lambda }{2}\right)+\sum _{n=1}^{\infty }\frac{-\lambda }{n(\lambda +n)},$$

$$\Psi (1+z)$$

$$\underset{q\to 0}{\lim }\frac{{\delta }_0}{2q}=-\frac{\pi }{2}\cot \left(\frac{\pi \lambda }{2}\right)-\gamma -\Psi (\lambda ).$$

Full text

I Introduction

Strongly interacting fermionic systems appear in many contexts, for example: superconductors, cold atomic Fermi gases, low-density neutron matter, and QCD at high baryon densities. Shedding light on properties of one of these systems may contribute to our comprehension of strongly interacting Fermi systems as a whole.

Cold atom systems provide an example where the interplay between experiments and theory led to rapid advances in the field. In these dilute systems, short-range interactions are characterized by a single parameter $k_{F}a$, the product of the Fermi wave number $k_{F}$ and the $s$-wave scattering length $a$. This interaction strength can be tuned using an external magnetic field near a Feshbach resonance, and the attractive interactions can span a continuum between the Bardeen-Cooper-Schrieffer (BCS) limit of superfluidity and the Bose-Einstein condensation (BEC) of dimers, passing through the unitary limit of infinite scattering length. Experiments with cold atoms can provide direct tests of quantities such as the equation of state and pairing gap, which are currently inaccessible to their neutron matter counterparts. For a review on the subject the reader is referred to Ref. Giorgini et al. (2008) and references therein.

On the other hand, if we compare cold gases to neutron matter, we find that the neutron-neutron interaction can be more complicated: short-range repulsion, two-pion exchange at a intermediate range, and one-pion exchange at large distances. However, this situation changes in the low-density regime, which is the case in the exterior of neutron-rich nuclei and neutron-star crusts. In these systems, the scattering length and effective range of the interaction are the most essential quantities for describing the physical properties, and properties of neutron matter and cold atoms are similar Carlson et al. (2012); Gandolfi et al. (2015).

A neutron matter model with a zero-range interaction Baker (2001) was presented as a “many-body challenge” proposed by Bertsch, 111 The challenge proposed to the participants of the Tenth International Conference on Recent Progress in Many-Body Theories can be stated as: what are the ground state properties of the many-body system composed of spin-1/2 fermions interacting via a zero-range, infinite scattering length contact interaction? , much before cold atom experiments could shed light on the properties of these systems. In dilute cold gases, the effective range $r_{e}$ between atoms is much smaller than the interatomic spacing $r_{0}$, and can be taken to be zero. The diluteness can guarantee that the scattering length $a$ is much larger than $r_{0}$. Comparison with other systems is meaningful if they also obey $|a|\gg r_{0}\gg r_{e}$. The scattering length of neutron matter, $a^{nn}\sim-18.5$ fm, is substantially larger than the interparticle distance and the effective range, $r_{e}^{nn}\sim 2.7$ fm, such that $|r_{e}^{nn}/a^{nn}|\approx 0.15$. However, only at very low-densities is the effective range much smaller than the interparticle distance. If we neglect the effects of a finite effective range in the neutron-neutron interaction, cold atoms and neutron matter are universal in the sense that properties depend only on the product $k_{F}a$.

Quantum Monte Carlo (QMC) methods have been successful at comparing the equation of state and pairing gap of cold atom systems and low-density neutron matter Gezerlis and Carlson (2008, 2010). In the present work we used a similar model to compare properties of vortices in low-density neutron matter and cold Fermi gases. One signature of superfluidity is the formation of quantized vortices, where the quantization of the flow is given in units of $h/(2m)$, $m$ being the mass of the fermion. The microscopic structure of a vortex line in neutron matter has been studied using Bogoliubov-de Gennes equations De Blasio and Elgarøy (1999); Elgarøy, Ø. and De Blasio, F. V. (2001), and nuclear energy density functional approaches Yu and Bulgac (2003). For cold atom gases there is an abundance of studies, for example $T=0$ results using Bogoliubov-de Gennes equations at unitarity Bulgac and Yu (2003) and throughout the BEC-BCS crossover Sensarma et al. (2006), and finite temperature calculations Simonucci et al. (2013).

Here we report results for a single vortex line in a cylindrical geometry for both low-density neutron matter and cold Fermi gases using QMC methods. We investigated the consequences of the finite effective range of the neutron-neutron interaction, in contrast to $r_{e}\approx 0$ for cold gases. We also studied effects that go beyond low-energy scattering by using two potential models for the neutron-neutron interaction, one based on phenomenology, and another that was tuned to reproduce the desired low-energy phase shifts. We calculated the equation of state for cold atoms and low-density neutron matter in the bulk systems. We show that it is possible to separate the energy contributions of systems in a cylindrical container with hard walls into bulk and surface terms. The excitation energy necessary to produce a vortex line was computed by using the energy difference between a system of pairs with angular momentum $\hbar$ and the ground state. We show that for very low densities there is an agreement between the excitation energies for vortex-line formation between cold gases and neutron matter. However, as the density increases (or as the interaction strength increases in absolute value) they differ. We also calculated density profiles, which allows us to determine the density depletion at the vortex core. We found that the depletion varies from 28% up to 47% for cold gases, whereas for neutron matter the depletion is approximately 25%, for the density range we studied in this work. Our results are compared to previous mean-field calculations.

This paper is structured as follows. In Sec. II we introduce our methodology. We discuss aspects of the cylindrical container in Sec. II.1, and low-energy two-body scattering in Sec. II.2. We present the wave functions we built in Sec. II.3, which describe properties of the bulk systems, and systems in a cylindrical container (both the ground state and systems with a vortex line). In Sec. II.4 we give a brief description of the QMC methods we employed. Section III presents our results, namely the ground state and vortex excitation energies in Sec. III.1, and density profiles in Sec. III.2. An outlook is provided in Sec. IV. Finally, in Appendix A we show how to obtain an exact relationship between scattering length and the parameters of the modified Poschl-Teller potential.

II Methods

II.1 Cylindrical container

The choice of which trapping potential (or geometry) to use in this problem is not unambiguous, as there is a trade-off for each possible candidate. A choice that minimizes surface effects is to have an array of counter rotating vortices with periodic boundary conditions. One drawback is that this state has zero total angular momentum, thus it can decay to the ground state of the system. Also, from the computational perspective, this choice is not feasible for fermionic systems. For example, 4He calculations of Ref. Sadd et al. (1997) used 300 particles and four counter rotating vortices in the simulation cell. In order to use the same number of fermion pairs we would require a system of 600 fermions. Another possible choice would involve harmonic traps, which are readily available in experimental setups, however, the density profiles of cold gases in harmonic traps can also differ substantially from what is expected in the thermodynamic limit Chang and Bertsch (2007).

Instead, we opted for using a cylindrical container of radius ${\cal R}$ and height ${\cal L}$, with hard walls, periodic in the axial direction. This choice is consistent with previous bosonic Vitiello et al. (1996) and fermionic Madeira et al. (2016) calculations. Also, this is the generalization of the two-dimensional (2D) disk geometry to 3D Ortiz and Ceperley (1995); Giorgini et al. (1996); Madeira et al. (2017), where we made the axial direction periodic. Throughout this work we use $(\rho,\varphi,z)$ to denote the usual cylindrical coordinates.

In the thermodynamic limit, ${\cal R,L}\to\infty$, the energy per particle is independent of the cylinder radius and height, and it should go to the bulk value. The relationship between thermodynamic properties of a confined fluid and the shape of the container is often expressed as a function of the various curvatures of the container König et al. (2004). For these reasons, we chose the following functional form for the energy per particle in the cylindrical geometry:

[TABLE]

where $E_{0}^{\rm cyl}$ represents a bulk contribution to the energy, and the second term on the right-hand side is a “surface” contribution. Corrections to this functional form would come in powers of ${\cal R}^{-1}$ and/or ${\cal L}^{-1}$, however, we found those do not improve the description of the results.

One of the complications of introducing hard walls is the presence of the so-called Friedel oscillations. Fermionic systems bound by hard walls display density profiles characterized by Friedel oscillations. Although they are present in three dimensions, they are more pronounced in low-dimension systems such as 1D McKenney et al. (2016), and 2D Madeira et al. (2017). We would like our system to exhibit some desirable features with respect to the energy and density distribution ${\cal D}(\rho)$. Regarding the density distribution as a function of the radial coordinate (see Sec. III.2 for the normalization and profiles in the interacting cases), besides a vanishing density at the walls, we want the profile to be flat close to the axis of the cylinder. This would be the behavior in the thermodynamic limit, but this is not always true for finite-size systems. If we fix the number density at $n=k_{F}^{3}/(3\pi^{2})$, the free Fermi gas density, we have freedom to choose either the cylinder radius ${\cal R}$ or the height ${\cal L}$. In making this choice we adopted the following procedure. We calculated analytically the energy and density profile for the free gas, and we looked for systems that obeyed the criteria established above, that is: (i) the energy of the system for different particle numbers $N$ is well described by Eq. (1); (ii) the slope of the density profile in the vicinity of the origin ($\rho\leqslant\rho_{0}$) is less than a prescribed tolerance, $|\partial{\cal D}(\rho)/\partial\rho|_{\rho\leqslant\rho_{0}}\leqslant\epsilon$; (iii) density oscillations are minimized.

Throughout this work we report energies per particle in units of the free Fermi gas energy per particle,

[TABLE]

We found that for $N=\{78,80,82,84,86\}$, with the radius ${\cal R}=9.3k_{F}^{-1}$ for $N$=78 and $9.4k_{F}^{-1}$ for the other systems, are well described by Eq. (1) with $E_{0}^{\rm cyl}=1.02(2)E_{FG}$ and $\lambda_{S}=108(8)E_{FG}k_{F}^{-2}$, and the maximum value of $\epsilon$ is of $\sim 5\times 10^{-5}k_{F}^{4}$ for $\rho_{0}\leqslant 2k_{F}^{-1}$. To illustrate how criteria (ii) and (iii) are not easily met, we plot in Fig. 1 the density profile for the free Fermi gas with $N=78$ and several different values of ${\cal R}$. Our ansatz takes into consideration only the free gas case, which corresponds to the $-k_{F}a\to 0$ limit. However, we show in Secs. III.1.1 and III.2 that our choices produced the desired results in the $0.5\leqslant-k_{F}a\leqslant 5.0$ range.

II.2 Scattering

Two-body scattering for a finite range potential $V(r)$ is described by the Schrödinger equation. We separate the solutions into radial and angular parts, with the latter being a constant for $s$-wave scattering. The scattering length $a$ and the effective range $r_{e}$ can be determined from the zero-energy solution of the radial equation and its asymptotic form. The low-energy behavior of the phase shift $\delta(k)$ can be related to $a$ and $r_{e}$ Bethe (1949),

[TABLE]

hence different potentials that reproduce the same scattering length and effective range yield the same low-energy phase-shift behavior. When simulating cold gases, we chose the modified Poschl-Teller (mPT) potential to describe interactions between antiparallel spins,

[TABLE]

where $v_{0}$ and $\mu$ are parameters that can be tuned to reproduce the desired $a$ and $r_{e}$. We restricted the parameters so that no bound state is supported. The quantities $a$, $\mu$, and $v_{0}$ are related through (see Appendix A)

[TABLE]

where $\gamma=0.577\dots$ is the Euler-Mascheroni constant, $\Psi$ is the digamma function, and $\lambda$ is such that $v_{0}=\lambda(\lambda-1)/2$. In the equation above, the requirement on the number of bound states, and a fixed $r_{e}$, completely determine the parameters of the potential for a given scattering length.

For the neutron matter simulations, we employed two different potential interactions. Our goal with this approach is to see if there are any relevant effects beyond the low-energy regime described by Eq. (3). The first interaction we considered is a modified Poschl-Teller potential, Eq. (4), tuned so that the scattering length is $a^{nn}=-18.5$ fm and the effective range is $r_{e}^{nn}=2.7$ fm. The other one is based on the AV18 nucleon-nucleon pairwise interaction Wiringa et al. (1995), which has been extensively used in QMC simulations of nucleon systems Carlson et al. (2015). We chose the neutron-neutron interaction between particles with antiparallel spins to be the $s$-wave part of AV18. We fixed the spin-isospin degrees of freedom such that we have a unpolarized gas of neutrons, hence the potential interaction becomes spherically symmetrical. The most important feature of the interaction is that the scattering length $a^{nn}=-18.5$ fm and effective range $r_{e}^{nn}=2.7$ fm are correctly described by the potential. In Fig. 2 we compare the potential interactions we use for cold gases and neutron matter for $-k_{F}a=1$.

II.3 Wave functions

The BCS wave function, which includes pairing explicitly, projected to a fixed number of particles $N$ (half with spin-up and half with spin-down), can be written as an antisymmetrized product Bouchaud, J.P. et al. (1988). Since neither the Hamiltonian or any operators in the quantities we calculate flip the spins, we adopt hereafter the convention of primed indices to denote spin-down particles and unprimed ones to refer to spin-up particles. Thus, the BCS wave function reduces to

[TABLE]

where R is a vector containing the particle positions $\textbf{r}_{i}$, $S$ stands for the spins $s_{i}$, and the antisymmetrization is over spin-up and spin-down particles only Carlson et al. (2003). This wave function can be calculated efficiently as a determinant. The $\phi$ are pairing functions, which have the form

[TABLE]

where we have explicitly included the spin part to impose singlet pairing. The assumed expressions for $\tilde{\phi}$ depend on the system being studied, see Secs. II.3.1, II.3.2, and II.3.3.

The BCS wave function accounts for the long-range behavior. Short-range correlations are included in the form of a two-body Jastrow factor $f(r_{ij^{\prime}})$, $r_{ij^{\prime}}=|\textbf{r}_{i}-\textbf{r}_{j^{\prime}}|$, which accounts for correlations between antiparallel spins. It is obtained from solutions of the two-body Schrödinger-like equation,

[TABLE]

where $V(r)$ is specified for cold gases and neutron matter in Sec. II.2, and the boundary conditions are $f(r>d)=1$ and $f^{\prime}(r=d)=0$, $d$ being a variational parameter, and $\lambda$ is adjusted so that $f(r)$ is nodeless. The total trial wave function is written as

[TABLE]

II.3.1 Bulk system

We employed the same pairing function for the bulk case as Ref. Carlson et al. (2003),

[TABLE]

where $\alpha_{n}$ are variational parameters, and contributions from momentum states up to a level $n_{c}$ are included. The $\tilde{\beta}$ function describes contributions with $n>n_{c}$,

[TABLE]

with

[TABLE]

where $r=|\textbf{r}-\textbf{r}^{\prime}|$ and $b$, $c$, and $d$ are variational parameters. We considered $b$ = 0.5 $k_{F}$, $d=5$, and $c$ is adjusted so that $\partial\tilde{\beta}/\partial r=0$ at $r=0$. This functional form of $\beta(r)$ describes the short-distance (high-momentum) correlation of particles with antiparallel spins.

II.3.2 Cylinder

The free-particle solution of the Schrödinger equation in a cylinder or radius ${\cal R}$, height ${\cal L}$, finite at $\rho=0$, and with periodic conditions along the $z$-axis is

[TABLE]

where ${\cal N}_{\nu p}$ is a normalization constant, $J_{\nu}$ are Bessel functions, $k_{\nu p}=j_{\nu p}/{\cal R}$, $j_{\nu p}$ is the $p$-th zero of $J_{\nu}$, and $k_{z}=2\pi n/{\cal L}$. The eigenvalues are $E_{n\nu p}=\hbar^{2}(k_{\nu p}^{2}+k_{z}^{2})/(2m)$. The quantum numbers $n$ and $\nu$ can take the values 0,$\pm 1,\pm 2,\cdots$, and $p=1,2,\cdots$.

The pairing function for the cylinder geometry is constructed using the single-particle orbitals of Eq. (13) coupled with their time-reversed counterparts. This ansatz has been used before in the unitary Fermi gas Madeira et al. (2016). We assume the pairing function to be

[TABLE]

where the $\tilde{\alpha}_{q}$ are variational parameters, and $q$ is a label for the cylinder momentum shells, such that different states with the same energy have the same variational parameter. The $\bar{\beta}$ function is a modification of $\tilde{\beta}$ such that the hard wall boundary condition is met,

[TABLE]

and $\beta$ is given by Eq. (12).

II.3.3 Vortex

The vortex line excitation is accomplished by considering pairing orbitals which are eigenstates of $L_{z}$ with eigenvalues $\pm\hbar$. This is achieved by coupling single-particle states with angular quantum numbers $\nu$ differing by one. Explicitly, we are considering $(n,\nu,p)$ paired with $(-n,-\nu+1,p)$, such that the pairing orbitals take the form

[TABLE]

where $q$ is a label for the vortex shells, and $\bar{\alpha}_{q}$ are variational parameters. Equation (II.3.3) is symmetric under interchange of the prime and unprimed coordinates, as required for singlet pairing.

II.4 Quantum Monte Carlo

The Hamiltonian of the two-component Fermi gas, or spin-up/spin-down neutron matter, is given by

[TABLE]

with $N=N_{\uparrow}+N_{\downarrow}$. The diffusion Monte Carlo (DMC) method projects out the lowest-energy state of $H$ present in a initial state $\psi_{T}$ (obtained from variational Monte Carlo simulations). The propagation, in imaginary time $\tau$, can be written as

[TABLE]

where $E_{T}$ is an energy offset. In the $\tau\to\infty$ limit, only the lowest energy component $\Phi_{0}$ survives

[TABLE]

The imaginary time evolution can be written in the integral form

[TABLE]

where $G(\textbf{R},\textbf{R}^{\prime},\tau)$ is the Green’s function associated with $H$. We solve an importance sampled version of Eq. (20) iteratively, using the Trotter-Suzuki approximation to evaluate $G(\textbf{R},\textbf{R}^{\prime},\tau)$, which requires the time steps $\delta\tau=\tau/N$ to be small. We circumvent the fermion-sign problem by using the fixed-node approximation, which restricts transitions across a nodal surface defined by $\psi_{T}$, making our estimates of energy expectation values upper bounds. For a detailed explanation of the algorithm, the reader is referred to Ref. Foulkes et al. (2001) and references therein.

The direct calculation of the expectation value of an operator $O(\textbf{R})$ from $\Phi_{0}(\textbf{R})$ corresponds to the mixed estimator

[TABLE]

which is exact only when $O$ commutes with the Hamiltonian $H$. There are several methods to compute expectation values of quantities, such as the density, that do not commute with $H$. One of them is the extrapolation method where the results of diffusion and variational simulations are combined. However, the accuracy of the extrapolation method relies completely on the trial wave function. Moreover, even in the case of accurate trial wave functions, the bias of the extrapolated estimator is difficult to calculate. For these reasons we used the forward walking method, which is discussed in detail in Ref. Casulleras and Boronat (1995), to evaluate the density profiles. This method relies on the calculation of the asymptotic offspring of walkers coming from the branching term to compute the exact estimator,

[TABLE]

The variational parameters in Eqs. (10), (II.3.2), and (II.3.3) were determined using the linear method Toulouse and Umrigar (2007). In this method, parameter variations are found by diagonalizing a non-symmetric estimator of the Hamiltonian matrix in the basis of the wave function and its derivatives with respect to the parameters. We also adopted the heuristic procedure of Ref. Contessi et al. (2017), which suppresses instabilities that arise from the non-linear dependence of the wave function on the variational parameters.

III Results

Comparison between cold atom systems and low-density neutron matter is achieved by expressing energies (per particle) in units of the free Fermi gas energy, see Eq. (2), and distances in units of $k_{F}^{-1}$. For the cold gases systems we keep the effective range fixed at $k_{F}r_{e}=0.05$, which is much smaller than the interparticle spacing and the scattering lengths involved in the simulations. The number density is kept constant at $n=k_{F}^{3}/(3\pi^{2})$. For bulk systems this corresponds to $n=N/L^{3}$, and for cylindrical containers $n=N/(\pi{\cal R}^{2}{\cal L})$. The interaction strengths we considered for cold gas systems are $-k_{F}a=\{0.5,1.0,2.0,3.3,5.0\}$, while we do not include the $-k_{F}a=0.5$ case for neutron matter because it is extremely dilute, and Friedel oscillations prevent any meaningful analysis of the density profiles.

III.1 Energy

III.1.1 Ground state energy

The ground state energy per particle of the bulk systems for several values of $k_{F}a$ was calculated using the pairing function of Eq. (10), and the results are shown in Table 1 and also in Fig. 3. The energy per particle of the cold atoms systems is lower than the neutron matter systems, for the same value of $k_{F}a$, in accordance with previous simulations. In fact, our results for cold atoms are lower than those reported in Ref. Gezerlis and Carlson (2008) because we chose a smaller effective range, $k_{F}r_{e}=0.05$, than the value employed by them. As for the bulk energies comparing the two models for the neutron matter interactions, the values obtained using the modified Poschl-Teller potential are slightly larger than the ones using the $s$-wave part of AV18, although the relative difference is 2% at most.

We used the pairing functions of Eq. (II.3.2) to calculate the ground state energy of the cylindrical systems for $N=\{78,80,82,84,86\}$. Then we fitted the results to the functional form of Eq. (1), and we report the parameters $E_{0}^{\rm cyl}$ and $\lambda_{S}$ in Table 1. Ideally we would like to have $E_{0}^{\rm cyl}$ match the bulk value for every interaction strength, meaning that we can separate the ground state energy of the fermionic systems into a bulk component and a surface term. For both cold gases and low-density neutron matter, and most interaction strengths, the results are within the error bars. For the $s$-wave part of the AV18 model the agreement is quite good. The relative difference does not exceed 11%, and most values agree within error bars. It is worth pointing out that the values of $\lambda_{S}$ are negative for these systems due to the repulsive core of the interaction, see Fig. 2, a feature that is absent in the purely attractive potentials employed in the other cases. For the modified Poschl-Teller potential, the fitting procedure yielded larger errors. Also, the results for $-k_{F}a=1$ do not follow the trend, most probably due to the diluteness of the system. In Fig. 3 we compare the bulk energies of cold gases and neutron matter with the corresponding values of $E_{0}^{\rm cyl}$.

In Sec. II.2 we presented the potential interaction used for neutrons of antiparallel spins, and we set the interaction between particles of the same spin to zero. In doing so, we neglected the interaction of the $M=\pm 1$ triplet states. Previous QMC simulations of bulk low-density neutron matter Gezerlis and Carlson (2008), using a similar formalism to ours, found that, perturbatively, corrections for the artificial attraction in the $M=0$ triplet state account for 10% of the total energy in the $-k_{F}a=10$ case. The corrections become even lower for lower densities, such that in the range considered in this work they are of order of a few percent. Later calculations Gezerlis and Carlson (2010) compared results using the pure $s$-wave interaction with the AV4’Wiringa and Pieper (2002) potential, which yielded $\sim$ 7% difference for $-k_{F}a=10$, $\sim$ 1% for $-k_{F}a=5$, and essentially the same results for lower densities. These results in the bulk neutron matter systems justify our approach because, besides vanishing small corrections to the total energy as the density is lowered, one of our goals is to calculate the vortex excitation energy, which is an energy difference, and thus the corrections is expected to cancel.

III.1.2 Vortex excitation energy

The energy of the systems with a vortex line was calculated using the pairing functions of Eq. (II.3.3). The excitation energy was computed using the energy difference between those systems and the ground state of the cylinder. The results were averaged for $N=\{78,80,82,84,86\}$. Figure 4 shows the excitation energy for low-density neutron matter and cold atoms as a function of $k_{F}a$.

Although several results are within error bars, we can see that for $-k_{F}a=1$ the vortex excitation energy for cold gases and neutron matter (both models) is comparable. As the interaction strength increases, we can clearly see that the excitation energy is higher for cold gases systems compared to low-density neutron matter. The results for neutron matter, using both models, seem to be much less dependent on the interaction strength for this density regime. As was the case in the previous section, the errors associated with the modified Poschl-Teller potential for neutron matter are larger than the other two cases, however, it is still possible to see that the results for the two neutron matter models are close.

III.2 Density profiles

The density profile ${\cal D}(\rho)$ was calculated averaging the angular ($\varphi$) and axial ($z$) directions. We chose a normalization such that

[TABLE]

where the integral is over the volume $V=\pi{\cal R}^{2}{\cal L}$ of the cylinder. We show our results for the ground state density of the cylindrical container in Fig. 5. The results for cold gases, Fig. 5(a), follow a similar trend, with the exception of $-k_{F}a=5.0$, the largest interaction strength considered. Nonetheless, the Friedel oscillations are much smoother than in the neutron matter systems, Figs. 5(b) and 5(c). The results using the $s$-wave of the AV18 model show a very pronounced oscillation near $\approx$ 2.0$k_{F}^{-1}$ for $-k_{F}a=1.0$, and it is less intense for stronger interactions.

The hard wall condition introduces a characteristic density behavior close to it as it was discussed in Sec. II.1, and as seen in Fig. 5. We were able to separate two contributions to the ground state energy of the cylindrical systems, which we identified as bulk and surface terms, Sec. III.1.1. However, this analysis requires that there is a sufficient number of particles in the central region of the cylinder, away from the walls. To this end, we define the particle number $\eta(R)$ a distance $R$ from the $z$ axis,

[TABLE]

such that $\eta(R={\cal R})=N$. In Fig. 6 we plot $\eta(R)$ for cold gases and neutron matter systems using $N=84$ particles, which show essentially the same behavior, independently of the interaction strength. Figure 5 suggests that the hard walls affect the systems at $\rho\gtrsim 6.0k_{F}^{-1}$. As we can see in Fig. 6, $\eta(\sim 6.0k_{F}^{-1})\gtrsim$ 45, meaning that we have approximately this number of particles in the “bulk” portion of the cylinder. For systems with a vortex line this number is lower, $\approx$ 42. Previous QMC simulations of bulk properties have employed $N=$38 Carlson et al. (2003) and $N=$40 Forbes et al. (2011), hence the number of particles we have in the center of the cylinder is larger than in those bulk calculations.

In Figs. 7 and 8 we plot the density profiles for $-k_{F}a=\{0.5,1.0,2.0,3.3,5.0\}$ of the ground and vortex line states for cold atoms and neutron matter. We compared the density profiles of cold gases for $-k_{F}a=1.0$ and 2.0, Fig. 7, with the Bogoliubov-de Gennes calculations of Ref. Sensarma et al. (2006). They used a different geometry than ours, so to compare the results we changed their normalization to match our number of particles in the $\rho\leqslant 6.0k_{F}^{-1}$ region of the cylinder. Their results are closer to ours in the $-k_{F}a=1.0$ case, as expected. In the low-density neutron matter case, we compared our results for $-k_{F}a=3.3,5.0$ with the mean-field results of Ref. Yu and Bulgac (2003), see Fig. 8. In a similar fashion to what we did in the cold atoms case, we matched the normalizations to ensure the same number of particles in the $\rho\leqslant 24.5$ fm region.

A direct comparison of the density profiles for cold gases and neutron matter, or the two models we used for neutron matter, is difficult due to the different position of the oscillations in the profiles. However, a quantity of interest in both rotating superfluid cold gases systems and neutron matter is the density depletion at the vortex core, which depends only on the density near the origin. In the BCS limit the density should be close to the ground state one, while in the BEC regime the core should be completely depleted. We found that, in the cold gases case, the ratio of the density at $\rho=0$ for the system with a vortex line and the ground state of the cylindrical container decreases from 72% ($-k_{F}a=0.5$) to 53% ($-k_{F}a=5.0$), with the values 71%, 62%, and 59% for $-k_{F}a=1.0,2.0,3.3$, respectively. The mean-field calculation of Ref. Sensarma et al. (2006) finds a much higher density close to the BCS limit, 94% at $-k_{F}a=1.0$. However, their result for $-k_{F}a=2.0$ is comparable with ours, 60%. For low-density neutron matter using the modified Poschl-Teller potential, we found a density at the core of approximately 75%, 71%, 76%, and 78% of the ground state density for $-k_{F}a=1.0$, 2.0,3.3, and 5.0, respectively. The $s$-wave part of AV18 model for neutron matter yields 75% of the ground state density for $-k_{F}a=1.0,2.0$. For $-k_{F}a=5.0$ this ratio is 65%, close to the value of 60% of Ref. Yu and Bulgac (2003). For $-k_{F}a=3.3$ we see a small depletion, but that is an artifact of the density oscillations near the origin, so we chose not to include this interaction strength in the density at the core discussion.

IV Summary and outlook

In this work we compare properties of vortices in low-density neutron matter and cold atomic dilute Fermi gases. Our goal is not to show that they are identical, but rather to draw a parallel between their properties such that measuring quantities in cold Fermi gases can help to constrain properties of vortices in neutron matter, as previously done for the ground state Gandolfi et al. (2015).

Although the ground-state energies per particle of the bulk systems were lower for cold gases than for neutron matter, for a given $k_{F}a$, the difference becomes smaller as we move toward more dilute systems. This was the main motivation to expect that vortices in the low-density regime show a duality between cold gases and neutron matter. The excitation energy for the formation of a vortex line is comparable when the density is low enough. However, it is higher for cold gases than in neutron matter, so that must be accounted for when comparing the two systems.

We chose to analyze the density depletion at the vortex core, because it only depends on the density behavior close to the axis of the cylinder, away from the hard walls. Again, we found an agreement between the values for very low densities, although the density at the vortex core tends to remain close to 75% of the ground-state density for neutron matter, whereas we can clearly see it dropping from 72% to 53% for cold gases.

We found an excellent agreement when comparing the two models we employed for the neutron-neutron interactions. It seems remarkable that two potentials of completely different shapes, see Fig.(2), give us the same physical properties. However, the fact that they have the same scattering length and effective range is the key feature. This indicates that the low-energy limit of Eq. (3) is also valid for low-density neutron matter.

Our results can help to relate cold atom experiments with properties of low-density neutron matter. The extraction of bulk properties from experiments is extremely difficult when they employ harmonic traps. However, box-like traps Gaunt et al. (2013) have been successfully implemented in Bose systems, and they can help pave the way to determining the equation of state for cold gases. That, in turn, could be contrasted with Fig. 3 to constraint the low-density neutron matter equation of state. Vortices in fermionic gases on both BCS and BEC sides of the crossover, and also at unitarity, have been observed Zwierlein et al. (2005).

Our approach is valid for the low-density regime of neutron matter. However, it would be interesting to investigate vortex properties at higher densities. In Sec. III.1.1 we discussed possible corrections to account for our choice of neutron-neutron interaction potential. We showed that they would be small in the bulk case, thus justifying our approach, but they increase with the density. Instead of carefully including corrections, it seems more promising to consider realistic nuclear Hamiltonians. There are calculations using auxiliary-field diffusion Monte Carlo (AFDMC) Gandolfi et al. (2008, 2009) where bulk properties of neutron matter are calculated, at higher densities than in this present work, using realistic nuclear Hamiltonians. A possible extension of our work is to generalize the wave functions we presented by including spin correlations, and perform AFDMC simulations. The comparison of the results using both methods should enlighten how important spin correlations are when describing low-density neutron matter.

Acknowledgements.

We thank Alessandro Lovato for the useful discussions. This work was supported by the São Paulo Research Foundation (FAPESP) under the grant 2018/09191-7, and by the National Science Foundation under the grant PHY-1404405. The work of S.G. was supported by the NUCLEI SciDAC program, by the U.S. DOE under contract DE-AC52-06NA25396, by the LANL LDRD program, and by the DOE Early Career Research Program. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) SuperMIC and Stampede2 through the allocation TG-PHY160027, which is supported by National Science Foundation grant number ACI-1548562. This research also used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001. We also used resources provided by NERSC, which is supported by the US DOE under Contract DE-AC02-05CH11231.

Appendix A Equation (5)

The two-body problem of three-dimensional scattering with the modified Poschl-Teller potential, Eq. (4), can be solved analytically. At low-energies, we have an expression for the phase shift Flugge (1994),

[TABLE]

where $q=k/(2\mu)$. We can use the following relations Abramowitz and Stegun (1964),

[TABLE]

to cast the Eq. (A) in the form

[TABLE]

Approximating $\delta_{0}\approx-ka$ yields Eq. (5).