# Thermodynamics, contact and density profiles of the repulsive Gaudin-Yang model

###### Abstract

We address the problem of computing the thermodynamic properties of the repulsive one-dimensional two-component Fermi gas with contact interaction, also known as the Gaudin-Yang model. Using a specific lattice embedding and the quantum transfer matrix we derive an exact system of only two nonlinear integral equations for the thermodynamics of the homogeneous model which is valid for all temperatures and values of the chemical potential, magnetic field and coupling strength. This system allows for an easy and extremely accurate calculation of thermodynamic properties circumventing the difficulties associated with the truncation of the thermodynamic Bethe ansatz system of equations. We present extensive results for the densities, polarization, magnetic susceptibility, specific heat, interaction energy, Tan contact and local correlation function of opposite spins. Our results show that at low and intermediate temperatures the experimentally accessible contact is a nonmonotonic function of the coupling strength. As a function of the temperature the contact presents a pronounced local minimum in the Tonks-Girardeau regime which signals an abrupt change of the momentum distribution in a small interval of temperature. The density profiles of the system in the presence of a harmonic trapping potential are computed using the exact solution of the homogeneous model coupled with the local density approximation. We find that at finite temperature the density profile presents a double shell structure (partially polarized center and fully polarized wings) only when the polarization in the center of the trap is above a critical value which is monotonically increasing with temperature.

###### pacs:

67.85.-d, 02.30.Ik## I Introduction

The Gaudin-Yang model Ga1 (); Y1 () which describes one-dimensional spin- fermions interacting via a delta-function potential has a long and distinguished history being the subject of theoretical investigations for more than fifty years GBL (). This experimentally realized model MSGKE (); LRPLP () not only presents an incredibly rich physics of which we mention Tomonaga-Luttinger Gia () and incoherent spin Luttinger liquids BL (); Ber (); ISLL1 (); ISLL2 (), Bardeen-Cooper-Schrieffer BCS1 (); BCS2 () and Fulde-Ferrel-Larkin-Ovchinnikov like pairing, spin-charge separation, Fermi polarons McGuire () and quantum criticality and scaling QC (), but is also amenable to an exact solution which allows for a parameter free comparison of theory and experiment. From the historical point of view the Gaudin-Yang model, also known as the two-component Fermi gas (2CFG), was the first multi-component system solved by the nested Bethe ansatz (BA) (the particular cases of one and two particles with spin up in a sea of opposite spins were considered in McGuire (); FL1 ()). The thermodynamics of the 2CFG in the framework of the thermodynamic Bethe ansatz (TBA) was derived independently by Takahashi and Lai Tak1 (); Tbook (); Lai (). The experimental advances of the last 20 years in the field of ultracold atomic gases which allowed for the creation and manipulation of 1D multicomponent systems renewed the interest in such models which were investigated further using a large variety of exact and approximate methods Tokatly (); Astrak1 (); Fuchs (); IidaW (); Orso (); HLD (); ColTat (); LGSB ().

At zero temperature the groundstate in the thermodynamic limit of the Gaudin-Yang model is characterized by a system of two Fredholm integral equations whose solution allows to derive the phase diagram Orso (); HLD (); ColTat (). The situation at finite temperature is much more complicated. This is due to the fact that the application of TBA produces an infinite system of nonlinear integral equations which are very difficult to investigate even numerically. This is the main reason why the vast majority of results regarding the thermodynamic behaviour of the 2CFG found in the literature are restricted to a small region of the relevant parameters (temperature, chemical potential, magnetic field and coupling strength). Even though numerical schemes to treat the TBA equations have been developed TS (); CKB (); KC (); CJGZ () the need for an efficient thermodynamic description of multicomponent 1D systems cannot be overstated.

In this paper we address this problem by deriving a system of only two integral equations characterizing the thermodynamics of the repulsive Gaudin-Yang model which is valid for all values of the physical parameters and from which physical information can be easily extracted numerically. We employ the same method we have used in the case of the two-component Bose gas KP1 (); KP2 () which utilizes the fact that the 2CFG can be obtained as the scaling limit of the Perk-Schultz spin chain with a specific grading PS1 (); Schul1 (); BVV (); dV (); dVL (); L () and the quantum transfer matrix MS (); Koma (); SI1 (); K1 (); K2 (). The largest eigenvalue of the quantum transfer matrix gives the free energy of the associated lattice model which means that performing the same scaling limit we obtain the grandcanonical potential of the continuum model. In the homogeneous case we present extensive results for the densities, polarization, susceptibility, specific heat, Tan contact and local correlation function of opposite spins for a wide range of coupling strengths, temperatures and chemical potentials. We find that, for and any value of the polarization, the contact Tan1 (); OD (); BP1 (); ZL (); WTC (); BZ (); VZM (); WC () is a nonmonotonic function of both coupling strength and temperature. The implication of this unusual behavior can be more easily understood if we take into account that for delta-function interactions the contact controls the tail of the momentum distribution via with . The local maxima or minima of the contact result in abrupt changes of , which can be experimentally detected, and can be used for accurate thermometry. The reconstruction of the momentum distribution as a function of the temperature for the balanced impenetrable system was first described by Cheianov, Smith and Zvonarev in CSZ ().

The experimentally relevant situation in which the system is subjected to an external harmonical potential is treated using the local density approximation KGDS2 (); Orso (); HLD (); ColTat () and the solution of the homogeneous system. Compared with the zero temperature case ColTat (), at low-temperatures the density profiles present a double shell structure with partially polarized center and polarized wings only if the polarization at the center of the trap is above a critical value which depends on the temperature. As we increase the temperature eventually the entire system becomes partially polarized.

The plan of the paper is as follows. In Section II we introduce the Gaudin-Yang model and present the zero and finite temperature Bethe ansatz solution. The new thermodynamic description of the model is introduced in Section III. Numerical data for the densities, polarization, susceptibility, specific heat and Tan contact are reported in Sections IV and V. The density profiles are investigated in Section VI. The derivation of the NLIEs from the similar result for the Perk-Schultz spin chain is performed in Sections VII, VIII and IX. The solution of the generalized Perk-Schultz model can be found in Appendix A.

## Ii The Gaudin-Yang model

The Gaudin-Yang model Ga1 (); Y1 () describes one-dimensional spin- fermions interacting via a delta-function potential and represents the natural extension in the fermionic case of the Lieb-Liniger model LL (). In second quantization the Hamiltonian reads

(1) |

where is the trapping potential, is the chemical potential, the magnetic field and we consider periodic boundary conditions on a segment of length . are fermionic fields which satisfy canonical anticommutation relations In (1) the symbol denotes normal ordering. The coupling constant can be positive or negative corresponding to repulsive or attractive interactions. In this paper we will consider only the repulsive case Due to the spin-independent interaction the Zeeman term in the Hamiltonian (1) is a conserved quantity. In the following it will be preferable to use units of and introduce the dimensionless coupling constant where is the density of the system and the scattering length. In terms of this coupling constant the strong coupling regime is defined as and the weak coupling regime as .

The Gaudin-Yang model is exactly solvable only when . However, in the case of a sufficiently shallow trap, the system can be considered locally homogeneous and the local density approximation (LDA) Orso (); HLD (); ColTat (); KGDS2 () can be used. The LDA coupled with the solution of the homogeneous system allows for accurate predictions of the density profiles which will be computed in Section VI. Therefore, we will first present the solution in the case. For a system with particles of which have spin up and have spin down the energy spectrum of the homogeneous system is Ga1 (); Y1 ()

(2) |

where the charge and spin rapidities satisfy the Bethe Ansatz equations (BAEs)

(3a) | |||

(3b) |

The solution at finite temperature is a very difficult task to accomplish even though the system is integrable. Assuming the string hypothesis, Takahashi and Lai Tak1 (); Tbook (); Lai () obtained for the grandcanonical potential ( with the temperature) with satisfying the following infinite system of nonlinear integral equations (NLIEs)

(4a) | |||

(4b) | |||

(4c) |

together with the asymptotic condition . In Eqs. (4), which are also known as the thermodynamic Bethe ansatz equations, , with . It is clear that extracting physically relevant information from the TBA equations is very hard even from the numerical point of view. In general, numerical implementations require the truncation of the equations after a certain level, an approximation which introduces uncontrollable errors especially in the intermediate and high-temperature regime.

## Iii Efficient thermodynamic description of the repulsive Gaudin-Yang model

A similar situation is encountered in the case of the one-dimensional two-component repulsive Bose gas (2CBG). Using a lattice embedding and the quantum transfer matrix in KP1 (); KP2 () we derived a simple system of two NLIEs characterizing the thermodynamics of the 2CBG which circumvents the problems associated with the TBA equations. The same method can be applied in the case of the Gaudin-Yang model (the derivation is presented in Sections VIII and IX) obtaining for the grandcanonical potential per length

(5) |

where are auxiliary functions satisfying ()

(6a) | ||||

(6b) |

and

(7) |

We would like to stress that this system of equation is exact and valid for all values of chemical potential, magnetic fields (including the balanced case ), positive coupling strengths and temperature. We also note that (6) differs from the result for the 2CBG derived in KP1 () by dropping the diagonal terms.

In the left panel of Fig. 1 we present results for the temperature dependence of the grandcanonical potential of a system with and showing how for increasing values of the coupling strength we approach Takahashi’s result for impenetrable particles Tak1 (); Tbook ()

(8) |

A comparison between the grandcanonical potential for a system with and computed using Eq. (5) and the TBA equations (4) truncated at the level can be found in the right panel of Fig. 1. The results show almost perfect agreement which not only confirms the validity of our results but also supports the string hypothesis in the case of the repulsive Gaudin-Yang model.

In the noninteracting limit, , using it is easy to show that (5) reduces to which, as expected, is the result for free fermions in a magnetic field.

## Iv Thermodynamic properties of the homogeneous system

Despite being the first multi-component model solved by Bethe ansatz at zero and finite temperature our knowledge of the thermodynamic properties of the Gaudin-Yang model is still very limited. Fueled by the interest in the FFLO pairing phase the case of attractive interactions has received more attention Orso (); HLD (); HJLPAD (); BAD (); GYFBLL () but even in this case a complete characterization of the thermodynamical properties is still lacking. The situation is even more dire in the repulsive case. The phase diagram at was derived by Colomé-Tatché ColTat () and the low-temperature strong-coupling regime was investigated by Lee, Guan, Sakai and Batchelor in LGSB ().

The main reason behind this information scarcity is the complexity of the TBA equations and other methods used in the investigation of thermodynamic properties, like lattice Monte-Carlo calculations, which require rather cumbersome numerical schemes in order to obtain accurate numerical data. In contrast, our equations (6), can be accurately and easily solved employing a simple iterative procedure with the convolutions treated using the Fast Fourier Transform and the convolution theorem. This scheme allows us to probe a wide region of the parameter space with the exception of the low polarization at low-temperatures regime ( and ). This regime requires a more careful treatment of the NLIEs which will be deferred to a future publication.

In Fig. 2 we present for several values of the magnetic field the dependence on the reduced temperature of the densities , polarization , magnetic susceptibility and specific heat which can be derived from the grandcanonical potential (5)

At zero temperature and zero magnetic field the groundstate of the Gaudin-Yang model is balanced (number of particles with spin up is equal with the number of particles with spin down) and the excitation spectrum is gapless which means that switching a magnetic field the system will become partially polarized. The monotonously decreasing polarization of the system as a function of the reduced temperature in the presence of a constant magnetic field is presented in the upper panels of Fig. (2). We can also see that for a given and temperature the polarization is an increasing function of the interaction strength. For large values of the dimensionless coupling strength the system becomes paramagnetic as it can be seen from Eq. (8) which describes free fermions at chemical potential .

The magnetic susceptibility is finite everywhere except at and . In this case we can see from Eq. (8)) that the magnetization is and the susceptibility diverges like at vanishing magnetic field. At low-temperatures and small magnetic fields presents a complex behavior as a function of . For the maximum susceptibility is obtained for but already for the susceptibilities for present more pronounced maxima at very low-temperatures while the maximum for moves at a higher temperature. For strong magnetic fields the susceptibility is zero at low-temperatures, presents a global maximum at a reduced temperature close to the value of the magnetic field and the dependence on the coupling strength becomes small. The specific heat is linear at low-temperatures and similar to the case of the susceptibility presents a maximum at a temperature . A similar behavior was noticed by Klauser and Caux KC () for the 2CBG but in our case the maxima are more pronounced for the same values of temperatures and magnetic field. At high temperatures the specific heat per particle approaches , the value for the ideal gas. As a function of the coupling strength and for low values of the magnetic field the specific heat reaches its maximum for but this is no longer true for strong magnetic fields where the dependence on becomes less pronounced.

## V Tan contact and local correlation functions

S. Tan Tan1 () discovered that the momentum distribution of 3D fermions with zero range spin-independent interactions exhibits a universal behavior at large momenta (a similar relation for the Lieb-Liniger model was derived earlier by Olshanii and Dunjko OD ()). The large momentum distribution is the same for both particle species and the constant , called contact, is given by the probability that two particles of opposite spin are found at the same position in space. The contact also appears in Tan’s adiabatic theorem Tan1 () which determines the change of the energy with respect to the interaction strength and a series of universal thermodynamic identities called Tan relations Tan1 () (see also BP1 (); ZL (); WTC (); BZ (); VZM (); WC ()). The analogues of Tan relations in 1D were derived by Barth and Zwerger BZ () using the Operator Product Expansion of Wilson and Kadanoff Wils (); Kad (). The contact is not only an important and valuable theoretical concept but is also experimentally measurable NNCS (); SGDJ (); KHLDMD (); KHHDHV2 (); KHDHHV1 (); SDPJ (); WMPCJ (); HLFHV ().

In the case of the Gaudin-Yang model with the Hamiltonian (1) the contact is defined as BZ () (remember and with the scattering length)

(9) |

where by we denote the zero or finite temperature expectation value. In our definition we can see that is an extensive variable and is closely related to the interaction energy and the local correlation function of opposite spins defined by

(10) |

A simple application of the Hellmann-Feynman theorem gives the 1D version of Tan’s adiabatic theorem BZ () which expresses the variation of the total energy with respect to the scattering length.

The situation is considerably simplified in the homogeneous case . The local correlation function and density no longer depend on and we can introduce the contact and interaction energy per length

(11) |

Another application of the Hellmann-Feynman theorem allows us to obtain the local correlation function from the derivative of the grandcanonical potential per length with respect to the coupling strength . Noting that with we have which proves the identity. For the homogeneous system the pressure and total energy per length satisfy the following Tan relation BZ () . We have checked numerically this relation and found perfect agreement providing an additional check for the validity of our NLIEs. Other analytical or numerical computations of the contact in 1D systems can be found in HJLPAD (); BAD (); RPLD (); YB (); ZWKGD (); DK (); BBF (); VM ().

### v.1 Results at T=0

At zero temperature the thermodynamic properties of the system can be
extracted from the following system of Fredholm type integral equations Y1 ()
^{1}^{1}1 An analytical derivation of Eqs. (12) from the low-temperature limit of our NLIEs, Eqs. (6), is still lacking.
The main difficulty lies in the fact that our (complex) auxiliary functions are not easily related to the dressed energies of the TBA
formalism.

(12a) | ||||

(12b) |

where In Eqs. (12) and are parameters which fix the values of the density of spin down particles and energy per length via and The balanced system is characterized by and the fully polarized system by . Even though in the general case an analytic solution of the system of equations (12) is not known, in certain limits approximate results can be derived. In the strong interaction limit () Guan and Ma GM1 () obtained (the leading order term was derived in Fuchs (); BCS2 ())

(13) |

with the Riemann function, and in the weak interaction limit

(14) |

It should be emphasized that Eqs. (13) and (14) are not exact. They represent numerical approximations which are valid only in the Tonks-Girardeau () and weak interaction () regimes. Using Tan’s adiabatic theorem (see above) the local correlation function of opposite spins, contact and interaction energy can be computed as

(15) |

In the following it will be preferable to work with the dimensionless contact density defined as with .

In Fig. 3 we present results for the zero temperature energy and dimensionless contact as functions of the coupling constant obtained from the numerical integration of Eqs. (12) and the asymptotic expansions (13) and (14). In the balanced case () the asymptotic expansions for energy and contact represent a good approximation in the weak and strong interacting regime but in the case of imbalanced systems () the results derived from Eq. (13) for the strongly interacting regime become more accurate as the polarization approaches one. The contact is a monotonously increasing function of and in the Tonks-Girardeau limit () approaches asymptotically a finite value. As a function of the polarization reaches a maximum for and is zero for the fully polarized system. For the balanced system in the strongly interacting regime, using (13) and (15) we obtain

(16) |

with a result which was first obtained in BZ ().

### v.2 Results at finite temperature

As we have mentioned before at finite temperature analytical results are restricted to the low-temperature strong-coupling regime LGSB () and therefore we will have to rely on numerical solutions of the NLIEs (6). In Fig. 4 we present the dependence on the interaction strength of the dimensionless contact, local correlation function and interaction energy for fixed polarization and several values of the reduced temperature.

An interesting feature revealed by our data is that the contact at low and intermediate temperatures (see upper panels of Fig. 4) presents a local maximum (for any value of the polarization) which does not appear at and gets suppressed at high temperatures. In the Tonks-Girardeau regime the contact approaches a finite value for all values of . Remembering that the contact governs the tail of the momentum distribution ( for ) this means that at and fixed polarization (or magnetic field) the width of the momentum distribution increases monotonically with reaching its maximum at In contrast, at low and intermediate temperatures the width of the momentum distribution reaches a maximum at a finite value of the coupling strength and then becomes narrower as This reconstruction of the momentum distribution as a function of the strength of the interaction should be in principle experimentally accessible.

As shown in the middle panels of Fig. 4 the local correlation function which quantifies the probability that two particles of opposite spin occupy the same point in space is a monotonically decreasing function of with a maximum in the noninteracting limit ( for and all ). In the strongly interacting limit drops to zero like . For a fixed value of as a function of the polarization the local correlator is maximal in the balanced system and vanishes in the fully polarized case. The interaction energy (bottom row Fig. 4) exhibits a local maximum at zero and finite temperature and vanishes for . This is due to the fact that in the strong interaction limit the short-range repulsive interaction acts as an effective Pauli exclusion principle also between fermions of opposite spins BZ (). As expected, the interaction energy is largest for the case and it is zero for the fully polarized system.

An alternative way of deriving the contact requires the analytical or numerical computation of the momentum distribution which is the Fourier transform of the static Green’s function i.e., In the general case this is an almost hopeless task due to the fact that our knowledge of the correlators is still rather limited BL (); Ber (); ISLL1 (); IP1 (); GIKP (); GIK (); GKa (). However, in the impenetrable case we can use the Fredholm determinant representations derived by Izergin and Pronko IP1 () which have the advantage of being extremely easy to implement numerically using the method presented in Bor (). In our case the Green’s functions have the following representation IP1 ()

(17) |

where the integral operators and act on an arbitrary function like and with kernels

(18) |

and is the Fermi weight. In addition we have The results for the contact obtained from the large analysis of the Fourier transform of (17) are presented in the panels of Fig. 4 as red disks (), blue squares (), black triangles () and green diamonds () and they are in perfect agreement with the results derived from the NLIEs.

The temperature dependence of the contact shown in Fig. 5 reveals another interesting phenomenon (in Fig. 5 the values of the contact at were computed from the integral equations (12) and for from the NLIEs (6)). For weak interactions the dependence of on temperature is very small, but as we increase the coupling strength, the contact develops a local minimum which becomes very pronounced in the Tonks-Girardeau limit. The nonmonotonic behavior is present for all polarizations and the temperature interval in which it manifests itself shrinks as the value of increases. Outside of this temperature interval and for strong interactions the contact presents a linear dependence on the reduced temperature for all values of The nonmonotonic behavior of the contact implies that the shape of the momentum distribution of the strongly repulsive Gaudin-Yang model suffers an abrupt change in a very small interval of temperature. This reconstruction of the distribution is also counterintuitive because as the temperature increases the distribution becomes narrower and not wider as one would expect from a system of weakly interacting or free particles. However, one should not forget that in this case we are dealing with a strongly interacting system where this type of picture is not correct.

This crossover of the momentum distribution was first discovered by Cheianov, Smith and Zvonarev in CSZ () for the balanced system and signals the transition from the Tomonaga-Luttinger liquid phase to the incoherent spin Luttinger liquid regime. Our results show that a similar crossover is present also in the imbalanced case and it might be possible to be detected even for moderate values of the coupling strength. Compared with the single-component case in the two-component system there are two temperature scales and (remember ) which define two different quantum regimes: characterized by the Tomonaga-Luttinger liquid theory and in which the incoherent spin Luttinger theory is applicable. The two temperature scales are well separated in the strong coupling limit and for the charge degrees of freedom are effectively frozen while the spin degrees of freedom are strongly “disordered”. It is important to note that in the computation of the correlation functions the limits and do not commute ISLL1 (). The determinant representation (17) has been derived by taking the limit in the wavefunctions followed by the summation of form factors at finite temperature. If we take the limit in (17) and compute the dimensionless contact from the tail of the momentum distribution we obtain which is considerably smaller than (see Eq. (16)) obtained from the integral equations at zero temperature (Eqs. 12). The momentum crossover for the balanced impenetrable system is presented in Fig. 1 of CSZ ().

## Vi Density profiles at finite temperature

In most experiments the system is subjected to a harmonic potential , a situation in which the Hamiltonian (1) is no longer exactly solvable. However, for slowly varying potentials we can apply the Local Density Approximation (LDA) and the solution to obtain reasonably accurate data for the trapped system. Under LDA each point in the trap can be seen as a locally homogeneous system with

(19) |

where and are the chemical potential and magnetic field in the center of the trap. An immediate consequence of (19) is that the density along the trap is monotonously decreasing and at zero temperature vanishes at a distance from the center of the trap where is called the Thomas-Fermi radius and is determined by . An inhomogeneous system can be characterized by three parameters: the dimensionless coupling strength , the polarization and the reduced temperature all evaluated at the center of the trap.

The density profiles at were computed in ColTat () and it was noticed that for all they present a two shell structure: an imbalanced mixture of spin up and down fermions in the center and fully polarized wings. At finite temperature the situation becomes more complex as it can be seen in Fig. 6 where we present the density profiles at for different values of the coupling strength and two polarizations. For and the entire system is in the mixed phase for all values of (the density is the same in all cases presented) even though the polarization along the trap increases away from the center. The double shell structure appears only for values of the center polarization above a critical value which is a monotonic increasing function of and temperature (see the lower panels of Fig. 6). At high temperatures the entire system will be found only in the mixed phase irrespective of the value of and .

## Vii The Gaudin-Yang model as the continuum limit of the Perk-Schultz spin chain with grading

The derivation of Eqs. (5) and (6) is performed in three steps. First, we identify a lattice embedding of the Gaudin-Yang model for which the quantum transfer matrix method MS (); K1 () can be employed. The need for this lattice model stems from the fact that the quantum transfer matrix method, which has the advantage of producing a finite number of NLIEs for the thermodynamics of the system, can be defined only for discrete systems. In the second step we calculate the free energy of this lattice model from the largest eigenvalue of the QTM. Finally, the thermodynamics of the continuum model is obtained by performing the scaling limit in the result for the lattice model.

As in the case of the two-component Bose gas KP1 (); KP2 () the appropriate lattice model is the critical Perk-Schultz spin chain dVL (); L ()

(20) |

where is the coupling strength, is the number of lattice sites, is the anisotropy, are chemical potentials and are the grading parameters. Also with and the canonical basis and the unit matrix in the space of -by- matrices. If in the case of the 2CBG we considered the grading for the Gaudin-Yang model the relevant grading is . The energy spectrum of the spin chain is (see Appendix A)

(21) |

with the Bethe ansatz equations

(22a) | |||

(22b) |

In order to prove that the Perk-Schultz spin chain with the grading is the lattice embedding of the 2CFG we are going to show that the spectrum and Bethe equations of the continuum model (2), (3) can be obtained in a scaling limit from the lattice analogues (21), (22). We start with the Bethe ansatz equations. Performing the transformation and with the lattice constant and a small parameter, Eqs. (22) take the form

Considering we obtain

(23a) | |||

(23b) |

If we take the limits , such that , , identifying with and using

we see that for even and odd Eqs. (23) transform in the Bethe equations of the Gaudin-Yang model (3).

Under this set of transformations the energy spectrum becomes

(24) |

If we denote by the inverse temperature of the continuum model and consider , such that is finite we have

(25) |

where is the canonical partition function of the spin chain and