A lattice pairing-field approach to ultracold Fermi gases

We develop a pairing-field formalism for ab initio studies of non-relativistic two-component fermions on a $(d\!+\!1)$-dimensional spacetime lattice. More specifically, we focus on theories where the interaction between the two components can be described by the exchange of a corresponding pairing field. The introduction of a pairing field may indeed be convenient for studies of, e.g., the finite-temperature phase structure and critical behavior of, e.g., ultracold atomic Fermi gases. Moreover, such a formalism allows to directly compute the momentum and frequency dependence of the pair propagator, from which the pair-correlation function can be extracted. For a first illustration of the application of our formalism, we compute the density equation of state and the superfluid order parameter for a gas of unpolarized fermions in $(0\!+\!1)$ dimensions by employing the complex Langevin approach to surmount the sign problem.


Introduction
Pair formation of fermions governs the dynamics in a tremendous variety of systems at vastly different length scales, ranging from electrons in solids and ultracold quantum gases to stronginteraction matter and the physics of neutron stars.In gases of interacting spin-1/2 fermions, for example, pairing of spin-up and spin-down particles at diametrically opposite points of the Fermi surface plays a very prominent role, as it leads to the opening of a gap in the quasiparticle spectrum accompanied by the formation of a superfluid state at sufficiently low temperatures [1,2].Interestingly, such a formation of pairs alongside the emergence of superfluidity even appears to hold for strongly correlated systems, as has been found in detailed theoretical and experimental studies of ultracold Fermi gases (see Refs.  and Refs.[25][26][27][28][29][30][31][32][33][34][35][36], respectively).Naturally, an understanding of the microscopic properties of fermion pairs is fundamentally relevant to gain an insight into the macroscopic properties of such fermion gases.
In the present work, for concreteness, we restrict ourselves to systems which can be described by a Hamilton operator of the form Ĥ = d d r − ψ † σ (r) where we have set ħ h = 1 and d determines the number of spatial dimensions.Unless stated otherwise, we always assume summation over repeated indices associated with the species σ =↑, ↓.
The field operators ψσ and ψ † σ associated with fermions of spin σ are assumed to obey the usual anti-commutation relations.We moreover assume that the interaction, parametrized by the coupling g > 0, is attractive.From a phenomenological standpoint, we mainly have the development of an alternative lattice formalism for the description of ultracold atomic Fermi gases in mind (which should be extendable to nuclear matter).Still, we would like to add that suitably discretized versions of the model described by Eq. ( 1) can be directly related to condensed-matter systems, most prominently to the Hubbard model, see, e.g., Refs.[37,38] for recent discussions.
The macroscopic properties of the many-body system described by the Hamilton operator (1) can be extracted from the grand-canonical partition function: where k B = 1 and β = 1/T is the inverse temperature.The particle number operators N σ = d d r ψ † σ (r) ψσ (r) couple to the chemical potentials µ σ associated with the two spin projections.
The partition function can in principle be computed in various ways.In the present work, we restrict ourselves to a path-integral formulation of the partition function.Unfortunately, for an operator of the form (1), exact results are not available for general spatial dimension d, temperature T , and chemical potentials µ σ , which makes the use of numerical approaches necessary.Instead of computing the path integral directly, however, stochastic path-integral approaches generally require to remove the Grassmann-valued fermion fields by introducing a suitably chosen auxiliary field.This basically corresponds to performing a Hubbard-Stratonovich transformation [39,40].This step is by no means unique.On the contrary, a variety of auxiliary fields have already been used in the literature, see, e.g., Refs.[41][42][43] for reviews.
Often, the physics of interest suggests the use a specific type of auxiliary field.For example, auxiliary fields closely related to the density are frequently used for computations of the density equation of state and related quantities, see Ref. [43] for a review.However, studies of the propagation of fermion pairs requires the construction of suitable (and in general nontrivial) operators.In such a situation, it may therefore be advantageous to employ an auxiliary field which can be identified with the field associated with such fermion pairs.In particular, for studies of spontaneous symmetry breaking associated with a superfluid state and critical behavior close to the finite-temperature phase boundary, such a pairing-field formulation may be appealing.In the continuum limit, such a formulation indeed exists for studies of, e.g., ultracold Fermi gases, which we shall briefly summarize in Sec. 2, see, e.g., Refs.[44][45][46][47][48][49] for reviews.This formalism turns out to be particularly convenient for studies of the phase diagram of ultracold Fermi gases, see, e.g., Refs.[6,[17][18][19][20]50].In Sec. 3, we develop a lattice formulation of this pairing-field formalism.The application of this formulation is then demonstrated in Sec. 4 with the aid of an exactly solvable zero-dimensional fermion model, i.e., the zero-dimensional version of the model described by the Hamilton operator of Eq. ( 1).We add that our lattice formulation of the well-known continuum formulation of the pairing-field formalism has a sign problem, even for spin-balanced systems, see also Refs.[51,52] for a discussion.This is not unexpected, as the introduction of the pairing field transforms the original purely fermionic model [6,[17][18][19][20]50] into a complex scalar field theory, see our discussion in Secs. 2 and 3.Such theories are known to have a sign problem on a spacetime lattice.In order to deal with this sign problem, we employ the complex Langevin (CL) approach, which has been employed before to study relativistic as well as non-relativistic complex scalar field theories [53][54][55].Moreover, the CL approach has been successfully applied to compute properties of ultracold Fermi gases described by the Hamilton operator (1), see Refs.[56][57][58][59][60][61][62].For reviews on the CL approach [63], we refer the reader to, e.g., Refs.[43,[64][65][66][67][68].In Sec. 4, we discuss this in more detail and also present results for the density equation of state as well as the expectation value of the pairing field.Our conclusions are presented in Sec. 5.

Pairing field formalism in the continuum limit
For reference, we briefly summarize the derivation of the pairing-field formalism in the continuum limit in this section.It is based on a specific Hubbard-Stratonovich transformation and is well-known in the literature.More detailed discussions can be found in, e.g., Refs.[44][45][46][47][48][49].
The starting point of our discussion is the path-integral representation of the partition function (2): The (fermionic) action S F reads where the fields ψ ↑ and ψ ↓ are associated with spin-up and spin-down fermions, respectively.We can now bosonize the theory by performing a Hubbard-Stratonovich transformation [39,40].To be more specific, we insert a suitably chosen constant into the path integral.
Here, we choose with a complex-valued auxiliary bosonic field φ = φ(τ, r).We refer to the field φ as the pairing field and assume that it carries the same quantum numbers as the fermion composite ψ ↑ ψ ↓ .Shifting this field and its complex conjugate according to we arrive at the paritially bosonized action: From this action we deduce that the complex scalar fields mediate the interaction between the fermions.A resonance in this channel indicates the formation of pairs of spin-up and spin-down fermions which may condense at sufficiently low temperatures.Note also that i.e., the expectation value of the auxiliary field φ is identical to the expectation value of ψ ↓ ψ ↑ associated with a pair of fermions, justifying the name pairing field.This expectation value is of particular interest as it represents an order parameter for U(1) symmetry breaking as associated with the formation of a superfluid state.Indeed, from Eq. ( 7), we deduce that a finite expectation value 〈φ〉 introduces a gap in the fermion spectrum, as it is the case in standard Bardeen-Cooper-Schrieffer (BCS) theory [1,2].Since the action ( 7) is only bilinear in the fermion fields, we can integrate them out to obtain the following path integral: where we assume that irrelevant normalization factors have been absorbed into the path integral measure.Hence, the bosonic action S B associated with this path integral reads The fermion matrix M [φ, φ * ] appearing in the functional determinant is given by Since the pairing field also appears within this matrix, we are left with a nonlocal bosonic theory.Properties of the pairs of spin-up and spin-down fermions can be studied by computing correlation functions of the pairing field, such the two-point function 〈φ(τ, r)φ * (τ ′ , r ′ )〉.
We close this section by noting that, at first glance, this fully bosonized theory appears well suited as a starting point for a stochastic evaluation of the path integral.Evaluated on a constant auxiliary field, the action S B yields nothing but the effective potential in the meanfield approximation, which is indeed real-valued and bounded from below.As we show below, however, the action becomes complex on a spacetime lattice when evaluated on general field configurations, i.e., with a nontrivial dependence on τ and r.
3 Pairing-field formalism on the lattice

Preliminaries
To obtain a lattice formulation of the pairing-field formalism, we begin by replacing the fermionic field operators ψσ (r) in Eq. ( 1) by corresponding lattice field operators ψσ,r i which are restricted to a d-dimensional (hyper-)cubic lattice of side length L and spacing a x with periodic boundary conditions.The number of lattice sites in each spatial dimension is given by N x = L/a x .The Hamilton operator (1) then reads Here, r i determines the position on the spatial lattice and the matrix D ′ ∆ represents a finite difference approximation of the continuum Laplace operator.One possible realization of D ′ ∆ in one spatial dimension with r j = ja x êx is given by the central second-order difference: Here and in the following, we only present non-zero entries of matrices.Note that the entries in the top-right corner and the bottom-left corner result from the implementation of periodic boundary conditions.With our concrete choice for the discretization of the Laplace operator specified on the right-hand side of Eq. ( 13), the sums over spatial lattice sites in Eq. ( 12) effectively reduce to a sum over all points and their nearest neighbors.However, because said choice for the discretization of the Laplace operator is not unique and may even be optimized such that the convergence towards the continuum limit is improved, we opt to keep the sums in Eq. (12) in their general form.At this point, we would also like to mention that, for the specific discretization (13) of the Laplace operator, the lattice Hamiltonian is identical to that of a Hubbard model, see, e.g., Refs.[37,38].However, exploring this aspect is beyond the scope of the present work.
Following the standard procedure for the derivation of a path-integral expression for the partition function by introducing coherent basis states and a discretization of the imaginary time interval [0, β) into N τ slices of length a τ = β/N τ (see, e.g., Ref. [69]), we arrive at the following lattice action: with the calligraphic S distinguishing the lattice action from actions in the continuum.Note that the action S F does not depend explicitly on the lattice scales L, a x , and a τ .These quantities as well as the physical parameters g, β, and µ σ have been absorbed in suitably chosen dimensionless quantities, which also causes the Grassmann-valued field variables ψ σ,τ i ,r j (fields evaluated at time τ i and position r j ) to be dimensionless.To be specific, we have rescaled the chemical potential with the temporal lattice spacing: The coupling g has been rescaled and rendered dimensionless as follows: where r = a τ /a 2 x is the dimensionless lattice spacing ratio and λ = β 1−d/2 g is the dimensionless coupling.Finally, the dimensionless Laplace operator is given by Note that, in our derivation, we are in principle free to choose the specific form of the discretization of the spatial derivatives.However, the discretization of the temporal derivative, specifically the appearance of a backward derivative, follows from the construction of the path integral.In case of relativistic theories, the situation is different as the form of spatial and temporal derivatives is constrained by Lorentz invariance.
We add that the chemical potentials enter our lattice action effectively in the form of constant temporal gauge fields.Indeed, in the continuum limit, a simultaneous transformation of the fermion fields and the chemical potentials exists which leaves the action invariant.This is known as the Silver-Blaze symmetry, see, e.g., Refs.[70][71][72][73][74][75] for detailed discussions.In our lattice theory, this symmetry is broken by the presence of the temporal lattice and would only be recovered in the continuum limit.In order to preserve this symmetry on the lattice and improve the convergence to the continuum limit, we replace (1 + μσ ) in the discrete temporal derivatives by e μσ , as advocated in Ref. [76].In practice, this is relevant to obtain accurate results in a regime of small chemical potentials.

Matrix notation
For our development of a lattice pairing-field formalism below, it is convenient to introduce a specific notation for the fields and operators.In this notation, field configurations are represented by column vectors that contain the field values at all lattice sites in an arbitrary but defined order, i.e., where x } is used to enumerate all spacetime lattice sites within the box.A scalar product of two such configuration vectors is associated with a spacetime integral over the product of the two fields in the continuum: With the aid of suitably defined matrix multiplications, we can now construct the various terms in the lattice action.For the spatial derivative, we define the matrix D ∆ , which extends the purely spatial matrix D′ ∆ : Note that this is not yet the term associated with the spatial derivative which appears in the action (14).In fact, the fermion fields are evaluated at two different points in time in the spatial-derivative term in Eq. ( 14), which is not the case in Eq. (20).To take into account the difference in the lattice sites in time direction, we define a (time) retarder operator R − and a (time) advancer operator A − , such that and This definition implies that A and that the matrices inherit the temporal boundary conditions of the field.The subscript of these operators refer to boundary conditions: antiperiodic (-) or periodic (+) boundary conditions in time direction.In the present work, we only need those associated with antiperiodic boundary conditions.In the special case of a 0+1-dimensional theory (i.e., d = 0), the retarder operator R ± has the following matrix representation:1 The retarder and advancer operators allow us to define discrete temporal derivatives as and It follows that where we have used Eq. ( 23).
For the interaction term, we introduce the vectors (ψ ↑ • ψ ↓ ) and (ψ * ↑ • ψ * ↓ ) with • being the element-wise Hadamard product.With this notation at hand, we can rewrite our discrete action (14) as follows This form of the fermionic action represents the starting point for our development of a lattice pairing-field formalism in the next subsection.

The discretized pairing field
Analogously to the continuum theory, we now perform a Hubbard-Stratonovich transformation to eventually remove the fermionic degrees of freedom.To this end, we again insert a constant factor into the path integral expression of the partition function.We choose where However, in contrast to the continuum, we cannot perform a shift like since due to the time shifting matrix R − in the interaction term in the action in Eq. ( 28) this shift will not create a term that cancels the four-fermion interaction.To mitigate this, one might be tempted to begin the Hubbard-Stratonovich transformation by inserting a different factor of one already containing the time shift, i.e., however, since the time shifting matrices are not positive-definite, such a Hubbard-Stratonovich transformation would not be well defined.Instead, we shift the starred and unstarred bosonic fields independently, i.e., and then obtain the partially bosonized lattice action: As in the continuum limit, the fermions can now be integrated out to obtain the bosonized lattice action: 2 where the fermion matrix M is given by Note that this matrix is in general not positive-definite and therefore the bosonized action is in general complex, which leads to a sign problem in conventional stochastic computations of the path integral, even in the limit of vanishing spin and mass imbalances. 3  2 This can be done by reordering terms of the form ψ † ↓ Oψ ↓ such that they appear as terms of the form ψ in the partially bosonized action which allows to rewrite the path integral as a Gaussian integral in the fermionic degrees of freedom.In the continuum limit, the determination of the new operator O ′ requires an integration by parts.In our lattice formalism, this is done by transposition.For example, we have The fermions can then be conveniently integrated out by eventually introducing Nambu-Gorkov spinors 3 There also exist Hubbard-Stratonovich transformations which do not suffer from a sign problem, as long as masses and chemical potentials of the two fermion species are identical.An example of such a formulation would be the density formulation discussed in, e.g., Ref. [42].This approach comes with a computationally less costly update process in numerical applications, at the cost of a more elaborate calculation of pairing observables, such as the superfluid order parameter and the spacetime-dependent pair propagator, see also Sec. 2. The latter quantities are directly accessible with our pairing-type Hubbard-Stratonovich transformation.As a side effect, it also creates a closer connection to continuum calculations where the use of a pairing field as auxiliary field is more common.
Let us finally add that our partially bosonized lattice action S PB and its continuum analogue in Eq. ( 7) look similar at first glance.However, we would like to emphasize that a naive discretization of the continuum theory does not lead us to our lattice theory.For example, a naive discretization of the continuum theory in Eq. ( 7) would lead to a Yukawa interaction term which differs from the partially bosonized lattice action S PB by This indicates that the difference in the interaction term vanishes in the limit a τ → 0, which also holds for the kinetic terms.In fact, in the continuum limit N τ → ∞ and N x → ∞ (such that N τ a τ = β and N x a x = L remain constant), the finite difference operators become derivatives and, loosely speaking, the effect of the retarder and advancer matrices vanishes.
The discrete theory then indeed becomes the continuum theory given in Eq. ( 7).
Finally, we would like to state again that the specific type of Hubbard-Stratonovich transformation underlying our formalism is not new but has been frequently employed in continuum studies, see, e.g., Refs.[44][45][46][47][48][49].Here, we have derived a lattice formulation of this approach.Let us also add that a large variety of Hubbard-Stratonovich transformations has been discussed in the literature to study models as described by the Hamilton operator (1).In this respect, we emphasize that our present lattice formalism should not be confused with the one considered in Ref. [77], where a density-type field is used as an auxiliary field and homogeneous pairing fields are only introduced via source terms.Moreover, it should be mentioned that, in contrast to the spacetime representation of the auxiliary field used in our present work, lattice formulations built on a purely spatial representation of the problem are also very popular.The latter also include a formulation, where the auxiliary field is introduced by coupling it to the pairing operators [51].In this context, it is worth mentioning that the sign problem has also been studied on more general grounds based on a symmetry classification of theories [78][79][80][81].However, an analysis of whether and how the classifications presented in these works can be applied to our spacetime formulation of theories with an even number of fermion species is beyond the scope of our present study.In any case, our analysis of the fermion matrix (36) indicates that our pairing-field approach comes with a sign problem, in accordance with, e.g., Ref. [52].For recent reviews on the sign problem and the role of Hubbard-Stratonovich transformations in the context of stochastic calculations, we refer the reader to Refs.[42,68,82].

Lattice parameters
Based on the required range and accuracy of the results, we can define criteria that determine the lattice parameters.In the following, we present criteria for the determination of the numbers of lattice sites N τ and N x based on the largest chemical potential βµ max of interest and three empirical constants δ SB , C I and C λ .The latter two control the numerical accuracy of the results.For a given βµ max , the numerically exact solution is approached in the limit C I → ∞ and C λ → ∞, wherein the vanishing deviation caused by the replacement term to restore the Silver-Blaze symmetry is ensured by the increasing temporal lattice size following the limit for C I .

Number of temporal lattice sites
By considering the energy scale set by the interaction and the Silver-Blaze symmetry, we first explore criteria which have to be satisfied by the temporal lattice (in terms of extent and spacing).For our numerical studies presented in Sec. 4, we have in general used the smallest even value of N τ which fulfills all the criteria.However, we have also checked the correctness of our criteria by studying the convergence of the density as a function of N τ for selected parameter sets.
On a spacetime lattice, the reciprocal temporal spacing defines a cutoff for the energies that can be resolved on the lattice.We therefore must ensure that the energy scale E I set by the interaction is sufficiently small compared to this cutoff, i.e., 1/a τ > C I E I with an empirically determined factor C I .Taking into account the dimension of the coupling parameter g, we define the interaction energy scale E I = λa −d x β d/2−1 .For d < 2 this leads us to the criterion In the special case of d = 2, the coupling is dimensionless and does not provide an energy scale.For d > 2, the scale set by the interaction energy imposes an upper bound for the number of temporal lattice sites.
In the derivation of our formalism, we replaced factors of (1 + μ) by factors of e μ to preserve the Silver-Blaze symmetry.Such an improved lattice theory leads to faster convergence towards the continuum theory relative to the unimproved counterpart.Specifically, the difference between improved and unimproved increases quadratically at leading order with increasing values of | μ|.Since we have we can reduce the deviation by simply increasing the temporal lattice extent N τ .In practice, we define a range absolutes of values of interest for the quantity βµ .The upper bound of this range, βµ max , is then used to determine an appropriate value of N τ by solving the equation: This relation sets an empirical upper bound δ SB for the difference.The results in this work are obtained with δ SB = 0.005.However, for the coupling strengths considered in Sec. 4, the temporal lattice spacing is determined by the scale set by the interaction energy rather than the Silver-Blaze correction.

Number of spatial lattice sites
In addition to the temporal lattice size, we have to determine a number of spatial lattice sites.
In this case, we have to ensure that the length of our box L = a x N x is sufficiently large compared to the thermal wavelength λ th,σ = (2πβ/m σ ) 1/2 .To this end, we introduce an empirical constant C λ and require Note that the (inverse) temperature is a phenomenological control parameter rather than an artificial parameter of our lattice theory, but, contrary to that, the spatial extent L of our (hypercubic) lattice is in general an artificial parameter which has been introduced to evaluate the path integral numerically.That being stated, the above inequality leads us to the criterion which determines the value of N x based on the number of temporal lattice sites N τ .In other words, this inequality relates the temporal lattice spacing to the spatial lattice spacing for given values of β and L. We also note, that the required number of spatial lattice sites per dimension scales with the square root of the number of temporal lattice sites, as N x ∝ N 1/2 τ .

Application: Fermions in zero dimensions
To illustrate and discuss the application of our lattice pairing-field formalism, we consider our model in the limit of zero spatial dimensions (d = 0).This is motivated by the fact that exact analytic results in closed form can be obtained for at least some physical observables in this case, providing a clean benchmark for our numerical calculations.

Exact analytic results
The partition function (2) of our model can be computed analytically for d = 0. We find From this formula, we can extract exact results for the spin-up and spin-down densities by taking a derivative with respect to µ ↑ and µ ↓ , such that The total density is defined to be the sum of the spin-up and spin-down density: Below, we often use the densities of the noninteracting system to normalize the densities of the interacting system.We have and n 0 := n 0,↑ +n 0,↓ .We add that the densities of the two fermion species become independent of the chemical potentials in the limit of an infinitely strong attractive interaction: Moreover, for µ ↑ =µ ↓ =0 and g > 0, the normalized densities n σ /n 0 are found to be bounded from below and above: 1 ≤ n σ /n 0,σ ≤ 2. Therefore, this ratio may be considered a macroscopic measure of the effective strength of the interactions in the system.Note that, for infinite repulsion, the densities still depend on the chemical potentials.
The main focus of our illustrative numerical studies will be on the computation of densities.Of course, the actual strength of our formalism lies in the computation of observables which can be constructed directly from the pairing field, such as the pair-correlation function, the quasiparticle gap as well as other order parameters for superfluidity, which can in principle also be computed analytically in d = 0. We will report on such computations elsewhere and only show results for the expectation value of the pairing field in this work.From an analytic calculation of the latter quantity, we obtain 〈φ〉 = 0 since and 〈φ〉 = 〈ψ ↓ ψ ↑ 〉.Note that 〈φ〉 can be used as an order parameter for U(1) symmetry breaking in our model.This allows to relate this quantity to the formation of a superfluid ground state in d ≥ 2 in the long-range limit.
We close this subsection with a comment on mean-field theory.It is also possible to study our zero-dimensional model in the mean-field approximation in the continuum limit.The effective potential U then reads where µ is the average chemical potential of the two species, and the chemical potential asymmetry is measured in terms of the so-called Zeeman field h, The fields φ and φ * represent constant pairing fields.Note that, in U, we dropped terms independent of the field φ and its complex conjugate.
A minimization of the effective potential U yields the ground state φ0 .The pressure equation of state can be extracted from an evaluation of U at the ground state φ0 = 〈ψ ↓ ψ ↑ 〉.From the first derivative of the pressure with respect to the chemical potential µ σ , we then obtain the densities.Notably, for the coupling strengths considered in the present work, the meanfield approximation yields that the density in the interacting system is identical to the one in the noninteracting system, even if our numerical studies (in accordance with the exact results) predict that this is not the case.Even more, for sufficiently attractive couplings (beyond those considered in our numerical studies below), the mean-field approximation indicates that the ground state is nontrivial, i.e., φ0 = 〈ψ ↓ ψ ↑ 〉 ̸ = 0, again in disagreement with the exact solution.Thus, the results from the mean-field approximation are not even correct on a qualitative level.However, this is not unexpected as fluctuation effects (which are missing in the mean-field approximation) are very relevant in low-dimensional systems.

Complex Langevin approach
Since we have integrated out the fermion fields in our formalism, we are left with a purely bosonic field theory.As such, the dynamics of the system is fully encoded in the pairing field φ and its complex conjugate φ * .In order to compute observables with the path integral, numerically it is convenient to generate a finite number suitable pairing-field configurations which allow us to evaluate the path integral in an efficient way by approximating its value with the average of the integrand at the chosen field configurations.Often this is done by using Monte-Carlo (MC) techniques.Unfortunately, standard MC techniques are not applicable in our case since the bosonized lattice action S B is complex, as already indicated above.Therefore, we employ the CL approach [63] to generate field configurations suitable for an efficient computation of observables, see Refs.[43,[64][65][66][67][68] for reviews.
The application of the CL approach requires a complexification of the fields in our bosonized theory.To this end, we first rewrite the pairing field φ and its complex conjugate φ * in terms of their real and imaginary parts which leaves us with two real-valued fields, φ 1 and φ 2 : where φ 1 := Re(φ) and φ 2 := Im(φ).On the lattice, these fields are parametrized by a finite set of complex numbers as described above in our derivation of the lattice action S B : where k = 1, 2. Within our CL approach, these real-valued fields are then complexified, i.e., the field φ 1 and φ 2 are considered to be complex fields.Formally, this is achieved by replacing these fields by a sum of their real and imaginary parts.For the field variables, this replacement reads φ (k) Note that this complexification of the fields implies that we have to deal with four real-valued fields when we compute physical quantities.
In addition to the complexification of the original field variables, the CL approach also introduces a fictitious time t CL .The complexified fields are then evolved along this so-called CL time according to the following set of discretized differential equations for our complexified fields: φ(k) where n is the CL time index of the fields and δt CL is the CL time step, t CL = n δt CL .The noise η is real and Gaussian with 〈η (k) and 〈η (k 1 ) The so-called drift term is defined as This set of coupled differential equations can now be used to generate field configurations for the evaluation of physical observables in the spirit of the original path-integral formulation.
For more details on the CL approach in the context of nonrelativistic many-body physics, we refer the reader to Ref. [43].

Computation of observables
The computation of a physical observable requires the definition of a corresponding operator in our formalism.For example, we can derive an operator O CL N σ for the calculation of the particle number N σ of species σ from the definition of the partition function: where M is the fermion matrix introduced in Eq. (36).From this, we deduce the following "CL representation" of the density operator: Computing the average of this operator with the field configuration samples obtained from the solution of the CL equations, we obtain the particle number of the species σ.To be specific, we have where N CL is the number of CL time steps.Since our conventions for all quantities appearing in the action S B are such that our formalism is completely free of any dimensionful scale, we can strictly speaking only extract dimensionless observables from our CL calculations.Since the density n = N /V depends on the volume of the box, V = (N x a x ) d , it can therefore not be computed directly.Instead, we define a dimensionless density by multiplying the density n with the thermal wavelength volume λ d th : Here, the thermal wavelength λ th = (2πβ) 1/2 is made species-independent by using a mass of one.Note that in 0+1 dimensions the physical density n and the dimensionless density n λ d th are identical as the theory has no spatial extent.Finally, we would like to add that the computation of expectation values of the pairing fields, such as the superfluid order parameter and the pair propagator, can be implemented straightforwardly since the pairing field is the fundamental field in our formalism.For example, we have Expectation values of more than one pairing field can be computed correspondingly.Note that the computation of fermionic correlation functions (e.g., the propagator of one of the fermion species) is more involved as it requires to introduce a source into the path integral.

Results for the density equation of state
For λ = 0, we simply expect to obtain the density distribution of a non-interacting Fermi gas for each of the species.Indeed, as Fig. 1 shows, the simulation reproduces our expectation precisely.No error bars are shown in this figure because, without interaction, the density is independent of the pairing field.Therefore, no sampling is needed, and the calculation of the observable is numerically exact.
For the interacting case we study the density in units of the non-interacting density n 0 , i.e., n/n 0 in the balanced case of µ ↑ = µ ↓ ≡ µ.In Fig. 2, we show densities for a range of interaction strengths, together with the corresponding analytical results.The error bars For the data points marked with an "x" lattice effects are extrapolated out, while those marked with an "o" are obtained using a maximum likelihood estimator.
represent an estimate of the standard deviation of the MC samples obtained through Jackknife resampling [83] to reduce effects of autocorrelation in the samples.For chemical potentials βµ ≳ −1, the calculated densities start to exhibit a significant dependence on the parameter C I which determines the size of the lattice.To remove these lattice artefacts, we have performed an extrapolation of the data in this regime.These lattice artefacts are extrapolated out in this βµ-regime.For βµ ≲ −1, the dependence of the calculated densities on the lattice size is negligible and a maximum likelihood estimator was used to determine the density instead.We refer the reader to App.A for a more detailed discussion.
We observe that that the calculated densities are in excellent agreement with the analytic solution.We emphasize that this is non-trivial.For example, in the mean-field approximation, the density equation of state agrees identically with the one of the non-interacting gas in this coupling regime, i.e., n/n 0 = 1 for all values of βµ. 4 Interestingly, we find that the observed dependence of the density on βµ is qualitatively very similar to the one found in one-and two-dimensional fermion models of this type, see, e.g., Refs.[58,[84][85][86].Note that, in the unitary limit (in three dimensions), the density equation of state increases monotonically as a function of βµ and does not develop a maximum around βµ = 0, see, e.g., Refs.[87,88].

Pairing
In principle, we can also search for non-trivial ground state configurations in terms of the pairing field.For example, in Subsec.4.1, we have computed the effective potential U in the mean-field approximation which develops a non-trivial ground state associated with 〈φ〉 ̸ = 0 for β g > 4. Such a non-trivial ground state would be associated with spontaneous symmetry breaking and the formation of superfluid gap.Of course, in the present zero-dimensional model, the appearance of a superfluid ground state is simply a failure of the mean-field approximation.In the following we are not interested in a demonstration of the breakdown of the mean-field approximation in this strong-coupling regime.We rather aim at a demonstration that expectation values of the pairing field (or even products of multiple pairing fields as associated with correlation functions) are indeed directly accessible in our approach.To be specific, we can study the behavior of the pairing field throughout the MC sampling process by evaluating the observable φ, φ = λ which can be calculated from the dimensionless discrete field configuration vector φ as φsample Since we have 〈φ〉 = ψ ↓ ψ ↑ , i.e., the expectation value of two annihilation operators, we expect to find a sample mean of zero for the observable φ.In Fig. 3, we show the sampling process for the observable φ for one simulation run.As it should be, we find that the mean of φ = 0 lies within the deviation of the result, already without correcting systematic sources of uncertainty (which originate from the CL time spacing and the lattice size).
Beyond the expected expectation value of φ, we can also see that the sampling process is well behaved and does neither "get stuck" nor exhibit "large jumps" that can occur when singularities are present in the drift term of the CL equation.Such singularities in the drift indeed form a potentially serious impediment for the CL method.A review of this aspect can be found in Ref. [43], including an illustration of some of the problems emerging from the presence of singularities in the drift term with the aid of a 0 + 0-dimensional theory (i.e., a theory where the fields depend neither on time nor space).In our present work, however, we consider a spacetime formulation of the problem, such that even the 0 + 1-dimensional case considered in our numerical studies already features high-dimensional integrals, rendering an exact analysis of singularities in the drift (as done for the 0+0-dimensional case in Ref. [43] and also in various ways in Refs.[89][90][91][92][93][94]) impractical.We therefore only analyze the sampling process, which we find to be localized to a region around the origin in the field space, and use the indicator observable defined in Eq. ( 63) to spot potential singularities.We find that the drift indeed seems to be free of singularities in our present study, as already indicated above.However, for studies with a finite number of space dimensions, the issues associated with singularities in the drift may have to be carefully analyzed on a case-by-case basis.Note that the localization of the sampling process in our formulation also aides us performing the sampling process more efficiently, indicating that the theory is a good candidate for treatment with the CL method.

Sign problem
Our system of interest exhibits a sign problem, even in the absence of a spin or mass imbalance, which follows directly from the properties of the matrix M in Eq. ( 36) when evaluated on nonuniform field configurations.To surmount this sign problem, we have chose the CL approach to MC sampling.As a consequence of that, however, both the real and imaginary part of the pairing field are themselves complex quantities.This results in a complex integrand of the path integral for all interacting systems.In this subsection, we would like to briefly study this aspect by analyzing the MC samples of the weight of the path integral e −S B .
In Fig. 4, we present distribution histograms of the complex argument of the path integral weight e −S B for different values for the chemical potential βµ and the interaction λ.We generally observe narrow distributions around an argument of zero.Thus, the weights are close to the real axis but exhibit a weak phase problem.The width of the phase distribution increases   65) and the fitted parameters are shown in the legend.The width of the argument distribution increases with increasing interaction with an exponent that depends on the chemical potential βµ. with increasing interaction strength, which may not come unexpected.This is illustrated in Fig. 5, where the standard deviations of the phase distributions from Fig. 4 are shown as a function of the interaction parameter λ together with a fit of the data to the following model: We have chosen this model because the system is known to be free of a sign problem in the noninteracting case associated with λ = 0 and f (λ) = 0 coincides with the expected distribution width of zero.

Conclusions
In this work we successfully developed a lattice formulation of the pairing field formulation of two-component Fermi gases interacting via a two-body interaction.We have shown that our formalism leads to the well-established continuum theory in the limit of large lattices.With our lattice formulation at hand, we employed the CL approach for the computation of densities for a 0 + 1-dimensional systems.We chose such a low-dimensional model since exact results are available to benchmark our approach.Moreover, studies of low-dimensional models may be viewed to be computationally less intense, at least at first glance.In any case, the use of the CL approach is convenient to deal with the sign problem which is present in our approach even in the absence of, e.g., spin and mass imbalances.In the non-interacting limit, our simulation reproduces the analytic solution exactly.Very importantly, in the phenomenologically relevant interacting case, the results from our stochastic approach are found to be in excellent agreement with the analytic solution, providing a proof of concept for this formalism.Note that we have also analyzed the sign problem in our studies which we found to be mild, at least for all coupling strengths considered in our present work.For a study of the phase structure and critical behavior of, e.g., ultracold Fermi gases, the computation of correlation functions is indispensable.As we have already demonstrated in terms of the expectation value of the pairing field, there is no conceptional issue within our present approach which would hinder the computation of such functions.Of course, a study of the formation of long-range order or at least quasi long-range order is only meaningful  in higher-dimensional models, d > 0. We have already presented the formalism, even for d > 0. However, explicit calculations of equations of state and correlation functions for these phenomenologically more models are deferred to future work.

Figure 1 :
Figure 1: Dimensionless density for a non-interacting system, i.e. λ = 0.The numerical results are in perfect alignment with the analytical solution.

Figure 2 :
Figure 2: Ratio of interacting and non-interacting density for interacting systems with λ ∈ {0.25, 0.5, 0.75, 1.0}.The errorbars represent the statistical uncertainty of the data points.For the data points marked with an "x" lattice effects are extrapolated out, while those marked with an "o" are obtained using a maximum likelihood estimator.

Figure 3 :Figure 4 :
Figure 3: MC samples points of φ in the complex plane.The color indicates at which CL time the sample point is calculated; the simulation ends at t CL,max .The sample mean of the cloud, indicated by the lines, lies at −0.22(03) − i 0.20(04) with the standard deviation of the mean obtained through Jackknife resampling and indicated by the shaded areas.The ellipse shows the principal standard deviations of the point cloud, i.e. the square roots of the eigenvalues of the covariance matrix of the cloud, with the semiaxes being the standard deviation values in length oriented along the eigenvectors of the covariance matrix.The simulation was carried out with λ = βµ = 1.0 with C I = 150.0resulting in a lattice with N τ = 150.

Figure 5 :
Figure 5: Sample standard deviation of the samples shown in Fig. 4. The dashed lines show fits of the data to the model in Eq. (65) and the fitted parameters are shown in the legend.The width of the argument distribution increases with increasing interaction with an exponent that depends on the chemical potential βµ.

Figure 6 :
Figure 6: Illustration of the extrapolation of finite-lattice effects in the density observable at the point λ = βµ = 1.The fit takes the uncertainties of the single simulation runs into account.This results in the shown uncertainty band around the extrapolated value and represents the standard deviation of the estimate.