{'fraction_non_alphanumeric': 0.03326711074753613, 'fraction_numerical': 0.0120369288301975, 'mean_word_length': 4.618734952943751, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ground-state properties of the one-dimensional attractive Hubbard model with confinement: a comparative study

We revisit the one-dimensional attractive Hubbard model by using the Bethe-ansatz based density-functional theory and density-matrix renormalization method. The ground-state properties of this model are discussed in details for different fillings and different confining conditions in weak-to-intermediate coupling regime. We investigate the ground-state energy, energy gap, and pair-binding energy and compare them with those calculated from the canonical Bardeen-Cooper-Schrieffer approximation. We find that the Bethe-ansatz based density-functional theory is computationally easy and yields an accurate description of the ground-state properties for weak-to-intermediate interaction strength, different fillings, and confinements. In order to characterize the quantum phase transition in the presence of a harmonic confinement, we calculate the thermodynamic stiffness, the density-functional fidelity, and fidelity susceptibility, respectively. It is shown that with the increase of the number of particles or attractive interaction strength, the system can be driven from the Luther-Emery-type phase to the composite phase of Luther-Emery-like in the wings and insulating-like in the center.


I. INTRODUCTION
The Hubbard Hamiltonian is a simplified but important prototype model for strongly correlated electrons in solid state. The one-dimensional (1D) version is usually used as a basic model examining the electronic correlation behavior in the quasi-one-dimensional systems exhibiting the itinerant character of electrons and superconductivity 1,2 like the organic superconductor (TMTSF) 2 X. 3 Due to its solvability by the Bethe-ansatz method, 4 this model serves as a popular benchmark to test the accuracy of different approximate methods in handling the ground-state (GS) properties of strongly correlated fermions especially in an intermediate coupling strength. [5][6][7] This model has recently arisen renewed interest, in ultracold atomic gases, where fundamental many-body phenomena can be investigated in a wellcontrolled way.
The ultracold atomic gases loaded on optical lattices have opened an exciting field in simulating condensedmatter model Hamiltonians like the boson and fermion Hubbard models. More recently a Tonks-Girardeau gas of bosonic 87 Rb atoms was realized experimentally in a 1D optical lattice. 8 In the case of fermionic atom gases, however, it is more difficult to cool them down because of the Pauli exclusion principle different from the bosonic atoms. Up to now, a two-component Fermi gas of 40 K atoms with tunable interacting strengths has been prepared in a quasi-1D geometry. 9 An interacting Fermi gas of 40 K atoms has been demonstrated in three-dimensional optical lattices. 10 Within present-day techniques, 11 it is possible to trap interacting gases of fermionic atoms in 1D optical lattices and then simulate the 1D Hub-bard model experimentally. At the same time, experimental achievements on population-imbalanced ultracold fermions 12,13 make it possible to systematically study the exotic paring states in 1D fermion systems with or without optical lattices. 14 Moreover, the square-well trap and optical box-like trap have been experimentally realized in the atomic system, which makes it meaningful to study the physics of 1D uniform atomic gases with a hard wall directly. 15 Considering the strongly correlated nature of 1D systems and the enhanced quantum fluctuations, it is crucial to have techniques beyond the mean field theory, both analytically and numerically. Among them, Bethe-ansatz solution, 4 quantum Monte-Carlo (QMC) simulation, [16][17][18][19][20][21] density-matrix renormalization group (DMRG), [22][23][24][25] and some sophisticated theoretical methods, such as the dynamical mean field theory combined with the numerical renormalization group 26 and a time-evolving block decimation numerical study 27 , are listed as well-established techniques. On the other hand, the exact diagonalization scheme 23,28 has also been used to obtain the exact numerical solution for the small size systems as an ultimate technique.
The inhomogeneity caused by the boundary conditions or the external potentials, normally invalidates a reliable analytical method that is usually used in the homogeneous system. On the other hand, the numerical schemes mentioned above, e.g., the exact Betheansatz solution incorporated with local density approximation is proved to be useful in obtaining the phase diagram and the collective oscillations of the atomic mass density in the presence of harmonic traps. 29 The lattice version of the Bethe-ansatz based density-functional theory (BADFT) 25,30-32 is also a good candidate in studying inhomogeneous strongly-correlated system. For the strongly interacting regime, some non-perturbative methods, like, the Fermi-Bose mapping and the Betheansatz method, have to be used. 4,[33][34][35] In general, the Bethe-ansatz solution provides reliable results beyond the mean-field theory. 36 BADFT has been used in studying the 1D Gaudin-Yang gases 37,38 and 1D Hubbard model, 31,39 which is different from the density-matrix functional theory with the basic variable to be the singleparticle density matrix. 40 The rich non-uniform phases induced by the trapping potential can be identified by the variance of the local density, the local compressibility, 21 the double occupancy, 41 or the thermodynamic stiffness. 42 The entanglement entropy is found by França et al. to be a powerful measure characterizing the phase-separated states in spatially inhomogeneous environment. 43 Recently, a new concept named fidelity in quantum-information and quantum-computation theory, is used to characterize quantum phase transitions by computing the overlap between two states of nearby points in parameter space. By measuring the similarity between these two states, it is expected to show a minimum in the fidelity due to a dramatic change in the ground-state wavefunctions around the critical point of the quantum phase transition. The studies are made both experimentally and theoretically on many model systems, such as on the 1D XY model, the Dicke model, and the 1D Hubbard model. 44 Gu further proposed a density-functional fidelity to measure the similarity between density distributions of two ground states in parameter space based on the observation that the density distributions 45 can reflect the change in the ground-state wave functions. An interesting question then arises: Can this method be applied to quantify the composite quantum phases transitions of the 1D Hubbard model in the presence of confinements?
The 1D Fermi-Hubbard model with attractive on-site interactions, has long served a "minimal" model that best describes superconductivity. For the fermionic atomic gases loaded on optical lattices, the strength or even the sign of the on-site interaction between different species can be easily tuned by using a technique named Feshbach resonance. 46 Marsiglio calculated the ground-state energy (GSE) and energy gap to the first excited state for the attractive Hubbard model using the variational canonical Bardeen-Cooper-Schrieffer (CBCS) and grand CBCS wave function. 47 This model for different concentrations of electrons is also studied within a selfconsistent field method. 7 In this work, motivated by earlier theoretical work and by ongoing experimental efforts, we revisit the 1D attractive Hubbard model in the weak-to-intermediate coupling strength under different boundary conditions with or without trapping potential. We present a fully microscopic theoretical study on the ground-state properties.
We begin with introducing the Bethe-ansatz solution of the model and solving the coupled equations. In Sec-tion III, we provide a Bethe-ansatz based spin-densityfunctional theory (SDFT) and its local-spin-density approximation. In Section IV, we show our main numerical results, compared to the exact Bethe-ansatz, the CBCS 48 , and DMRG ones. Our comparisons are done on general band fillings and a wide range of coupling strength in three different systems, i.e., 1D lattices with periodic boundary condition, hard-wall, and harmonic confinement. The energy, the binding gap, and the pairbinding energy of the ground state are in a good agreement with the exact results. We also provide a simple and direct diagnostic of the quantum phase transition by calculating the thermodynamic stiffness, density-functional fidelity, and fidelity susceptibility. At last, our results are summarized.

II. THE 1D ATTRACTIVE FERMI-HUBBARD MODEL
We consider a 1D attractively interacting twocomponent Fermi gas with N f atoms in a lattice with unit lattice constant and N s lattice sites, which can be described by a one-band Fermi-Hubbard model, 49 where σ =↑, ↓ is a pseudospin-1/2 label for two internal hyperfine states,n i = σn iσ = σĉ † iσĉiσ is the total site occupation operator, t is the tunneling between nearest neighbors, U is the strength of the on-site attractive interaction, and V iσ describes the trapping potential. We consider in this paper 1D lattices with periodic boundary condition, in a simple box-shaped trapping potential (V iσ = 0 inside the box and V iσ = +∞ elsewhere, namely, a hard-wall) and a parabolic confine- amplitude of the confining potential). The numbers of atoms with spin up and spin down are N ↑ and N ↓ , respectively.
The 1D attractive Fermi-Hubbard model in the absence of confining potentials (i.e.,Ĥ s =Ĥ ref ) belongs to the universality class of Luther-Emery liquids, characterized by a massive spin gap and zero charge gap. At zero temperature, the properties ofĤ ref in the thermodynamic limit (N σ , N s → ∞ with finite N σ /N s ) are determined by the filling n σ = N σ /N s and by the dimensionless coupling constant u = U/t. According to Lieb and Wu 4 , the ground state ofĤ ref for different fillings in the thermodynamic limit is described by the coupled integral equations for the momentum distributions ρ(k) and σ(λ), with The non-negative parameters Q and B are determined by the normalization conditions The ground-state energy per site of the system is written in terms of the momentum distributions as, which is a function of filling n, magnetization s, and interaction strength u. The GS properties of the 1D attractive Hubbard model have been analyzed by many authors. 7,50

III. LATTICE SPIN-DENSITY-FUNCTIONAL THEORY
Density functional theory (DFT) 51,52 is a powerful tool to calculate the GS properties of a many-body Hamiltonian including inhomogeneity. Through Kohn-Sham mapping, the true system of interacting particles is mapped onto a fictitious one in which the particles do not interact but reproduce the true particle density. The central problem in DFT is the exchange-correlation (xc) functional, which remains elusive and can be simply treated by the so-called local-density approximation (LDA). 53 The essence of LDA is to locally approximate the xc energy of the interacting inhomogeneous system based on that of an interacting homogeneous reference system. In discrete lattice systems, there is a latticeversion DFT, the so-called site-occupation functional theory (SOFT), introduced in the pioneering papers by Gunnarsson and Schönhammer to study the band-gap problem in the context of ab initio theories of fundamental energy gaps in semiconductors and insulators. 30 Within SOFT the GS site occupation can be obtained by solving self-consistently the lattice Kohn-Sham (KS) equations together with the closure relation Here the sum α runs over all the occupied orbitals, and the effective KS potential is given by The first term in the Kohn-Sham potential comes from the Hartree mean-field contribution, while V xc iσ [n iσ ] is the derivative of the xc energy E xc [n ↑ , n ↓ ] evaluated at the GS site occupation. The total GSE of the system is given by The local-spin-density approximation (LSDA) for the spin-polarized system has been shown to provide an excellent account of the GS properties over a large variety of inhomogeneous systems. In this work we employ the following Bethe-ansatz based LSDA functional (BALSDA), where the xc potential of the 1D homogeneous Hubbard model is defined by 39 Here, ǫ GS (n, s, u) is the GS energy per site of the 1D homogeneous system defined in Eq. (6) and ǫ GS (n, s, 0) = −(4t/π) sin(πn/2) cos(πs). Making use of the chain rule, we obtain an explicit expression for the exchangecorrelation potential in terms of the chemical potential µ(n, s, u) and magnetic field h(n, s, u), where the upper (lower) sign refers to the spin-up (spindown) atomic species. The chemical potential and magnetic field are defined by the first derivative of the GSE of the homogeneous system with respect to filling n and magnetization s as, h(n, s, u) = ∂ ∂s ǫ GS (n, s, u) . In Fig. 1, we illustrate the numerical results for the xc potential as a function of filling n for different magnetization s. Obviously, the xc potential is continuous in the whole range of n and there is no charge gap for the system of attractive interactions, which is different from the repulsive case where a discontinuity occurs and a charge gap opens at n = 1. 39 In Fig. 2, we plot the xc potential as a function of s for various values of n with u = −4. We see that there is always a discontinuity at s = 0, which corresponds to a gap opened in the spin channel signalling a Luther-Emery liquid phase of paired particles. We define the xc gap ∆ xc,σ as the discontinuity in V hom xc,σ (n, s, u) at s = 0, ∆ xc,σ (n, u) = lim s→0 + V hom xc,σ (n, s, u) − lim s→0 − V hom xc,σ (n, s, u) .
The results are shown in Fig. 3 for different fillings. We notice there is a symmetry between ∆ xc,↑ (n, u) and ∆ xc,↓ (n, u), as a result of the symmetries of the GSE ǫ GS (n, s, u) = ǫ GS (n, −s, u) and our choice for the Hartree energy. As a consequence, the xc potential in Eq. (10) possesses a discontinuity at s = 0 and naturally the information of the Luther-Emery nature and the spin energy gap in the reference system is transferred to the inhomogeneous system through Eq. (10).

IV. NUMERICAL RESULTS
In this section we present the numerical results obtained from the self-consistent solution of Eqs. (7)-(8) based on the BALSDA of Eq. (10). We illustrate our main numerical results, which are summarized in Figs. 4-10. We compare the results obtained with the SDFT/BALSDA to results of DMRG and CBCS. In the use of DMRG, the number of states kept is set maximum 1200 and the truncation error is less than 10 −10 .
First, in Figs. 4-6, we show the results for systems of periodic boundary condition. In Fig. 4, the GSE per site (absolute value) is plotted as a function of the particle density for u = −2 and u = −4. Good agreement is seen in the whole range of fillings between the SDFT results and the exact Bethe-ansatz ones, while CBCS calculations deviate in higher density. In Fig. 5, we show the GSE as a function of interaction strength for the system of even and odd number of particles. The comparison shows that SDFT results give poorer agreement for small system size with even number of particles than CBCS but better agreement in other situations, such as for larger system size or stronger coupling strength. We should keep in mind that the reasons for which the SDFT results deviate from the exact Bethe-ansatz ones are due to both finite size effects and the xc functional built in the KS potential, while for the CBCS results they are affected by both finite-size effects and lack of correlation. The ground-state energy and the energy gap (in Fig. 6, see below) are less accurate by SDFT in small size systems and by CBCS at intermediate coupling strengths.
For the 1D Hubbard model of attractive interaction, one defines the gap or binding energy as, 47,54 which describes the energy difference between two systems, the one of two subsystems with N f particles each, and the other of N f + 1 and N f − 1, respectively. From the viewpoint of DFT, this gap has two contributions, the one from the KS gap, and the other from the xc gap, related to the discontinuity of GSE as the number of fermions varies across an integer. 55 This definition is also explained as the negative of the excitation gap, which equals to the half of the energy it takes to break a pair in the system. 56 In Fig. 6 Quantitatively we find that a better agreement is achieved by SDFT for larger system size and stronger attractive interaction. The gap is negative for odd number of particles, and positive for even number of particles, which reflects the fact that the system with negative interactions favors pairing. 47,48,57 This is the even-odd effects discussed by Tanaka and Marsiglio. 48,58 From the right panel of Fig. 6, we notice that the energy gap oscillates as a function of even particle numbers (depending on whether the number of particles is 4m or 4m + 2), which is the super-even effect (for details, see Ref. 48).
The confined system composed of the hard wall, which  can be realized in experiments with the atom chips, focused laser beams and trapped impurity atoms, 59 breaks the inhomogeneity but still can be solved through the Bethe-ansatz method. 60 In Figs. 7 and 8, we show the SDFT results for the GSE compared to DMRG ones. The system is of N s = 32 lattice sites with u = −1 and −4. We find that, in almost all the range of particle density and weak-to-intermediate interaction strength, SDFT results give an excellent agreement with DMRG calculations. Small deviations happen around the half filling. Next, we discuss a 1D system of Fermi gases trapped in a harmonic potential, which is a system that breaks the homogeneity and the integrability of the model and can not be exactly solved. A quasi-classical approach called the local density approximation is usually used to study the density distribution profile when the external potential is slowly varying. 61 Here and in the following the DFT based on LDA is used. We would like to mention that for the system of even number of particles and spinunpolarized, the DFT is used. But for the system of odd number of particles, the SDFT has to be used. Our results are shown in Table I and II and Figs. 9 and 10. In Table I, we show the GSE for small systems of lattice size N s = 100 and large systems of N s = 300 for the weakto-intermediate interaction strength. The corresponding pair-binding energies are calculated accordingly, which are defined as the energy difference between two systems, the one with two subsystems of particles N f and N f + 2, the other with two subsystems of N f particles each, Our calculation shows that the pair-binding energies are negative, which illustrates that the system likes to accom- modate pairs and has the tendency to form stable dimers. Furthermore, we calculate the four-particle binding energies, 62 A bound state of two particles is confirmed by observing that, in Table II, E p (N f ) < 0 and E 4 (N f ) > 0, which shows the evidence of a paired superfluid. We conclude that, the system induces effective attractions between particles leading to pairs formation (E p (N f ) < 0) and effective repulsions between pairs (E 4 (N f ) > 0) avoiding clusters formation. The same calculation is done for the 1D Fermi-Hubbard model with repulsive interactions. It is found that E p (N f ) > 0 and E 4 (N f ) > 0 and no paired superfluid is formed. The comparison of the GSE between SDFT and DMRG shows that the agreement is up to the second decimal digit, except for few cases, such as for N f = 30 and u = −2 with the relative error of 1.2%. The harmonically trapped system described by Eq. (1) allows for two different spatial coexistent quantum phases. The one is the phase characterized by the Luther-Emery (LE) liquid combining with the appearance of the atomic density wave (n i < 2), and another is characterized by a composite phase with a plateau in the trap center of tightly bound spin-singlet dimers (n i = 2) surrounded by Luther-Emery layers (n i < 2). We call the latter the band-insulating (BI) phase. These two different phases can be well classified by a scaled dimensionless variable, the characteristic densityρ = N f V 2 /t and the interaction strength u in theρ-u phase diagram 21 . Thus for certain u, the quantum phases in the presence of confinement can be driven from one to another by increasing ρ. Here we study this phenomena by investigating the thermodynamic stiffness, which is defined as the inverse of the global compressibility, where µ is the global chemical potential. The stiffness has been successfully used in the confined system to characterize the quantum phase transitions. The different coexisted phases appeared in the harmonically confined system are well separated by the nonanalyticity point in the stiffness. The increase of stiffness is related to the incompressible nature of the insulating phases. 42 For the system we studied, the phase is of Luther-Emery type for smallρ 29 and a composite phase of Luther-Emery-liquid type in the wings and insulting type in the bulk for largẽ ρ. We test for a system in a confinement of strength V 2 /t = 10 −3 and N s = 200 lattice sites. The comparison of SDFT with DMRG shows that, in Fig. 9, the stiffness calculated from SDFT is quantitatively correct for weak interaction strength of u = −1 and qualitatively correct for u = −2. The large deviation happens for the system of large number of particles. We find by DMRG,  . We notice that the bulk compressibility or the stiffness calculated here is not a real singular quantity since it includes the sizable contributions from the liquid wings of the density profile. 63 Finally we study the phase transition driven by the attractive interaction strength through the densityfunctional fidelity, which is defined as 45 here n i (u) is the density of particles in lattice site i with interaction strength u. If we fix the distance between u and u ′ , δu = u − u ′ , the density-functional fidelity is expected to show a drop around the critical point due to the maximum distance between the two ground states in different quantum phases. 45 In the top panel of Fig. 10, the density-functional fidelity is plotted as a function of the coupling strength u with the distance between u and u ′ chosen as 0.002 and 0.004. We notice that the fidelity drops abruptly at the critical point u = −1.848, which reflects the fact that the two ground states in the different quantum phases have the maximum distance. In other words, when we fix the number of fermions and the amplitude of the external potential by varying u, this critical value of the attractive interaction is obtained. At that point, the system is in the BI phase with the localized site occupancy of n i = 2 in the center of the trap. The fidelity susceptibility, the leading term of the fidelity, is found to play a central role in the fidelity approach in precisely characterizing the quantum phase transition, which is defined as, where the fidelity susceptibility (FS) χ F takes the form In the bottom panel of Fig. 10, the FS is shown as a function of the coupling strength u, which manifests itself as a marked peak at the same critical point as that in the fidelity function. The sharp peak in the FS implies principally that the density distribution evolves dramatically around the critical point, which is therefore taken as another important quantity to detect the quantum phase transition. 45,64 The peak in the FS is free from the arbitrariness of the small parameter δu. To prove this, we calculated the FS with two different δu in Fig. 10(b). Both cases have good agreements and support the conclusion that the fidelity susceptibility does not depend on the value of δu, but the fidelity does. 65

V. SUMMARY
In summary, in this paper we have studied the groundstate properties of two-component Fermi gases of attractive interactions in a 1D system, by performing SDFT and DMRG simulations. Three different systems have been tested, for fermions on a lattice line with periodic boundary condition, in a hard wall or in the presence of a harmonic confinement. For the system of periodic boundary condition or hard-wall, we showed the calculations for the GSE and the energy gap in accordance with the exact Bethe-ansatz, CBCS, and DMRG results. We found that the ground-state energy is faithfully reproduced by SDFT in the large range of interaction strength and fillings. For the system under the harmonic confinement, we calculated the pair-binding and four-particle binding energies to show that the system of attractive interactions induces effective attractions between particles leading to pairs formation and effective repulsions II: Pair-and four-particle binding energies Ep and E4 per site (in units of t) are calculated with the energies from Table  I for  between pairs avoiding clusters formation. The system is either in the Luther-Emery type phase or in the composite phase of Luther-Emery-like in the wings and insulating like in the bulk. In order to characterize the quantum phase transition driven by the particle number or the attractive interaction strength, the thermodynamic stiffness, the density-functional fidelity, and the fidelity susceptibility are studied, respectively. We found that the thermodynamic stiffness is able to locate the quantum critical point driven by the particle number through the nonanalyticity point, and F (χ F ) is able to locate the quantum critical point driven by the attractive interaction through the abrupt drop (the sharp peak) at the critical point.