Structure Factors of The Unitary Gas Under Supernova Conditions

TL;DR

This paper calculates the static structure factors of the unitary gas at various temperatures and fugacities using lattice field theory, providing insights beyond traditional virial expansion methods relevant for supernova physics.

Contribution

It introduces a lattice field theory approach to compute structure factors of the unitary gas at high fugacities, surpassing virial expansion limitations.

Findings

Qualitative differences in structure factors at high fugacity compared to virial predictions

Calculations applicable to neutron gases in supernova environments

Potential to improve neutrino scattering rate estimates in astrophysics

Abstract

We compute with lattice field theory the vector and axial static structure factors of the unitary gas for arbitrary temperature above the superfluid transition and for fugacities 0.1 < z < 1.0. Using the lattice formulation, we calculate beyond the validity of the virial expansion, a commonly used technique in many-body physics. We find qualitative differences in the behavior of the structure factors at high fugacity compared to the predictions of the virial expansion. Due to the large scattering length of neutrons, we expect the unitary gas structure factors to approximate the structure factors of hot neutron gases, and we therefore expect our calculations to be useful in supernova simulations, where neutron gas structure factors are needed to compute in-medium neutrino-neutron scattering rates.

Figures

Tables

Table 1. Table 1 : Three sets of parameters used in finite volume, space and time continuum extrapolations.

Volume extrapolation	Space continuum extrapolation	Time continuum extrapolation
$\Delta x=0.4,\Delta t=0.01,N_t=40,$	$\Delta x=0.4\frac{5}{N_x},\Delta t=0.01,N_t=40,$	$\Delta x=0.4,N_x=7,\Delta t=0.01\frac{40}{N_t},N_t=40,$
$N_x=5,7,9$	$N_x=5,7,9$	$N_t=40,60,80$

Equations

$$\frac{d\Gamma }{d\cos \theta }=\frac{G_F^2{E}_{\nu }^2}{4{\pi }^2}\left[c_V^2(1+\cos \theta )S_V(q)+c_A^2(3-\cos \theta )S_A(q)\right]$$

$$S_V(q)$$

$$S_A(q)$$

$$H-\mu N=\underset{K}{\underbrace{\sum _{\sigma ;xy}{\psi }_{\sigma ;x}^{†}{k}_{xy}{\psi }_{\sigma ;y}}}\underset{V}{\underbrace{-\mu \sum _{\sigma ;x}{\psi }_{\sigma ;x}^{†}{\psi }_{\sigma ;x}-\frac{G}{\Delta x^3}\sum _x{\psi }_{1;x}^{†}{\psi }_{1;x}{\psi }_{2;x}^{†}{\psi }_{2;x}}}.$$

$$k_{xy}=\frac{1}{2M\Delta x^2}\sum _p{p}^2{\text{e}}^{-ip\cdot (x-y)},$$

$$H=\int{d^{3}x\Big{[}\psi^{\dagger}_{\sigma}\big{(}-\frac{\nabla^{2}}{2M}-\mu\big{)}\psi_{\sigma}-\frac{G}{2}(\psi_{\sigma}^{\dagger}\psi_{\sigma})^{2}\Big{]}}~{},$$

$$Z=\mathrm{Tr}(e^{-\beta (H-\mu N)})=\mathrm{Tr}\left(\prod _{t=1}^{N_t}{e}^{-\Delta tK}{e}^{-\Delta tV}\right)+\mathcal{O}(\Delta t^2),$$

$$e^{-\Delta tV}=\prod _x\exp \left(\Delta t\mu (n_{1x}+n_{2x})+\frac{\Delta tG}{\Delta x^3}{n}_{1x}{n}_{2x}\right)=\prod _x\frac{1}{f_0}\int d{A}_x{e}^{-\frac{1}{\hat{g}}(\cosh A_x-1)}{e}^{(n_{1x}+n_{2x})(A_x+\hat{\mu })},$$

$$f_{\alpha }(\hat{g})\equiv {\int }_{-\infty }^{\infty }dA\exp \left(-\frac{\cosh A-1}{\hat{g}}+\alpha A\right).$$

$$\mathrm{Tr}\left(e^{-\sum {\psi }_x^{†}(M_1{)}_{xy}{\psi }_y}{e}^{-\sum {\psi }_x^{†}(M_2{)}_{xy}{\psi }_y}\cdots e^{-\sum {\psi }_x^{†}(M_N{)}_{xy}{\psi }_y}\right)=\det \left(1+e^{-M_N}\cdots e^{-M_1}\right)$$

$$Z\approx \mathrm{Tr}\left(\prod _{t=1}^{N_t}{e}^{-\Delta tK}{e}^{-\Delta tV}\right)=f_0^{-N_s{N}_t}\int DA{e}^{-S_g(A)}\det {\left[1+e^{N_t\hat{\mu }}{B}^{-1}C(A_{N_t-1})\cdots B^{-1}C(A_0)\right]}^2$$

$$B_{x{x}^{\prime }}=\sum _p{e}^{-ip\cdot (x-x^{\prime })}{e}^{\hat{\gamma }\frac{p^2}{2}},C(A_t{)}_{x{x}^{\prime }}={\delta }_{x{x}^{\prime }}{e}^{A_{x,t}},\hat{\gamma }\equiv \frac{\Delta t}{M\Delta x^2},$$

$$S$$

$$\int DAD\hat{\psi }D{\hat{\psi }}^{†}{e}^{-S[A,{\psi }^{†},\psi ]}=\int DA{e}^{-S_g[A]}\det (D{)}^2,$$

$$\det D=\det B^{N_t}\det \left[1+e^{N_t\hat{\mu }}{B}^{-1}C(A_{N_t-1})...B^{-1}C(A_0)\right]=\det (1+U).$$

$$\begin{array}{rl}& {\int }_{-\infty }^{\infty }dA\exp \left(-\frac{\cosh A-1}{\hat{g}}\right)\exp [-{\hat{\psi }}^{†}{e}^{A+\hat{\mu }}\hat{\psi }]={\int }_{-\infty }^{\infty }dA\exp \left(-\frac{\cosh A-1}{\hat{g}}\right)[1-{\hat{\psi }}^{†}{e}^{A+\hat{\mu }}\hat{\psi }+\frac{1}{2}({\hat{\psi }}^{†}{e}^{A+\hat{\mu }}\hat{\psi }{)}^2]\\ & =f_0-f_1{\hat{\psi }}^{†}{e}^{\hat{\mu }}\hat{\psi }+\frac{1}{2}{f}_2({\hat{\psi }}^{†}{e}^{\hat{\mu }}\hat{\psi }{)}^2=f_0\exp \left[\frac{f_1}{f_0}{\hat{\psi }}^{†}{e}^{\hat{\mu }}\hat{\psi }+\left(\frac{f_2}{f_0}-\frac{f_1^2}{f_0^2}\right)\frac{1}{2}({\hat{\psi }}^{†}{e}^{\hat{\mu }}\hat{\psi }{)}^2\right].\end{array}$$

$$S_f=\sum _{\sigma ;x,y,t}{\psi }_{\sigma ;x,t}^{†}{B}_{xy}{\psi }_{\sigma ;y,t}-e^{\hat{\mu }}\frac{f_1}{f_0}\sum _{\sigma ;x,t}{\psi }_{\sigma ;x,t+1}^{†}{\psi }_{\sigma ;x,t}-e^{2\hat{\mu }}\left(\frac{f_2}{f_0}-\frac{f_1^2}{f_0^2}\right)\sum _{x,t}{\psi }_{1;x,t+1}^{†}{\psi }_{1;x,t}{\psi }_{2;x,t+1}^{†}{\psi }_{2;x,t}.$$

$$S_f=$$

$$-e^{2\hat{\mu }}\left(\frac{f_2}{f_0}-\frac{f_1^2}{f_0^2}\right){\int }_{-\pi }^{\pi }{d}^4{p}_1{d}^4{p}_2{d}^4{p}_3{d}^4{p}_4{\delta }^4(p_1+p_3-p_2-p_4)e^{i{\omega }_1+i{\omega }_3}{\psi }_1^{†}(p_1){\psi }_1(p_2){\psi }_2^{†}(p_3){\psi }_2(p_4),$$

$$\text{vbox}\text{fmfgraph}*(65,50)\text{fmfleft}i1\text{fmfright}o1\text{fmf}fermioni1,o1=\frac{1}{{\text{e}}^{\hat{\gamma }\vec{p}\cdot \vec{p}/2}-{\text{e}}^{i\omega +\hat{\mu }+\ln f_1/f_0}}\equiv G(\omega ,\vec{p})$$

$$\begin{gathered}\fmfgraph*(65,50)\fmfleft{i1,i2}\par\par\fmfright{o1,o2}\par\fmflabel{$1$}{i1}\fmflabel{$2$}{o1}\fmflabel{$3$}{i2}\fmflabel{$4$}{o2}\par\par\fmf{fermion}{i1,v,o1}\fmf{fermion}{i2,v,o2}\end{gathered}=\text{e}^{2\hat{\mu}}\big{(}\frac{f_{2}}{f_{0}}-\frac{f_{1}^{2}}{f_{0}^{2}}\big{)}\text{e}^{i\omega_{1}+i\omega_{3}}$$

$$-\mathcal{A}=\frac{1}{\Big{[}e^{2\hat{\mu}}\big{(}\frac{f_{2}}{f_{0}}-\frac{f_{1}^{2}}{f_{0}^{2}}\big{)}e^{iE}\Big{]}^{-1}-L(E)}~{},$$

$$L(E)={\int }_{-\pi }^{\pi }\frac{d{q}_0}{2\pi }{\int }_{-\pi }^{\pi }\frac{d^{3}q}{(2\pi {)}^3}G(q_0+E/2,\vec{q})G(-q_0+E/2,-\vec{q})$$

$$\frac{f_2{f}_0}{f_1^2}-1=\frac{1}{L(0)},\text{with}L(0)={\int }_{-\pi }^{\pi }\frac{d^{3}q}{(2\pi {)}^3}\frac{1}{e^{\hat{\gamma }{q}^2}-1}.$$

$$\frac{MG}{\Delta x}=\underset{\hat{\gamma }\to 0}{\lim }\frac{1}{\hat{\gamma }L(0)}=5.14435\dots ,$$

$$-\mathcal{A}=\text{parbox}85.35826pt\text{fmfgraph}\ast (100,70)\text{fmfleft}i1,i2\text{fmfright}o1,o2\text{fmf}fermioni1,v,o1\text{fmf}fermioni2,v,o2+\text{parbox}85.35826pt\text{fmfgraph}\ast (100,70)\text{fmfleft}i1,i2\text{fmfright}o1,o2\text{fmf}fermioni1,v\text{fmf}fermioni2,v\text{fmf}fermion,left=1v,v2\text{fmf}fermion,right=1v,v2\text{fmf}fermionv2,o1\text{fmf}fermionv2,o2+\text{parbox}85.35826pt\text{fmfgraph}\ast (100,70)\text{fmfleft}i1,i2\text{fmfright}o1,o2\text{fmf}fermioni1,v\text{fmf}fermioni2,v\text{fmf}fermion,left=1v,v2\text{fmf}fermion,right=1v,v2\text{fmf}fermion,left=1v2,v3\text{fmf}fermion,right=1v2,v3\text{fmf}fermionv3,o1\text{fmf}fermionv3,o2+...$$

$$\begin{split}\langle n\rangle&=\frac{1}{V}\frac{\partial}{\partial\mu}(T\log Z)\\ &=\frac{1}{\Delta x^{3}}\frac{1}{N_{s}}\frac{1}{Z}\int{DA~{}\text{e}^{-S_{g}(A)}\det(\mathbf{1}+U(A))^{2}\operatorname{Tr}\Big{[}\mathbf{1}-(\mathbf{1}+U(A))^{-1}\Big{]}}\\ &=\frac{1}{\Delta x^{3}}\frac{1}{N_{s}}\langle\operatorname{Tr}\Big{[}\mathbf{1}-(\mathbf{1}+U)^{-1}\Big{]}\rangle\,,\end{split}$$

$$\begin{array}{rl}{S}_V(q)& =\sum _x\Delta x^3{\text{e}}^{-iq\cdot x}⟨\delta n(x,0)\delta n(0,0)⟩,\\ S_A(q)& =\sum _x\Delta x^3{\text{e}}^{-iq\cdot x}⟨\delta s_z(x,0)\delta s_z(0,0)⟩.\end{array}$$

$$\begin{array}{rl}\delta n(x,t)& =n(x,t)-⟨n⟩,\\ \delta s_z(x,t)& =s_z(x,t)-⟨s_z⟩.\end{array}$$

$$\begin{array}{rl}\Delta x^6⟨n_1(x,0)n_1(0,0)⟩& =\frac{1}{Z}\int DA{\text{e}}^{-S_g(A)}\det (1+U(A){)}^2\left[{\left(1-(1+U(A){)}^{-1}\right)}_{00}{\left(1-(1+U(A){)}^{-1}\right)}_{xx}\right]\\ & +\frac{1}{Z}\int DA{\text{e}}^{-S_g(A)}\det (1+U(A){)}^2\left[{\left(1-(1+U(A){)}^{-1}\right)}_{0x}{\left(1+U(A)\right)}_{x0}^{-1}\right]\\ & \\ \Delta x^6⟨n_1(x,0)n_2(0,0)⟩& =\frac{1}{Z}\int DA{\text{e}}^{-S_g(A)}\det (1+U(A){)}^2\left[{\left(1-(1+U(A){)}^{-1}\right)}_{00}{\left(1-(1+U(A){)}^{-1}\right)}_{xx}\right].\end{array}$$

$$S_V(q)/⟨n⟩$$

Full text

I Introduction

The unitary gas, a gas of spin $1/2$ fermion with short range interactions tuned so a two-body bound state exists at threshold, is a classic example of a conformal system. The challenge in studying it stems from the fact that it is a strongly coupled theory for which no perturbative method applies. At the same time, a number of physical systems in Nature, varying over a large range of energy scales, are very close to being a unitary gas. Besides dilute atomic fermionic gases in atomic traps, the dilute and warm neutron gas, as found in core collapse supernova, is close to being a unitary gas. Indeed, the conditions required for a neutron gas to be close to the unitary gas is that the typical interparticle distance $l$ should lie between the scattering length $a$ and the range of the forces $R$, while, at the same time, the thermal wavelength $\lambda_{T}$ should be larger than $R$, so that the details of the potential are not probed during collisions. For a neutron gas where the scattering length is $a\simeq-23\,{\tt fm}$ and the range of the nuclear forces $R$ set by the inverse pion mass , $R\approx 1/m_{\pi}\approx 1.3\,{\tt fm}$, the neutron gas will approximate well a unitary gas over the range of densities and temperatures $0.01n_{0}\leq n\leq 0.1n_{0}$ ($n_{0}\approx 0.16/\,{\tt fm}^{3}$ is the density of nuclear matter) and $T\leq 10\,{\tt MeV}$. This is similar to the range of densities and temperatures found in the neutrinosphere of core collapse supernovae Lin and Horowitz (2017); Horowitz and Schwenk (2006a, b), the region where neutrinos decouple from matter during the explosion, and which is essential to the modeling of the core collapse supernova process Burrows (2013); Mezzacappa (2005).

The interactions between the neutrinos and the neutron matter are controlled not only by the cross section of neutrinos scattering off individual neutrons, but also by many-body effects. As a particularly dramatic demonstration of this, recent 3D supernova simulations have shown that small changes in the scattering rate, of the order that could be caused by in-medium effects Melson et al. (2015), can turn a dud into an explosion. The in-medium differential scattering rate of neutrinos due to neutral current scattering processes in the elastic (zero energy transfer) approximation is given by Horowitz and Schwenk (2006c):

[TABLE]

where $E_{\nu}$ is the energy of the neutrino, $\cos\theta$ is the scattering angle and $q=2E_{\nu}\sin\theta/2$ is the momentum transfer to the medium in which the neutrino is scattering. The functions $S_{V}(q)$ and $S_{A}(q)$ are the vector and axial static structure factors, respectively, and they encode the properties of the medium. In the non-relativistic limit relevant here 111The conditions of the neutrinosphere, $T\sim 5-10\text{ MeV}$ and $n\sim 10^{-3}-10^{-1}n_{0}$ are such that the neutrons composing the medium can be approximated by a non-relativistic spin-1/2 field. In the non-relativistic limit, neutrinos couple to neutron density and spin., the structure factors are defined as Fourier transforms of equal time correlation functions

[TABLE]

where $n$ is the density of particles, and $s_{z}$ is the density of z component of spin. In a supernova, these structure factors are computed in a high temperature medium composed of primarily neutrons; we calculate here $S_{V}$ and $S_{A}$ in a hot medium composed of unitary fermions as a first step toward calculating the structure factors of the neutron gas.

Various properties of the unitary gas have been computed in the past, with particular interest in the low temperature regime. Thermodynamic functions in the have been calculated by several groups Gilbreth and Alhassid (2013); Jensen et al. (2018); Bulgac et al. (2008); Drut et al. (2012); Magierski et al. (2009); Wlazłowski et al. (2013). These calculations are done in the vicinity of the superfluid phase transition and in the postulated “pseudogap” regime, both of which occur at low temperature. The contact has also been calculated by several groups Jensen et al. (2019); Drut et al. (2011); Goulko and Wingate (2016); Rossi et al. (2018), and the vector static structure factor has been calculated by relatively fewer groups, with a zero temperature Quantum Monte Carlo calculation performed in Hoinka et al. (2013).

Due to the unique conditions of the neutrinosphere, we explore here a higher temperature regime than in previous studies; as a result we probe novel regions of parameter space of the unitary gas. Furthermore, we present here the first non-perturbative calculation of both the axial and vector structure factors of the unitary gas at the thermodynamic parameters found in the neutrinosphere. In particular, we calculate the structure factors of the unitary gas beyond the domain of validity of the commonly used virial expansion, which has been used to calculate structure factors perturbatively in the past Bedaque et al. (2018) 222The virial expansion has also been used to compute dynamic properties of the unitary gas. For example the shear and bulk viscosities have very recently been computed in Hofmann (2019)..

This paper is organized as follows: In Section II we describe the methods we used to numerically simulate the unitary gas. In Section III we describe how we tune lattice parameters of the hamiltonian to satisfy the unitary condition. In Section IV we derive the observables of interest. In Section V we show the results of our lattice calculation and in Section VI we state our conclusions.

II Formalism

We will use the following lattice hamiltonian to describe the unitary gas:

[TABLE]

Here $\psi_{\sigma;x}$ annihilates a fermion with spin $\sigma$ at site $x$ on a 3D spatial lattice, the fermionic fields are normalized such that $\{\psi_{\sigma;x},\psi^{\dagger}_{\sigma^{\prime};x^{\prime}}\}=\delta_{\sigma\sigma^{\prime}}\delta_{xx^{\prime}}$ and the matrix $k_{xy}$ is given by

[TABLE]

where $p$ is a lattice momentum. The non-local hopping term generates the correct, continuum dispersion relation, instead of a lattice approximation that would be generated if a local hopping term had been used. The non-local structure of the hamiltonian will not cause any difficulty in the simulations as this term is treated in momentum space. Of course, the hamiltonian in Eq. (2) is a discretized version of

[TABLE]

which, for a properly renormalized $G$, describes non-relativistic fermions at unitarity. As we will show latter, if we regularize the theory using the lattice discretization in Eq. (2), the unitarity point is given by $G\approx 5.14435\Delta x/M$.

Our object of study is the partition function (and its derivatives with respect to local spacetime sources)

[TABLE]

where $\Delta t=\beta/N_{t}$ 333That this trotterization is correct to $\mathcal{O}(\Delta t^{2})$ can be seen by noting that

$\operatorname{Tr}(\prod_{t=1}^{N_{t}}{e^{-\Delta tK}e^{-\Delta tV}})=\operatorname{Tr}(\prod_{t=1}^{N_{t}}{e^{-\frac{\Delta t}{2}K}e^{-\Delta tV}}e^{-\frac{\Delta t}{2}K})$ by cyclicity of the trace.. The numerical Monte Carlo method we will use is very similar to the one in White et al. (1989); Blankenbecler et al. (1981). The four-fermion interaction can be written in terms of fermion bi-linears by noting that

[TABLE]

with $n_{\sigma}=\psi^{\dagger}_{\sigma}\psi_{\sigma}$, provided $\mu,G$ are related to $\hat{\mu},\hat{g}$ through $e^{\mu\Delta t}=\frac{f_{1}}{f_{0}}e^{\hat{\mu}}$ and $e^{\frac{G\Delta t}{\Delta x^{3}}}={\frac{f_{0}f_{2}}{f_{1}^{2}}}$, where the functions $f(\hat{g})$ are

[TABLE]

Finally, using the identity

[TABLE]

for fermionic operators $\psi^{\dagger}_{x},\psi_{y}$, we see that

[TABLE]

where

[TABLE]

$N_{s}$ is the number of spatial sites, and $DA=\prod_{xt}dA_{xt}$.

The same partition function, up to a constant multiplicative factor, can also be derived from the path integral over the Euclidean lattice action:

[TABLE]

where we define the fermion matrix $D_{xt,x^{\prime}t^{\prime}}=B_{xx^{\prime}}\delta_{tt^{\prime}}-e^{A_{xt}+\hat{\mu}}\delta_{t+1t^{\prime}}$ (we implement anti-periodic boundary conditions in time by taking $\delta_{N_{t},N_{t-1}}\equiv-\delta_{0,N_{t-1}}$.) Indeed, the integration over Grassmann variables $\hat{\psi}$ leads to

[TABLE]

where $S_{g}(A)\equiv\frac{1}{\hat{g}}\sum_{x,t}{(\cosh(A_{x,t})-1)}$ and $D$ is the fermion matrix. The determinant $\det D$ can be written, by repeated use of the Schur complement identity, as

[TABLE]

where $U\equiv e^{N_{t}\hat{\mu}}B^{-1}C(A_{N_{t}-1})...B^{-1}C(A_{0})$. Therefore, the expression for the derived from the action  Eq. (12) is equivalent to the one derived from the Hamiltonian in Eq. (9).

A few words about the Monte Carlo calculation of expression in Eq. (9) (or, equivalently, of expression Eq. (12)) are in order. Firstly, since $D$ is a real matrix $\det D$ is real, $\det D^{2}\geq 0$ and the theory is sign-problem free. Secondly, the matrix $B^{-1}$ is independent of $A_{xt}$ and can be computed once at the beginning of the calculation and used thereafter. The matrices $C(A_{t})_{xx^{\prime}}$, which depend on $A_{xt}$, are diagonal and their multiplication has a computational cost of order ${\cal O}(N_{s}^{2})$, where $N_{s}$ is the number of spatial sites (instead of the ${\cal O}(N_{s}^{3})$ cost for dense matrices). Finally, we use the Hybrid Monte Carlo algorithm of Alexandru et al. (2018) to sample fields.

III Unitarity Condition

To obtain the unitary gas in the continuum and infinite volume limits, we tune the parameters of the lattice theory such that there is a zero energy two-particle bound state in the ${}^{1}S_{0}$ scattering channel at $\mu=0$. The binding energy of the two-particle state can be determined from the position of the appropriate pole in the scattering amplitude. For this model, conservation of particle number allows us compute exactly the scattering amplitude by summing over all Feynman diagrams that contribute to two-particle scattering processes.

To compute the Feynman rules we start from the action Eq. (12), and integrate over the auxiliary fileds:

[TABLE]

The action in terms of Grassmann variables is then:

[TABLE]

For an infinite lattice in both space and time, the momentum space action reads

[TABLE]

leading to the Feynman rules in Fig. 1.

In terms of diagrams, the scattering amplitude is given by the sum of graphs shown in Fig. 2. Choosing kinematics $(E/2,\vec{p})$ and $(E/2,-\vec{p})$ on the external legs, the scattering amplitude has the analytic expression

[TABLE]

where

[TABLE]

is the fermion-fermion loop iterated in the bubble sum. To work at $\mu=0$, using the map derived above, we set $\hat{\mu}=\log\tfrac{f_{0}}{f_{1}}$ and the unitarity condition, that the amplitude has a pole at $E=0$ at zero density, becomes:

[TABLE]

Given any $\hat{\gamma}=\Delta t/(M\Delta x^{2})$, this condition fixes the lattice coupling $\hat{g}$. In principle we could enforce this condition for any value of $\Delta x$ and $\Delta t$ and then choose any sensible path to take the limit $\Delta t\to 0$, $\Delta x\to 0$. In practice we find it advantageous to fix $\Delta x$, fix $G$ from the unitarity condition:

[TABLE]

and take the limit $\Delta t\to 0$ while $\hat{g}$ is fixed by $\frac{\Delta tG}{\Delta x^{3}}=\log(\frac{f_{2}f_{0}}{f_{1}^{2}})$. This has the advantage that for fixed $\Delta x$, the discretized partition function has ${\cal O}(\Delta t^{2})$ errors, since for this parameter choice the partition function derived from the action is the same as the approximation introduced in Eq. (5). This is the condition we use to fix couplings in Monte Carlo calculations. Of course, a subsequent spatial continuum limit $\Delta x\to 0$ is required to get the final results.

IV Observables

We compute three observables: the density, the vector structure factor and the axial structure factor. The density is defined through the standard thermodynamic relation, and has the path integral expression

[TABLE]

where $U(A)=e^{N_{t}\hat{\mu}}B^{-1}C(A_{N_{t}-1})...B^{-1}C(A_{0})$. Our results are derived using cubic boxes with $N_{x}$ sites in each dimension, so the total number of sites in a time slice is $N_{s}=N_{x}^{3}$. The vector and axial structure factors are defined as Fourier transforms of equal time correlation functions of the density $n=n_{1}+n_{2}$ and spin $s_{z}=n_{1}-n_{2}$, respectively:

[TABLE]

Here we have defined the Heisenberg operators

[TABLE]

It is possible to derive path integral expressions of $S_{V}(q)$ and $S_{A}(q)$ by taking derivatives with respect to local sources. Coupling the theory to a spacetime dependent source on the $t=0$ timeslice, then setting the source to a constant after taking derivatives, one derives the following path integral expressions:

[TABLE]

Since the action is invariant under $SU(2)$ spin rotations on the lattice (and in the continuum) $\langle n_{1}(x,0)n_{1}(0,0)\rangle=\langle n_{2}(x,0)n_{2}(0,0)\rangle$ and $\langle n_{1}(x,0)n_{2}(0,0)\rangle=\langle n_{2}(x,0)n_{1}(0,0)\rangle$, therefore the two correlation functions above are the only two independent density-density correlation functions. We compute these in our lattice calculations, then form the necessary linear combinations to compute $S_{V}$ and $S_{A}$. We also take advantage of translation invariance and hypercubic rotation invariance by summing over all such transformations in our calculation. This process reduces stochastic fluctuations of observables.

V Results

Since the interactions in a unitary gas do not have any scale (the scattering length is infinite) the static structure factors $S_{V,A}$ of the unitary gas are function of the momentum $q$, chemical potential $\mu$, temperature $T$ and mass $M$. The dimensionless quantity $S_{V,A}/n$ (where $n$ is the density of the system) is therefore a dimensionless function of only two variables

[TABLE]

with $q_{T}=\sqrt{6MT}$ and $z=e^{\mu/T}$. The presence of a lattice introduces three unphysical scales, the finite volume and the lattice spacings $\Delta t$ and $\Delta x$. However, we will demonstrate below that our calculations are performed at large enough volumes and small enough $\Delta t,\Delta x$ that the dependence of $F_{V,A}$ on $q_{T}\Delta x$ and $T\Delta t$ can be neglected.

Since the supernova collapse simulations are an important motivation for our calculations, we will concentrate on a subset of parameters relevant for the neutrino propagation on the neutrinosphere, namely, densities in the range of $0.001n_{0}<n<0.03n_{0}$ (where $n_{0}$ is the density of nuclear matter) and $3\,{\tt MeV}<T<10\,{\tt MeV}$, corresponding to fugacities in the $0.01\lesssim z\lesssim 3$ range. This range of values is motivational only. A realistic evaluation of structure factors to be used in supernova modeling should include a more realistic neutron-neutron interaction (with a non-zero scattering length and effective range) as well as a small amount of protons and light nuclei. This corrections will be investigated in a future study.

In order to explore the uncertainties arising from the finite volume and lattice spacing employed in our calculations we first compute the structure factors for two values of the fugacity, $z=0.1$ and $z=1.0$ using three sets of parameters (and $M=1$, $T=4.13$) shown in Table 1.

We use lattices with odd number of points in the spatial direction because such lattices converge to the infinite volume limit faster than those with even number of sites. For nuclear astrophysics applications it is useful to think of the systems of units we use as corresponding to $\Delta x=1=5\,{\tt fm},M=1=938\,{\tt MeV},\Delta t=0.01=1.19\,{\tt fm}$. The results we find are shown on Fig. 3 and Fig. 4. They show that the calculations performed on the smallest volume and coarsest lattices are already sufficient to determine the infinite volume and continuum results, within the statistical errors. The only exception are the results for the densest systems in our range ($z=1.0$). Only the two finest lattice results for $S_{A}$ ($\Delta x=0.4\times 5/7,0.4\times 5/9$), agree well within the statistical uncertainties while the results $S_{V}$ are different for all three values of $\Delta x$. Even in this case, however, signs of convergence are clear, with the difference between the $\Delta x=0.4$ and $\Delta x=0.4\times 5/7$ results being much smaller than between $\Delta x=0.4\times 5/7$ and $\Delta x=0.4\times 5/9$. Due to this dependence on $\Delta x$, we assign a systematic uncertainty of $5\%$ on the results for the vector structure factors at the largest values of $z$.

The asymptotic limit of the structure factors are known on general grounds:

[TABLE]

Since we are limited by the lattice cutoff to values of $q<\pi/\Delta x$ we can see that this limit is reached in less dense system ($z=0.1$) but larger values of the momentum would be needed to verify it at more dense systems ($z=1.0$), as can be seen on Fig. 5.

In the same figure we compare our results to some known analytic results. The $S_{V,A}(q\rightarrow 0)$ limits are determined by thermodynamics

[TABLE]

These have been computed up to fourth order in the virial expansion Lin and Horowitz (2017) and are shown in blue. The second and fourth order virial result is shown as small segments on the vertical axis. The second order in the virial expansion of a free gas can also be easily computed and is shown in Fig. 5 as a grey dotted line. Virial results for non-zero value of the momentum are not yet available. However, a calculation up to second order in $z$ using a pseudo-potential interaction (tuned to unitarity) and finite values of $q$ is available Bedaque et al. (2018) and shown as dashed lines in Fig. 5 444The main motivation of the calculation in Bedaque et al. (2018) was the computation of dynamical structure factors which cannot be computed with Monte Carlo methods of the kind discussed in this paper. . As expected, all virial calculations agree well with our data at small ($z\lesssim 0.5$) but differ significantly at larger values of $z$.

Finally, we have compared with a recent calculation by Jensen et al Jensen et al. (2018) of the thermodynamics of the unitary gas, including the $S_{A}(q\rightarrow 0)$ limit. The authors explore the parameter space $0.1<T/T_{F}<0.4$, where $T_{F}=(3\pi^{2}n)^{2/3}/2M$. We extended our calculations to lower temperatures to be able to compare with their values. Jensen et al. report $S_{A}(q\rightarrow 0)/S_{0}=0.403(1)$ at $T/T_{F}=0.353$ (where $S_{0}=\frac{3n}{2}\frac{T}{T_{F}}$) on $9^{3}$ lattices, while we obtain $S_{A}(q\rightarrow 0)/S_{0}=0.39(1)$ at $T/T_{F}=0.371(5)$. Our results therefore appear consistent with those of Jensen et al. Another calculation Wlazłowski et al. (2013) reports values about $20\%$ higher.

Our main results are summarized in Fig. 6. The error bar include the statistical errors but it should be kept in mind that an additional systematic for the vector structure factor reaching up $5\%$ for the highest $z=1.0$ should be added to the uncertainties.

VI conclusions

We computed the static structure factors of the unitary gas in the $0.6T_{F}\lesssim T\lesssim 6T_{F}$ using Monte Carlo methods, corresponding to fugacities in the $0.1\lesssim z\lesssim 1.0$ range. The unitary gas at low densities is a reasonable model for a dilute neutron gas and this range of fugacities is the one relevant to the neutrinosphere (the region around a proto-neutron star where neutrinos decouple) in core collapse supernova. The structure factors are the relevant many-body information to describe the interactions between the neutrinos and the neuron matter that drives the explosion. We compared our lattice calculations of the static structure factors with an $\mathcal{O}(z^{2})$ calculation using the pseudopotential approach. We find reasonable agreement with the second order virial expansion at $z=0.1$, while we find qualitative differences between the lattice and virial expansion at $z=1.0$. In particular, the dramatic bump in the vector structure factor at small $q/q_{T}$ predicted by the virial expansion is in fact not present.

We demonstrate that, for all densities explored, both the infinite volume and the $\Delta t\rightarrow 0$ limits are very well behaved, with negligible variation in both structure factors. At small density, the $\Delta x\rightarrow 0$ limit also yields negligible variation in the structure factors. At $z=1.0$, the densest system we consider, we see some variation in the structure factors at the coarsest lattice spacings. We expect, however, that our calculations are within $5\%$ of the exact result at our finest lattice spacing.

We also leave for future work an extension of this calculation to hot neutron matter that will require a more realistic interaction between neutrons, at least one leading to finite scattering lengths and non-zero effective range. We also leave for future work a more careful $\Delta x\rightarrow 0$ extrapolations. Such extrapolations require calculations on larger lattices, which are more expensive, but there is no fundamental difficulty in doing so.

Acknowledgements.

A. A. is supported in part by the National Science Foundation CAREER grant PHY-1151648 and by U.S. Department of Energy under Contract No. DE-FG02-95ER-40907. A.A. acknowledges the hospitality of the University of Maryland where part of this work was performed. P. B. and N.W. are supported by the U.S. Department of Energy under Contract No. DE-FG02-93ER-40762. P.B. and N.W. acknowledge the hospitality of The George Washington University where part of this work was performed.