Superfluidity of fermionic pairs in a harmonic trap. Comparative studies: Local Density Approximation and Bogoliubov-de Gennes solutions

Experiments with ultracold gases on the lattice give the opportunity to realize superfluid fermionic mixtures in a trapping potential. The external trap modifies the chemical potential locally. Moreover, this trap also introduces non-homogeneity in the superconducting order parameter. There are, among other approaches, two methods which can be used to describe the system of two-component mixtures loaded into an optical lattice: the Local Density Approximation (LDA) and the self-consistent Bogoliubov–de Gennes equations. Here, we compare results obtained within these two methods. We conclude that the results can be distinguishable only in the case of a small value of the pairing interaction.

( ) , where r 0 is the center of the trap. The interaction term is treated within the mean-field approximation:

Local density approximation solution
The LDA method is based on an assumption that, locally, the system behaves like a uniform gas [1,2]. The parameters of the system depend solely on the value of the particle density = +   n n n i i i (i.e. the average number of particles), at each point i of space [13,14]. Hence, the local SOP is defined only by is the single-particle dispersion relation of the one-dimensional system. Here, it should be noticed that the presented derivations do not depend on the lattice dimensionality. Using anticommutation relation, the first term can be rewritten as: This allows us to write down the Hamiltonian  MF in the Nambu notation: are the Nambu spinors. The system can be solved using the Bogoliubov transformation: where g † and γ are the quasiparticle operators. The spectrum of the Hamiltonian is given by l =  + D   k 2 2 , while Figure 1. Schematic representation of fermionic gases on the lattice in a trapping potential. The atoms with spins  and  (orange and blue color, respectively) can tunnel between nearest-neighbor lattice sites with the hopping amplitude t. The on-site interaction strength U between atoms can also be tuned via Feshbach resonances. The trapping potential introduces an effective site-dependent chemical potential and non-homogeneous distribution of atoms in the system.
Hence, the coefficients u k and v k satisfy the condition: . The energy of the system is determined from the grand-canonical potential W º -- k T ln Tr exp k T

B MF B
( [ ( )]), while the SOP value is determined from its global minimum. In the absence of a discontinuous phase transition, the minimum of energy corresponds to ¶W ¶D = 0 or, equivalently: The above equation is called the gap equation.
can be found as: . As a result, one can find values of m D( ) and m n( ) within the LDA method.

Bogoliubov-deGennes solution
In the case of the non-homogeneous Hamiltonian(3), the system can be solved within the Bogoliubov-Valatin transformation: where g s n and g s n † are the quasi-particle fermionic operators and s s = - . The above transformation can be treated as an extended version of the transformation(7) in real space. If the Hamiltonian is diagonalized under the transformation(11), i.e., g g = å s s s s   MF n n n n † , the relations [17]: hold such that one straightforwardly finds Bogoliubov-deGennes equations [18] (called also Blonder-Tinkham-Klapwijk equations): is the single-particle Hamiltonian and d D = D U ij i ij are the on-site SOPs. Due to the symmetry of the BdG equations with respect to spin inversion, we can find an identity relation between eigenpairs: he SOPs can be found self-consistently as:

Numerical implementation
As has been pointed out in the previous paragraph, the chemical potential m 0 of the system in the trapping potential is fixed by the total number of particles N. Next, the solution of the system within the LDA method is defined by the effective potential μ i given by the Thomas-Fermi equation. Hence, the corresponding SOP Δ can be found by minimization of W D ( )at fixed μ i [19][20][21]. In fact, this way is equivalent to the one in which equation (9) should be solved [22]. The huge advantage of this method is based on the fact that it is relatively simple to implement and quite quickly solvable even for a trapped system, due to an efficient procedure of finding the solution for a homogeneous system. For instance, in the two-dimensional (2D) case, the solution (the minimum of the grand canonical potential) can be found for the system with around a thousand sites [21].
In turn, the self-consistent solution Δ i (15) and s n i (16) of the BdG equations (13), in the presence of the trapping potential, can be found only when all eigenpairs are known. Hence, the BdG equations are associated with the exact diagonalization of the Hamiltonian in matrix form, which can be found numerically using one of the available numerical libraries, e.g. LAPACK [23], Intel MKL [24] or MAGMA [25,26]. However, this method is very computational time-consuming ∼N 3 and leads to a limitation in the system size N. This problem can be partially solved by using the sparse matrices storage implemented in some modern softwares, for instance ARPACK [27]. Moreover, there exist some iterative methods for solving the BdG equations [28][29][30][31][32], based on the kernel polynomial method [33]. For instance, in the Chebyshev-Bogoliubov-de Gennes method, the Chebyshev expansion of Green's functions is implemented. This type of methods allows to study the large size systems, even with hundreds or thousands sites [32]. Moreover, in some physical problems, it is reasonable to limit the number of eigenpairs to only few nearest ones to 'zero energy' (i.e. the Fermi level). This issue is implemented, e.g., in FEAST [34] and it can be applicable to topological edge states investigations, for instance.

Spin-balanced system
Here, we describe the balanced system case, which means the average number of particles with opposite spin in each site is equal ( =   n n i i ). Due to the fact that the average number of particles n i in each site is determined by the effective chemical potential μ i , to perform the calculation within both methods, the chemical potential is discretized with the parametrization: m m = + V r i i ( ). Figure 2 shows n versus μ and Δ versus μ dependences, obtained within LDA calculations. In the case of a one-dimensional system, the band edges are at m  = t 2 . At U/t=0, the chemical potential changes fromt 2 to t 2 , which corresponds to the filling n=0 and n=2 at the band edges. The SOP is equal to zero for the noninteracting system. With increasing on-site interaction, the SOP increases. For intermediate couplings, the chemical potential drops below the lower band edge, which, according to the Leggett criterion [35], indicates the BCS-BEC crossover [36]. Found values of the n(μ) and Δ(μ) can be used as a solution of the system for a given U<0 from the LDA approach. Figure 3 shows results in the case of the system in a harmonic trap, We fix =´-V t 15 10 0 6 . The average number of particles n i , obtained within these two methods, is comparable for the chosen parameters. A similar behavior is observed for the SOP Δ i in the center of the trap. However, the situation looks differently at the edge of the harmonic trap (insets). Namely, LDA gives significantly different solutions for the SOP from those obtained within the BdG method at relatively small attractive interactions U  (top panel). A similar behavior has been observed in the case of a spin-imbalanced system in which unconventional pairing [37][38][39] can be realized at the edge of the system.
It is worth emphasizing that the LDA is a reliable approach when the gradient corrections are relatively small. In turn, the BdG solutions in each site depend on the rest of the system via equations (13) and (15). Hence, the BdG method is a better approach to describe the system under consideration. Moreover, the LDA can also fail in the description of the system in a strong anisotropic trap [40].

Other applications and limitations
The theoretical predictions within the LDA method are in relatively good agreement with experimental results [41][42][43][44]. This method was already successfully implemented in the context of trapped ultracold quantum gases [37,40,[45][46][47]. However, the LDA method has some limitations which can be lifted by applying the BdG technique. Below, we give two examples of the systems for which the LDA method fails.

Spin-imbalanced systems
Both methods under consideration can be applied for a spin-imbalanced system, i.e., when ¹   n n i i [13,48]. This corresponds to the experimental setup with a two-component Fermi mixture, with unequal number of atoms with opposite spins [49][50][51]. This leads effectively to the substitution  å - where h is the Zeeman magnetic field inducing the spin-imbalance. In the case of the LDA method, this is equivalent to the substitution . Hence, an additional condition is needed in the numerical calculations, i.e. d = å - Here, the -  n n i i corresponds to the magnetization at i-th site. It is worth mentioning, that the symmetry of the BdG equations, given by equation (14), is preserved.
In the weak coupling limit, at a large spin imbalacne (i.e. polarization) [22,38], the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase [52,53] can occur. In this case, the SOP is characterized by oscillations in real space and at non-zero local polarization [13,54]. If the polarization is equal to zero or low, the connection between the phase diagrams of the system without a trapping potential and in the presence of the harmonic trap [55] exists, in the LDA method [22]. Unfortunately, the LDA approach cannot account for the oscillating character of the SOP as well as the oscillations of the average number of particles in the FFLO phase [13,56]. A similar problem has been observed for the spin-imbalanced gap [57] estimated within the LDA method and the Density Matrix Renormalization Group (DMRG). In this case, the oscillating character of n i at relatively low polarization has not been reproduced within the LDA method. Moreover, the LDA calculations underestimate critical parameters of the system, with respect to the results obtained within the DMRG [57] or the BdG approach [22].  The situation can be even more complicated in the two dimensional systems. In the case of the harmonic trap, the additional SOP oscillations along the radial direction have been reported [37,[58][59][60]. However, for the toroidal trap, the oscillations of the SOP can be function of angle, due to rotational symmetry of the system [61][62][63]. The combination of these two types of oscillations can be expected [64] as well. This case shows an advantage of the BdG technique in the comparison to LDA, i.e. the spin-imbalanced systems with unconventional superconducting phases.

Non-homogeneous systems
The strongly inhomogeneous systems, e.g. systems with long-range magnetic or charge orderings or disordered systems are examples for which the LDA method fails. Apart from the above systems, one should mention ultracold fermionic gases or the condensed matter systems [65]. Here, the impurities [66][67][68] or vortices [69][70][71] play the role of inhomogeneity. It is worth mentioning the superconducting systems with spin-density (SDW) or charge-density (CDW) waves coexistence [72][73][74][75]. The 'local' disorder in the system influences on the properties of the system. Because of the presence of inhomogeneity with 'long-range' impact on the system, the LDA method can not be used. In such cases, only the BdG technique is suitable [65].

Summary
In this paper, we studied an ultracold fermionic mixture loaded into a one-dimensional lattice, put in a harmonic trap. We briefly compared the Local Density Approximation and the self-consistent Bogoliubov-deGennes equations-the main assumptions, implementations and limitations of both of these methods. We have shown that both approaches give comparable results in the case of relatively large paring interactions U. For smaller U, the results for the superconducting order parameter Δ can be significantly different, especially in the case of a small density of particles, i.e. small filling at the edge of the system. Concluding, the results obtained within the LDA method can be treated as a starting point for the more accurate BdG technique. One should be aware of the LDA method limitations. This method is strongly unrecommended to use for spin-imbalanced as well as strongly inhomogeneous systems.