Complex Langevin study for polarons in an attractively interacting one-dimensional two-component Fermi gas

We investigate a polaronic excitation in a one-dimensional spin-1/2 Fermi gas with contact attractive interactions, using the complex Langevin method, which is a promising approach to evade a possible sign problem in quantum Monte Carlo simulations. We found that the complex Langevin method works correctly in a wide range of temperature, interaction strength, and population imbalance. The Fermi polaron energy extracted from the two-point imaginary Green's function is not sensitive to the temperature and the impurity concentration in the parameter region we considered. Our results show a good agreement with the solution of the thermodynamic Bethe ansatz at zero temperature.


I. INTRODUCTION
The quantum Monte Carlo method [1,2] is widely used in various fields of physics as a non-perturbative tool of analysis.In a path integral formalism using Lagrangian, a partition function is written in terms of the integral of the Boltzmann weight e −S over field variables, where S is an action.When the action is a real-valued function, the Boltzmann weight is regarded as a probability density function.This ensures that quantum expectation values of physical observables can be estimated by importance sampling of the Boltzmann weight.
However, the positivity of the Boltzmann weight is violated in many physically interesting systems: Hubbard model, finite density quantum chromodynamics (QCD), QCD with a θterm, matrix superstring models and any systems defined by Schwinger-Keldysh formalism which describes real-time dynamics, for instance [3][4][5][6][7][8][9][10]. In these cases, the number of samples becomes exponentially large as the system size grows in order to obtain statistically significant results.In non-relativistic fermionic systems, a frequently used way to apply the quantum Monte Carlo method is introducing bosonic auxiliary field through the Hubbard-Stratonovich transformation [11][12][13].After integrating out the fermion fields, we will obtain an effective action of the auxiliary field.Since the effective action involves a logarithm of a fermion determinant, the positivity is not guaranteed except in a few cases where the action has the particle-hole symmetry [14], the Kramers symmetry [15], or Majorana positivity [16][17][18], for instance.
A promising approach to evade the sign problem is the complex Langevin method [19,20], which is an extension of the stochastic quantization to systems with complex-valued actions.
An advantage of this method is that it is scalable to the system size, and thus, the computational cost is similar to the usual quantum Monte Carlo method without the sign problem.
On the other hand, it is known that this method sometimes gives incorrect answers even when the statistical average of a physical observable converges.In the recent decade, a way to judge the reliability of the complex Langevin method is extensively studied [21][22][23][24][25][26][27][28][29][30][31], and proposed criteria which are able to compute in actual simulations using the boundary terms [21,22,29,30] and the probability distribution of the drift term [23,25].While it is still difficult to predict when the complex Langevin method fails without performing numerical simulations, we can eliminate wrongly convergent results thanks to these criteria.
In this paper, we consider a spatially one-dimensional spin-1/2 polarized fermions with contact attractive interactions which is known as the Gaudin-Yang model [40], and compute the single-particle energy of spin-down fermions in a spin-up Fermi sea, which is referred to as the Fermi polaron energy.Recently, the single-particle excitation spectra of Fermi polarons were experimentally measured in higher dimensional atomic systems [41][42][43][44][45][46][47][48] (also see a recent review [49] for Fermi polarons).While an analytic formula for the polaron energy in one dimension is obtained exactly at zero temperature based on the thermodynamic Bethe ansatz method [50], no analytical solutions are known at finite temperature (note that the numerical results for finite-temperature properties of polarized gases within the thermodynamic Bethe ansatz were reported in Ref. [51]).The one-dimensional Fermi polarons were studied with several theoretical approaches such as Bruckner-Hartree-Fock [52], T-matrix [53,54], and variational [55,56] approaches.In this study, we demonstrate that a microscopic quantity, that is, the polaron energy is efficiently computed by the complex Langevin method.This paper is organized as follows.In Sec.II, we derive a lattice action of the Gaugin-Yang model.In Sec.III, we review how to compute physical quantities using the complex Langevin method.In Sec.IV, we show a way to extract the ground state energy in the spin-down channel from a two-point imaginary time Green's function.In Sec.V, we present the numerical results.Section VI is devoted to the summary of this paper.In this work, k B and are taken to be unity.

II. THE GAUDIN-YANG MODEL
We consider a one-dimensional two-component Fermi gas with contact attractive interactions which is known as the Gaugin-Yang model [57,58].The Hamiltonian is given by where ĉp,σ and ĉ † p,σ are fermionic annihilation/creation operators with momentum p and spin σ =↑, ↓, respectively.In this work, the atomic mass is taken to be unity.The coupling constant g is related to a scattering length in one dimension a as g = 2 a > 0 [40].The chemical potential of spin-σ fermions are represented by µ σ .For convenience, we introduce an average chemical potential µ = (µ ↑ + µ ↓ )/2 and a fictitious Zeeman field h = (µ ↑ − µ ↓ )/2.The grand canonical partition function is given by Z = Tr e −β( Ĥ− σ µσ Nσ) with β being an inverse temperature and a number operator Nσ = p ĉ † p,σ ĉp,σ .The path-integral representation of Z reads where action S is given by Here ψ σ (x, τ ), ψ * σ (x, τ ) are a Grassmann field and its complex conjugate.
While the action (3) is given in a continuous spacetime, one should perform a lattice regularization appropriately to carry out numerical simulations.We write lattice spacing of temporal and spatial directions by a τ and a x , respectively, and their ratio by r = a τ /a 2 x .We also introduce lattice quantities by where n and j are integers that satisfy 0 ≤ n < N τ and 0 ≤ j < N x .The inverse temperature and the spatial length of the lattice is given by β = T −1 = N τ a τ and L = N x a x .With these notations, we consider a lattice action: where φj,n is a bosonic auxiliary field.As shown in Ref. [59], the lattice action (5) correctly converges to the continuum one as long as the matching conditions are satisfied, where f k is a ḡ-dependent constant given by In practice, it is sufficient to use an approximated form of the matching conditions which are obtained as the first order approximation in the expansion in terms of a τ .After integrating out the fermion fields, the partition function and the effective action of the auxiliary field read where the effective action of the auxiliary field is given by where I is the N x × N x identity matrix.Since we consider a naive finite difference as an approximation of the second order derivative with respect to x, the eigenvalues of B are . It has been argued in Ref. [11,59] that this naive lattice action converges too slow to the continuum limit, and the behavior can be improved by replacing the eigenvalues of B by After this replacement, the form of B σ is given by A notable point is that the effective action (10) involves a logarithm of the fermion determinant which can be complex in general.Therefore, this term may cause the sign problem if the Zeeman field h is not zero and then the Monte Carlo simulation can be difficult to apply to this system.

III. COMPLEX LANGEVIN METHOD
The complex Langevin method (CLM) [19,20] is an extension of the stochastic quantization which is usually applicable to real-valued actions.In the CLM, we first consider a complexified auxiliary field φn,k and extend the domain of definition of S eff to the complex space.For such a complex field, we consider a fictitious time evolution described by the complex Langevin equation: where η n,k (t) is a real Gaussian noise.When we assume that the system described by the complex Langevin equation reach equilibrium at t = t eq , an average of a physical observable O( φ) can be defined as with the average over the noise O( φη (t)) η being We note that η n,k (t)η n ,k (t ) η = 2δ nn δ kk δ tt , in particular.Although, we expect that the mean value O( φ) is equivalent to the quantum expectation value calculated in an original action, i.e., n,k d φn,k O( φ)e −S eff / n,k d φn,k e −S eff in the limit ∆t → 0, it is not correct in general.There are extensive studies [21][22][23][24][25][26][27][28][29][30][31] to understand when the CLM is justified, and criteria for determining whether a CLM is reliable or not have been proposed.One of a practical criterion which can be relied on in actual numerical simulations is discussed from a view point of a probability distribution of a drift term [23,25].In our case, it is sufficient to consider a magnitude of the drift term given by and its distribution.According to the criterion, the CLM is reliable if the probability distribution of v η shows an exponential decay.

IV. OBSERVABLES
The number density of spin-σ fermions is given by The particle number density on a lattice unit is defined by nσ = n σ a x .From below, we assume that the spin-down fermions are regarded as minority.Typical temperature and momentum scales are given by the Fermi scales which are determined by the density of spin-up fermions: In lattice simulations, we can compute dimensionless combination T /T F and p F a as follows: In order to calculate the polaron energy, we consider two-point Green's function: where T τ is the imaginary-time-ordered product.Hereinafter we restrict τ > 0. We write the eigenvalue and eigenstate of K by K |n = K n |n .In particular, We also assume that the ground state |0 is not degenerated.Expanding the trace by the eigenstates, the correlation function reads where ∆K n ≡ K n − K 0 .In the low temperature limit β → ∞, only the ground state contributes to the summation over n.Thus, we find Since the matrix elements appeared in the above expression does not depend on τ , it behaves like where A 0 , A 1 , . . .are τ -independent constants, and E 0 , E 1 , . . .are energies of the ground state and excited states.In particular, the energy of the ground state can be extracted by keeping τ β.The polaron energy U is defined by The polaron energy is the shift of single-particle energy from that in the case of free fermions due to the interaction between majority (spin-up fermions) and minority (spindown fermions).We note that the polaron energy at zero temperature is calculated exactly based on thermodynamic Bethe ansatz [50]: In lattice calculations, the polaron energy is obtained as follows.From the form of the effective action, the lattice expression of the inverse Green's function reads From a straightforward algebra, each component of G reads The momentum representation of G ij is calculated by the discrete Fourier transformation.
Therefore, if the temporal lattice size N τ is sufficiently large, G(p, τ ) is approximately given by V. NUMERICAL RESULTS We performed complex Langevin simulations on (N τ , N x ) = ( 40 In every Langevin step, the magnitude of the drift term ( 17) is calculated and stored, and finally the probability distribution P (v η ) of drift term can be drown.In Fig. 1, a typical result of the probability distribution P (v η ) is shown.It is normalized so that the integral of the distribution is 1.In each simulation, we confirmed that P (v η ) shows a linear falloff in the log-log plot.This means that P (v η ) shows an exponential fall-off, and then our calculation of CLM was reliable.We also investigated the eigenvalues of the matrix which is the reduced matrix of inverse Green's function Eq. ( 28) appearing in the effective action on the lattice Eq. ( 10) as an effective fermionic matrix.We calculate the eigenvalues w σ of the matrix from one configuration in the case of several values of βh = 0 to βh = 12 and other fixed parameters, N τ = 40 and βµ = −1.2.The imaginary part of w σ is negligibly small and hereinafter we discuss the real part of the eigenvalues.The numerical results of eigenvalues log Re[w σ ] of the matrices G −1 ↑,red and G −1 ↓,red are shown in Fig. 2. Red circles and blue squares correspond to the eigenvalues of G −1 ↑,red and G −1 ↓,red , respectively.In the case of βh = 0, w ↑ is exactly same as w ↓ because G −1 ↑,red = G −1 ↓,red .While the range of the eigenvalues of G −1 ↑,red tend to be broad, the range of the eigenvalues of G −1 ↓,red tend to be narrow when βh increases.
It is notable point of this eigenvalue-analysis that the eigenvalues of G −1 σ,red are always larger than 1 because of log(Re[w]) > 0 even in the case of the large βh, corresponding to the large population imbalance.This result indicates that the integrand of the partition function ( 9) is always positive and no sign problem happens in the parameter region of our calculations.Note that this is a numerical finding in our setup, and we do not prove that the sign problem never occurs in the Gaudin-Yang model with population imbalance.We In Fig. 3, we show dimensionless quantities T /T F , 1/p F a and n ↓ /n ↑ , which are typical indicators of the temperature, the interaction strength and the population imbalance, respectively.The ratio of particle numbers n ↓ /n ↑ becomes significantly small when βh 1 (βµ ↑ βµ ↓ ) as expected.In that case, T /T F and 1/p F a are also small since T F and p F are proportional to n ↑ .For each parameter, we computed the ratio of Green's functions R(0, na τ ) defined in Eq. ( 25) at zero momentum.Numerical results on an N τ = 40 lattice at βµ = 0 for a single configuration are shown in Fig. 4. Qualitative behavior of R(0, na τ ) at other N τ and βµ are same as these results.In the parameter region we swept, R(0, na τ ) has a plateau at intermediate imaginary time, which suggests that the energy spectrum is gapped from any possible excited states.In our analysis, we extract the single-particle ground-state energy E 0 (0) by because the long-time limit as Eq. ( 25) cannot be taken on a lattice.
After calculating the single-particle energy (33), the polaron energy U is obtained by Eq.For a fixed N τ , the numerical results show similar behavior to Eq. ( 27) as a function of 1/p F a despite they also depend on T /T F and n ↓ /n ↑ .Moreover, the numerical results tend to be close to the exact result at zero temperature when we take the continuum limit.Our result suggests that the polaron energy is insensitive to the temperature and the impurity concentration.

VI. SUMMARY
We have studied excitation properties of Fermi polarons at finite temperature for the attractive Gaudin-Yang model with large population imbalances using the complex Langevin method, a non-perturbative approach free from the sign problem.We have performed numerical simulations for several chemical potential βµ and Zeeman field βh, the dimensionless control parameters of the model, and found that our simulation covers wide range of temperature T /T F , strength of the coupling 1/p F a and population imbalance n ↓ /n ↑ .We have computed the polaron energy as a function of 1/p F a. While our result is still away from the zero temperature and single-polaron limit, the computed polaron energy shows similar (1/p F a)-dependence to the exact result at those limits.
The complex Langevin method well works in Gaudin-Yang model even in the presence of the population imbalance.Practically, within our setup, the probability distribution of the drift term always show an exponential fall-off, which means that the problem of wrong convergence does not occur.Moreover, the integrand of path integral is always positive within our simulation from the eigenvalue-analysis.However, it is known that the sign problem is severe in the case of higher dimension [60].Thus the behavior of the probability distribution of the drift term and the eigenvalues in higher dimension will be investigated as future study.
One interesting application of the complex Langevin method is to study the transition from degenerate Fermi-polaron regime to classical Boltzmann-gas regime of a unitary spinimbalanced Fermi gas which is found to be a sharp transition by a cold-atom experiment using 6 Li Fermi gases in a three-dimensional box potential [61].Also, it is interesting to explore an inhomogeneous pairing phase [35,36] and in-medium bound states [62], which cannot be addressed by quantum Monte Carlo simulation due to the sign problem in the mass-and population-imbalanced systems.In order to discuss such phenomena, we need more elaborate estimation of systematic errors.The work in this direction will be presented elsewhere.

( 26 )
. In Fig. 5, we show the polaron energy on N τ = 40 and 80 lattices.As the temporal lattice size N τ becomes large, the system is close to the continuum limit.The color of each point represents the statistical average of T /T F .The lowest temperature is T /T F 0.08 for N τ = 40 and T /T F 0.07 for N τ = 80, respectively.The ratio of particle numbers n ↓ /n ↑ varies from 1.0 × 10 −5 to 1.5 × 10 −1 for N τ = 40 and 8.4 × 10 −8 to 1.0 × 10 −1 for N τ = 80, respectively.The solid line indicates the exact result at zero temperature shown in Eq. (27).

12 FIG. 4 .
FIG.4.The τ -dependence of the ratio of Green's functions R(0, τ ) on an N τ = 40 lattice at βµ = 0. n denotes the number for the discretized imaginary time τ = na τ .Each point in this plot is obtained for a single configuration.

FFIG. 5 .
FIG. 5.The polaron energy computed by the CLM.The solid line is the exact result at T = 0 shown in Eq. (27).