Spin Polarized Non-Relativistic Fermions in 1+1 Dimensions

We study by Monte Carlo methods the thermodynamics of a spin polarized gas of non-relativistic fermions in 1+1 dimensions. The main result of this work is that our action suffers no significant sign problem for any spin polarization in the region relevant for dilute degenerate fermi gases. This lack of sign problem allows us to study attractive spin polarized fermions non-perturbatively at spin polarizations not previously explored. For some parameters values we verify results previously obtained by methods which include an uncontrolled step like complex Langevin and/or analytical continuation from imaginary chemical potential. For others, larger values of the polarization, we deviate from these previous results.


I. INTRODUCTION
Strongly coupled dilute fermi gases are rich systems which stand at the crossroads of various fields of physics. In the realm of nuclear physics, it is believed that the properties of dilute fermi gases with large scattering lengths are related to the properties of neutron rich matter at the edge of a neutron star [1]. In the context of core collapse supernovae, the dynamic properties of dilute neutron gases must be analyzed to understand the propagation of energy which drive the explosion [2,3]. In the atomic physics community, dilute fermi gases at unitarity have been studied to in the context of the BEC-BCS crossover [4]. Spin polarizing these gases adds a layer of richness to observable phenomena. For instance, unequal density for the two spin states may either lead to an instability of Cooper pairing [5] or the formation of the exotic superconducting phase predicted by Larkin, Ovchinnikov, Fulde, and Ferrell [6,7], characterized by formation of Cooper pairs which spontaneously break translation invariance. High magnetic fields in neutron stars can destroy S-wave pairing in neutron matter through spin polarization, which has observable consequences on spectroscopic data gathered from rapidly rotating neutron stars [8].
All these systems are strongly coupled and to investigate their properties requires non-perturbative methods. One general approach to analyzing strongly coupled systems with controllable systematic uncertainties are ab initio Monte Carlo calculations. This approach runs into difficulties for spin imbalanced non-relativistic systems due to the sign problem which causes statistical uncertainties to grow exponentially as the thermodynamic and zero temperature limits are taken. The source of the sign problem in fermionic systems are sign oscillations of the fermion determinant.
A few suggestions for dealing with the sign problem in spin polarized non-relativistic systems have been put forward in the past. One is to compute observables at imaginary chemical potential, which renders the fermion determinant positive definite, and hence eliminates the sign problem and analytically continue them to real chemical potentials. The process of analytic continuation from restricted, noise data is necessarily done with the help of an ansatz which introduces an uncontrolled source of error. Still, mean field results suggest that interesting features of the phase diagram lie within the region where extrapolation works [9]. For the one dimensional case, this avenue has been tested numerically [10]. Another approach is the Complex Langevin method where the configuration space is augmented by complexifying the field variables and stochastic process is used to sample configurations in this complex space [11,12]. Complex Langevin has also been used to analyze non-relativistic fermions with a mass imbalance, which is a system that also suffers from the sign problem [13]. The potential problem of complex Langevin calculations is that it may converge to the wrong result and/or lead to runaway solutions [14][15][16]. The problem of runaway solutions was cleverly taken care, in the one-dimensional non-relativistic case, by modifying the action with an extra term whose strength is taken to zero at the end [17,18]. It would be important to verify and validate this method by independent means.
In the present work we show that, the sign problem is very mild in the parameter region investigated using other methods [10,11] and direct investigations of the polarized systems can be performed. We study the simple case of fermions interacting through an attractive s-wave interaction in 1 + 1 dimensions, with an eye toward extending our work to atomic and nuclear systems 2+1 and 3+1 dimensions. Preliminary analysis indicate that the sign problem is also very mild in higher spatial dimensions as well.
This paper is organized as follows: In Section II we discuss the model and its lattice discretization. In Section III we describe the details of our Monte Carlo simulations. In Section IV we present the results of our simulations. Lastly, in Section V we summarize our findings and discusses future prospects.

II. THE MODEL
We study the thermodynamics of a dilute system of non-relativistic spin-1/2 fermions interacting via a contact interaction. The Hamiltonian of the system is where σ runs over the fermion polarizations and n σ (x) = ψ † σ (x)ψ σ (x) is the fermion density. We use nonrelativistic natural units = M = k B = 1 with M being the mass of the fermion. We take the interaction to be attractive so G > 0. In one spatial dimension, the s-wave coupling is cutoff independent since fermion-fermion loops are ultraviolet finite. G is related to the binding energy of the unique two-particle bound state by E = −G 2 /4 (in the continuum, infinite volume limit).
The grand canonical partition function for the system can be written using the path integral formulation: with N σ = dx n σ (x). In the path integral formulation the integrand involves the Euclidean action 3) The chemical potentials are defined as µ ↑ = µ + h, µ ↓ = µ − h. µ controls the density of particles while the spin chemical potential h controls the spin polarization. Note that in the path integral the fields ψ are Grassman fields, but we use the same notation as the operator fields appearing in the Hamiltonian to keep the notation simple. For numerical simulations the quartic fermion term associated with the contact interaction is transformed into a quadratic term via a Hubbard-Stratonovich transformation [19]. The discretized action we start with is where A xt is a real, bosonic auxiliary field andμ σ =μ ±ĥ. All hatted quantities denote lattice quantities. The discretized action converges to the right continuum limit when (2.6) The matching conditions above can be derived by integrating out the auxiliary fields using the relation and then take the continuum limit ∆t, ∆x → 0. The continuum time and space limits can be taken separately.
If we take the continuum time limit first we get the partition function for a discretized fermionic system with on-site interactions, an attractive version of the Hubbard model, with Hamiltonian Note that aboveψ's are operators normalized to satisfy the canonical commutation relations {ψ † σ;x ,ψ σ ;x } = δ σσ δ xx . The discretized partition function converges linearly, as O(∆t), to the time continuum limit. We will show this explicitly and discuss how to improve the ∆ → 0 convergence later. When taking the spatial continuum limit, the discretized Hamiltonian above contain errors of order ∆x 2 . In any case, the convergence rate does not change if we replace the matching conditions above with their first order approximation in the ∆t → 0 limit:γ We will use these matching conditions for our simulations.

III. NUMERICAL METHOD
The fermionic part of the action is quadratic in the fields and we can compute analytically the integral over the fermionic variables. The partition function resulting is then Above, A t indicates the slice of the lattice field A at time t. Note that the dependence on spin appears only as a factor multiplying the temporal hopping matrix C. Since the determinant is diagonal in spin space positive quantity since D σ 's are real matrices. For the spin polarized case,ĥ = 0, the determinant is no longer always positive raising the possibility of a sign problem. The sampling will then be done using the positive probability distribution The sign of the determinant is then included in the measured observable Above · pq indicates averages with respect to the phase quenched probability distribution p pq (A). An important observation in this paper is that for all simulations described here we did not encounter a negative determinant. This is not to say that there are no A configurations that lead to negative determinants, but that for the parameter region used in this study the probability of such configurations is extremely small. By for the most expensive step of the calculation is the computation of the fermion determinant, so we have optimized its computation. The fermion matrix, when the entries corresponding to a time-slice are grouped together, has the following structure Using the sparsity of the temporal blocks, the calculation of the determinant det D σ can be reduced to a smaller problem using Schur's complement methods, similar to the procedure used for Wilson fermions in lattice QCD [20] (note that a similar expression was derived earlier for the Hubbard model using operator methods [21,22]). We have The cost of calculating the fermion determinant is reduced from a complexity of O( . Further simplifications appear because the matrices C(A) are diagonal, so their multiplications has only a cost linear in N x and because matrix B does not depend on the fields A and its determinant and inverse are calculated only once. For periodic boundary conditions the eigenvectors of B are plane waves (ψ k ) x = exp(ikx) with k an integral multiple of 2π/N x . The inverse is then As it turns out the partition function in the reduced form as it appears in Eq. 3.6 is very similar to the one for the Hubbard model [22]: the fermionic contribution in given by a similar product of matrices with a diagonal matrix including the auxiliary field contribution and a non-diagonal matrix, similar to B, that encodes spatial hopping. The main difference is that the off-diagonal matrix is exp(∆tK), where K xx ∝ δ x,x +1 + δ x,x −1 is the spatial hopping matrix, whereas for us B ∝ 1 + ∆tK. The Trotter time discretization used for the Hubbard model has errors of the order O(∆t 2 ) [22,23], whereas our discretization has errors O(∆t) (this can be seen by considering that the two discretizations differ at order ∆t 2 for one step, but the full expression requires β/∆t steps.) A simple way to improve our calculation is to replace the B matrix above with the exponentiated expression. This amounts simply to replacing the eigenvalues in the precalculation of the inverse of the B matrix with λ k = exp[2γ sin 2 (k/2)] rather than 1 + 2γ sin 2 (k/2). In Fig. 1 we show the continuum time limit using the two versions of the B matrix, for two observables used in this study. We see that they both converge to the same limit, but the exponentiated version of the B matrix converges faster, as expected.
In a similar vein, we can improve the behavior of the space continuum limit. Note that using the exponentiated form the eigenvalues of the B matrix are λ k = exp[2γ sin 2 (k/2)] = exp[∆tÊ(k)] withÊ(k) the lattice dispersion relations which for small k approximate the continuum onesÊ(k) ≈ E(k) = (k/∆x) 2 /2. A simple improvement is to replace the dispersion relations in the eigenvalues with the continuum form, that is use λ k = exp[γk 2 /2] with k ∈ {− N x /2 , . . . , −1, 0, 1, . . . , N x /2 } in units of 2π/N x . In the left panel of Fig. 2 we show the continuum limit using the lattice and continuum dispersion relations. We can see that both converge to the same limit but the convergence is much faster for the calculation using the continuum dispersion relation. In the right panel of Fig. 2 we look at the thermodynamic limit. For this calculation we use the exponentiated form for the B matrix and the continuum dispersion relations. We set ∆t = 0.025 and ∆x = 1.0 which produce values indistinguishable from the continuum limit at the level of stochastic error bars. We see that for this rather high density, the thermodynamic limit can be achieved with L ≥ 20.
Our simulations were carried out using two sampling methods: a straightforward Metropolis sampling and Hybrid Molecular Dynamics (HMC). In the Metropolis method we generated new fields by proposing global random changes in the auxiliary field A, with the shift at every point limited by a step-size parameter adjusted to produce an acceptance rate around 50%. The probability distribution used in the accept/reject step is the one given in Eq. 3.3, with det D as given by Eq. 3.6 but with the matrix B replaced withB associated with the continuum dispersion relations For the HMC we have to evolve the field A and its conjugate momentum π according to the canonical equations of motion induced by the Hamiltonian H = xt π 2 xt /2 + S(A). The force term is derived from the integrand in Eq. 3.3π where U t (A) =B −1 C(A t+Nt−1 )...B −1 C(A t ). The HMC sampling requires the determinant to be positive, but this seems to be the case in all our simulations; we will discuss this point below in more detail. We carried out simulations for many different parameters using both Metropolis and HMC sampling and the results agreed in all cases.
A point to note is that in the evaluation of matrix U t (A) for large N t , required for both Metropolis and HMC sampling, the product matrix U exhibits numerical instabilities. This problem is encountered in Hubbard model simulations and it was found that the instability can be overcome by using a factorization of the intermediate results [22]: the product is split into subproducts of matrices that can be computed directly intermixed with factorizations. A similar instability appears when computing the force term for the HMC algorithm: A simple optimization is to compute [1 + exp(N tμσ )U 0 (A)] −1 and then evaluate the other time-slices using the property (3.10) The matrices C(A t ) are diagonal and trivial to invert and theB matrix and its inverse are precomputed. This iteration can be used safely for a small number of time-slices, but for large N t we need to recompute U t (A) from time to time. In our runs we found that subproducts of up to 20 matrices can be performed safely when using double precision arithmetic. One final issue to discuss for our numerical simulations is the possibility of trapping in the positive determinant region of the configuration space. This issue arises even for non-polarized simulations where det D = (det D ↑ ) 2 and the determinant is positive. While the determinant is positive, the individual spin determinant det D ↑ changes sign and the positive and negative regions are separated by borders where the determinant is zero. Since both Metropolis and HMC sampling rely on small (or smooth) changes in the configuration, the edges of the same sign regions act as potential barriers and it is possible that the simulation becomes trapped. This lack of ergodicity was studied in the context of Hubbard model for both Metropolis sampling [24] and HMC [25]. For polarized systems the same-sign regions of det D ↑ and det D ↓ are not identical: as the polarization increases the borders between these regions are broadened and the determinant det D has fluctuating signs. Since we do not see any sign fluctuation, we believe that the probability distribution we sample has small (almost vanishing) overlap with these regions.
As a further test, we changed the integration manifold to avoid trapping, a method inspired by our work on thimbles [26]. The basic idea is that a generalization of the Cauchy theorem guarantees that the path integral remains unchanged for rather large set of deformations of the integration domain in the complex field space. We consider here a shift of the integration domain in the imaginary direction, which amount to adding a constant imaginary component to each field variable: A xt → A xt + iδ. While the integral remains the same, the integrand becomes complex and exhibits phase fluctuations. For small shifts, the fluctuations will be mild except close to the regions where the determinant changes sign. As we move away from the real integration domain the potential barriers associated with the zero determinant are lowered, removing trapping, but the sign fluctuations become larger. If the results are biased by trapping we should see a change in the value of the observables as we vary the imaginary shift δ away from zero. In Fig. 3 we plot the density and the average phase as a function of δ for one ensemble. The parameters used for this ensemble are G = 1/ √ 8, L = 61 as for the other tests, but we move to lower temperature β = 8 and higher densities and polarizations  [10] by an analytic continuation from imaginary spin chemical potential. This method is expected to work well for βh ≤ π, the region indicated by solid lines, but it quickly becomes unreliable as we move to higher polarizations, as we can see by comparing our results with the dotted lines. βµ = 4 and βh = 4 as trapping is expected to appear more readily at low temperature and high densities. The lattice parameters are set to ∆x = 1 and ∆t = 1/20. We see that the density remains unchanged as we vary δ while the sign fluctuations become significant. The spin density is also unaffected by the shift. We conclude that it is unlikely that our simulations are trapped.

IV. RESULTS
The observables that we focus on here are the number and spin density: where the average is taken over Monte Carlo configurations. For our simulations we set G = 1/ √ 8 and β = 8. The results presented are computed using ∆x = 1 and ∆t = 1/8 which are indistinguishable from the continuum results at the level of the error-bars (see Fig. 1 and Fig. 2). The volume is set to L = 61 which is close to the thermodynamic limit. We set the parameters to these values to compare our results with those from imaginary polarization studies [10] and complex Langevin [17]. All of our results are computed using 500 statistically independent measurements.
Our main result is presented in Fig. 4: for three values of βµ at −2.0, 0.0, 2.0 we sweep the spin chemical potential between 0 ≤ βh ≤ 4.5. In addition to our data, we also include in Fig. 4 the result of analytically continuing the results of calculations performed at imaginary chemical potential, a study carried out in [10]. One finds at small polarization reasonable agreement between our method and analytic continuation, however at large polarization there is clear disagreement. The imaginary polarization method is expected to fail for large polarizations, when βh π, but we see large deviations for polarizations as small as βh 2. This may invalidate the expectations that interesting features of the phase diagram lie in the region where analitycal continuation is reliable [9].
In Fig. 5 we compared our results with analytical continuation and the complex Langevin methods [17]. We use a polarization of βh = 2 and focus on the low chemical potential region where there seems to be a tension between the analytical continuation and the complex Langevin results. In this region our results seem to agree at low density with the imaginary-polarization method but as we increase the density the results seem to agree better with complex Langevin results. Since we do not have error estimates for the other methods, we cannot make a sharp statement. Our results for the density compared with results from analytical continuation from imaginary polarization [10] and complex Langevin simulations [17]. The parameters for these simulations are G = 1/ √ 8, β = 8, and βh = 2.0. The right panel indicates the difference between our results and the other methods. The error-bars indicate the statistical uncertainties of our results since we did not have information about the error bars on the results from the other approaches.

V. CONCLUSIONS
We study using Monte Carlo methods the physics of a system composed of two species of fermions in one spatial dimension, including unbalanced systems with unequal densities of the two species where a sign problem exists. We find that the sign problem is extremely mild for a vast range of parameters corresponding to densities small compared to the lattice cutoff scale which is the relevant parameters region for the continuum limit. This is in contrast to the parameter values corresponding to densities close to half-filling (but not exactly equal to it), relevant to the physics of the Hubbard model, where a severe sign problem occurs. The simulations in this work were carried out at G = 1/ √ 8, but we also find no negative determinants in simulations with couplings as high as G = 4/ √ 8 and it's possible that even higher couplings have mild sign problems.
The range of density and spin studied in this work is wide. For some parameter values we are able to compare our results to previous calculations using the complex Langevin or the imaginary chemical potential methods. Detailed comparisons are not possible as the published results lack proper error bars but the general lesson is that out calculations are very close to the complex Langevin results (lending credence to its convergence to the right result) and agree with the imaginary chemical potential results at small enough asymmetries. For larger asymmetries our results are very different from the imaginary chemical potential ones. This is hardly surprising giving the nature of the analytic continuation from imaginary to real chemical potentials.
One potential problem with our calculations is the possibility that the Monte Carlo chain is trapped in a region of positive determinants. Besides producing the wrong result this would give the false impression that the sign problem is milder that it is in reality. Many pieces of evidence were offered to argue against that. The first is that calculations done with the Hybrid Monte Carlo method (more prone to trapping) agrees with the one done with the Metropolis algorithm (which is not likely to be trapped). We also used ideas from our work on the "thimble" approach to perform a calculation of the path integral not over the real variables, but over a hyperplane slightly shifted on the imaginary direction. Since now the determinants are complex the regions with opposite signs of the determinant are not separated by a repulsive barrier where the determinant vanishes. We observed no difference in the results. Finally, the agreement with other methods for parameter values where these methods are more reliable lend further confidence in our calculation.
It is unclear at this point whether similar methods can be easily extended to higher dimensional systems. Work in this direction is being pursued and results will appear elsewhere.