Particle Projection Using a Complex Langevin Method

Using complex stochastic quantization, we implement a particle-number projection technique on the partition function of spin-1/2 fermions at finite temperature on the lattice. We discuss the method, its application towards obtaining the thermal properties of finite Fermi systems in three spatial dimensions, and results for the first five virial coefficients of one-dimensional, attractively interacting fermions.


Introduction
Calculating the thermodynamics of any quantum many-body system is a challenging but generally interesting problem for a variety of applications. By far the main roadblock in Monte Carlo calculations of such systems is the sign problem, whose impact affects fields ranging from ultracold atoms [1,2] to QCD [3]. In the last few years, a concerted effort to tackle this issue has been taking place, as several techniques have been put forward to that effect (see e.g. [4]). In this work, we make use of one of those ideas, namely that of complex stochastic quantization (i.e. we use the complex Langevin approach) to explore the feasibility of particle-number projection in the calculation of the thermal properties of finite clusters and in the numerical determination of virial coefficients.
We begin by recalling that the equilibrium thermodynamics is encoded in the grand canonical partition function, whereĤ is the Hamiltonian,N is the particle number operator, β is the inverse temperature, and µ is the chemical potential. This in itself is nearly impossible to solve for essentially all interacting cases.
To gain a better understanding, in particular in dilute regimes, we may expand about z = e βµ = 0, where z is the fugacity. Such a low-fugacity expansion is what we call the virial expansion. The corresponding coefficients accompanying the powers of z in the expansion of the grand-canonical partition function Z are the N-particle partition functions; specifically, for spin-1/2 fermions, where z ↑ = e βµ ↑ and z ↓ = e βµ ↓ are the fugacity for the spin-up and spin-down particles, respectively, and Q n ↑ ,n ↓ are the canonical partition functions for n ↑ + n ↓ particles.
To further analyze the thermodynamics of many-body systems we can study the grand-canonical potential, −βΩ = ln Z. Again, expanding in powers of z ↑ and z ↓ we achieve an expansion about the dilute limit. This yields, Where b n ↑ ,n ↓ are the "virial coefficients" of order N = n ↑ + n ↓ . These are of particular interest as they give access to the equation of state (EoS), and therefore to observables for these many-body systems (i.e. density, pressure, compressibility, etc.). Moreover, current experiments with ultracold atoms have gained increasing access to finely controllable parameters such as temperature, coupling, polarization, mass imbalance, and trap shape (harmonic, hard wall, lattices), which allows direct comparison of thermodynamic predictions in the dilute limit, which is the focus of this work.
In this contribution, we show our progress towards accessing Q n ↑ ,n ↓ for the universal regime called the unitary Fermi gas [5], with the aim of characterizing properties of finite systems, i.e. those derived from the Helmholtz free energy −βF = ln Q n ↑ ,n ↓ of finite clusters. Our second main goal is to use a similar approach to compute virial coefficients, b n ↑ ,n ↓ , at various orders for the one-dimensional Fermi gas with attractive, short-range interactions. To that end, we use a variant of the complex Langevin method as a means of stochastic quantization.

Preliminaries
Below, we will make use of the path integral representation of the grand-canonical partition function, namely where σ is a Hubbard-Stratonovich field representing a two-body interaction and the system is assumed to be non-relativistic. The specific dynamics of the system at hand is encoded in the matrix M[σ]. In particular for a two-species non-relativistic system, it is not difficult to show [6] that where U[σ] does not depend on the chemical potential. The above equation is the source of the so-called sign problem in polarized non-relativistic matter, as U[σ] is real for attractive interactions but z ↑ z ↓ leads to different determinants for each species (generally of different, fluctuating sign).
In the sections below, we implement a method to compute canonical partition functions and virial coefficients, both based on the idea of particle projection, which is in essence a Fourier transform on a variable φ if the fugacity is replaced by taking z → ze iφ . Such a transformation, which makes z a complex variable, introduces a phase problem in conventional lattice Monte Carlo approaches.

Complex Langevin and Stochastic Quantization
The idea of stochastic quantization is that a random process that obeys the Langevin equation will yield configurations σ that are distributed according to e −S [σ] . Here, θ is Langevin time and η is gaussian random noise (see e.g. [3]). We assume our Hubbard-Stratonovich field σ is real. However, when S [σ] is complex, the Langevin equation naturally implies that σ must be complexified into σ = σ R + iσ I . Thus, stochastic quantization requires complexified fields and the (spacetime local) Langevin equation now has both real and imaginary parts, where the Langevin time has been discretized as θ = n , and is the time step. Again, η, is the real Gaussian noise such that η(n) = 0. The drift terms are now specified as,

Particle Projection Method
As a starting point, we use the auxiliary-field representation of the grand-canonical partition function on the lattice, which reads where σ is the Hubbard-Stratonovich auxiliary field, and U[σ] encodes all the properties of the Hamiltonian of the system. Expanding in terms of the low-fugacity limit we rewrite Particle projection consists in implementing delta functions that select n ↑ and n ↓ in the above sum, that is δ n ↑ ,n and δ n ↓ ,m . To that end, we first take , such that our fugacity expansion of Z becomes, We may thus project out the desired individual partition functions Q n ↑ ,n ↓ using a Fourier integral in the usual way, namely In practice, we may define Φ = (σ, φ ↑ , φ ↓ ) as a new auxiliary field variable and perform the Fourier integral as part of stochastic process (see below). For this purpose, we include φ ↑ and φ ↓ as part of our action and integral as follows, where From this point on, thermodynamic properties of finite clusters can be obtained in the usual way by formally taking derivatives of Q n,m with respect to the desired source and computing expectation values in the extended ensemble defined by the variable Φ.

Virial Projection
By definition, the virial expansion of the partition function Z can be written as where z is the fugacity, b n are the virial coefficients and Q 1 is the single-particle partition function. Based on the above expression, we may use Fourier projection to obtain the virial coefficients: More generally, we may define the function where the last identity is easy to see from our first definition above. This last equation, however, does not help in computing any of the above: while we have a pathintegral form for Z, the presence of the logarithm makes a direct evaluation impractical. Nevertheless, inserting the path-integral form of the partition function and differentiating with respect to ξ we obtain a more useful expression: where P[σ, z] ≡ det M[σ, z]/Z[z]. Using shorthand notation for the expectation value with respect to the above weight P, In practice, we do not need a continuous representation for the Fourier integral, particularly because we cannot include this as part of the stochastic process as in our particle projection method due to the logarithm between these integrals (i.e. the integral over the projecting angles does not appear in the normalization of P). Rather, we may use an exact discrete representation via the associated angle values φ k = 2πk/N k , where k = 0, . . . , N k − 1 and N k is the number of discretization points, i.e.
What is interesting about Eq. (23) is that, by using the complex Langevin method to determine the complex expectation value as a function of φ k and ξ, followed by summing over φ k with appropriate phases, we may have access to b n itself.
In this approach, varying ξ could be very advantageous, as we know the exact ξ dependence of the final result via the last identity of Eq. (23). To enhance the sensitivity of the approach, it is preferable to take ξ to be as large as possible, since the b n tend to be small numbers. However, we anticipate that taking ξ ∼ 0.1 or larger could introduce considerable noise in the calculation. Thus, a priori (and ideally), we expect to find a window of values of ξ makes the approach viable. Our knowledge of the functional dependence mentioned above should connect the data points for all ξ with a single curve. In fact, the expected dependence can be accounted for by plotting as a function of ξ and fitting a constant. If the number of Fourier discretization points N k is large enough, we expect to see such a constant behavior beyond the statistical precision. We pursue such an approach below.

Particle Projection Results
Figure 1: Running average over Langevin time, τ, where τ 0 is the total number of Langevin time steps, of the total particle numbers, n ↑ and n ↓ . Separate runs are seen to converge to the expected particle number for the 2 + 1 case (blue) and the 6 + 3 case (red). The data corresponds to V = N 3 x , N x = 6, β = 3.0, tuned to the unitary limit.
The simplest and most important verification of this method is computing the projected particle numbers for each flavor. The observable in the grand canonical ensemble is given by which is to be averaged over Langevin time. As Langevin time progresses, the running average for the particle numbers of each flavor, n ↑ and n ↓ , should converge to the selected values. This is demonstrated in Figure. (1), where we indeed verify that the method converges to the correct solutions. Figure. (1) is computed in a relatively small spatial volume, V = 6 3 , but at β = 3.0, which implies N τ = 60 lattice sites in the time direction. Computations for larger volumes (i.e. N x = 8, 10, 12) are needed to limit finite-size effects. Likewise, larger values of β are needed to extrapolate to the continuum limit (see below).

Virial Projection Results
As mentioned above, the goal of this method is to access the virial coefficients b n by carrying out calculations for varying values of ξ. A first step toward extracting b n is to show that the results over the range of ξ yield a constant value, at least in some window. Fig. 2, shows results for b n (ξ) at β = 8.0 and 10.0, both at √ βg = 0.1. We find consistent results for up to order n = 5 for ξ = 0.05−0.1, but also find statistical noise at smaller ξ, which is particularly evident in the higher order virial coefficients; the latter correspond to progressively higher frequencies in the Fourier projection. Nevertheless, the roughly constant behavior allows us to extract estimates of b n up to n = 5. Systematic effects in the number of Fourier points N k and the Langevin time T are currently under study. The data shown here corresponds to N k = 100 and T = 20. Computing the above results for varying values of the coupling, λ ∈ [0.1, 1.0], and averaging the results over the consistent values (excluding outliers, which are dominated by noise), we obtain a first estimate of b n vs. λ. An extrapolation of our answers to λ = 0 allows us to compare with the exact solutions for the non-interacting case. Extrapolating in the other direction, i.e. to larger λ it may be possible to obtain approximate virial coefficients at strong coupling. These results for b n vs. λ at β = 8.0 and β = 10.0 are shown in Fig. 3, together with the exact solution for b 2 (λ) (see Ref. [7]).
Notice in Fig. (3) at low coupling our results match extremely well the analytical solution for b 2 (λ), but that agreement deteriorates as λ increases. More statistics may be needed at larger coupling to obtain more consistent results. Fig. (3) shows continuous results for most orders of b n , but again with inconsistencies at higher order as λ increases. Both are encouraging results that point to the possibility of computing high-order virial coefficients and extrapolating to larger coupling with fair accuracy. This approach, in that sense, complements the conventional, high-precision way to obtain virial coefficients, based on solving the n-body problem and computing canonical partition functions.

Summary and conclusions
In this contribution we have presented two calculations that use the complex Langevin approach to extract information about finite systems from finite-temperature, grand-canonical calculations. We focused on particle-number projection to obtain thermodynamics via canonical partition functions, and to obtain virial coefficients. We have shown that both methods converge to reasonable results within some constraints on various parameters. More data are needed in order to estimate the free energy for various systems and particle numbers, and to further understand systematic effects. To extrapolate to the continuum limit, we must take 1 λ T N x (where λ T = 2πβ is the thermal wavelength), which amplifies the computation time for these calculations as both β and N x must be taken to be large. Due to this limitation, we are exploring ways to optimize these calculations. We have shown, however, that particle projection using complex Langevin does converge to the correct state, and with the correct number of particles for each flavor.
Furthermore, our method for virial projection has yielded promising results for one-dimensional systems, demonstrating possibilities for computing higher order coefficients (i.e. b 3 , b 4 , and above) as well as extrapolation to stronger values of the coupling. We are currently exploring finite-size effects of this method by varying the inverse temperature, β, and the volume. In addition, we are investigating the effects of varying the total number of Fourier points N k , as decreasing this number likewise decreases the total run time of the simulation (linearly). Increased statistics (longer total Langevin time) is needed for improved accuracy for larger virial coefficients, and particularly for larger values of the coupling.
Both particle projection methods have produced promising results that, with more in-depth study, may provide access to various observables in finite systems that remain unknown and may be directly compared with ultracold atom experiments. Looking beyond the latter, the ultimate goal is to perform full-fledged, finite-temperature calculations of finite nuclei with realistic interactions, for which the present work provides a proof of principle.
This material is based upon work supported by the National Science Foundation under Grant No. PHY1452635 (Computational Physics Program).