Effective Nonlinear Schr\"odinger Equations for Cigar-Shaped and Disk-Shaped Fermi Superfluids at Unitarity

In the case of tight transverse confinement (cigar-shaped trap) the three-dimensional (3D) nonlinear Schr\"odinger equation, describing superfluid Fermi atoms at unitarity (infinite scattering length $|a|\to \infty$), is reduced to an effective one-dimensional form by averaging over the transverse coordinates. The resultant effective equation is a 1D nonpolynomial Schrodinger equation, which produces results in good agreement with the original 3D one. In the limit of small and large fermion number $N$ the nonlinearity is of simple power-law type. A similar reduction of the 3D theory to a two-dimensional form is also performed for a tight axial confinement (disk-shaped trap). The resultant effective 2D nonpolynomial equation also produces results in agreement with the original 3D equation and has simple power-law nonlinearity for small and large $N$. For both cigar- and disk-shaped superfluids our nonpolynomial Schr\"odinger equations are quite attractive for phenomenological application.


Introduction
In the last few years several experimental groups have observed the crossover [1] from the weakly paired Bardeen-Cooper-Schrieffer (BCS) state to the Bose-Einstein condensation (BEC) of molecular dimers with ultra-cold two-hyperfine-component Fermi vapors of 40 K atoms [2][3][4] and 6 Li atoms [5,6]. The unitarity limit of the Fermi superfluid was attained by manipulating an external background magnetic field near a Feshbach resonance which allows an experimental realization of infinitely large value of the s-wave scattering length a [7].
The most interesting feature of the BCS-BEC crossover, noted in experiments [8] on a BCS superfluid as well as demonstrated in theoretical model calculations [9] is that, due to a dominance of the Pauli repulsion among fermionic atoms over the interatomic attraction, the Fermi superfluid remains essentially repulsive in the unitarity limit, i.e. the gas does not collapse and its properties are quite regular. Elaborated Monte Carlo calculations have confirmed this effect [10]. A similar conclusion follows from an examination of the compressibility of a Fermi gas [11]. The Fermi system should exhibit universal behavior in the unitarity limit which should limit the maximum attractive force to a finite value as a → ±∞ [8]. This phenomenon has greatly enhanced the interest in the theoretical study of a Fermi gas in the unitarity limit [12].
Close to the critical temperature a Fermi superfluid can be studied by using the Ginzburg-Landau theory with a complex order parameter [13][14][15]. Recently a nonlinear Schrödinger (NLS) equation for the complex order parameter has been proposed at zero-temperature to study the BCS-BEC crossover [16][17][18][19]. This equation has a simple nonlinear term in the weak-coupling BCS limit (a negative and small), at unitarity (a = ±∞), and the BEC regime (a positive and small) where it becomes the Gross-Pitaevskii equation for Bose-condensed molecules. However, in the full BCS-BEC crossover the 3D NLS equation has a complicated nonlinear term [7,10,17,20]. In addition, this 3D NLS equation must be solved numerically [16][17][18]. Hence an effective one-dimensional (1D) and two-dimensional (2D) reduction of this equation under convenient trapping conditions of cigar and disk-shaped traps is well appreciated. The present paper addresses this important question of 1D and 2D reduction of the original 3D NLS equation at unitarity. Through numerical studies we also investigate the validity of the 3D-1D and 3D-2D reduction of the 3D NLS equation.
In Sec. 2 we consider a confined 3D Fermi superfluid at unitarity and discuss the 3D NLS equation [18,21]. We derive it from an energy functional which takes into account the bulk properties of the system and also the inhomogeneities in the density profile due to the external potential. We show that the gradient term of our energy functional is practically the same of the one recently deduced [22] with an epsilon expansion of energy density in 4 − ǫ dimensions [22,23]. In addition, we find that the total energy calculated from this equation for a 3D system of spin 1/2 fermions, trapped in a spherically-symmetric harmonic potential, is in good agreement [24] with those obtained from accurate Monte Carlo calculations [25,26]. In Sec. 3 we introduce a generic axially-symmetric harmonic potential which models the confining trap of the superfluid fermionic system. In Sec. 4 we consider a cigar-shaped Fermi superfluid. Using a variational ansatz we minimize the 3D energy functional and derive an effective 1D nonpolynomial Schrödinger (1D NPS) equation for this system. Simple analytic forms of the 1D NPS model are derived for small and large number of atoms. The variation of the nonlinear term from small to large number of atoms is illustrated. In Sec. 5 we consider a disk-shaped Fermi superfluid. Using a variational ansatz we minimize the 3D energy functional and derive an effective 2D nonpolynomial Schrödinger (2D NPS) equation for this system. Simple analytic forms of the 2D NPS model are derived for small and large number of atoms. In Sec. 6 we present a numerical study of the different model equations in one and two dimensions for cigar-and disk-shaped systems and compare the results with the full 3D NLS equation. The results of the approximate models are found to be in good agreement with exact calculations. Finally, in Sec. 7 we present concluding remarks.

Fermi superfluid at unitarity
We consider a dilute Fermi gas of N atoms with two equally populated spin components and attractive inter-atomic strength. At zero temperature the gas is fully superfluid and the superfluid density coincides with the total density. A description of the superfluid state can be obtained by using a complex order parameter Ψ(r) conveniently normalized to the total number of superfluid pairs N p [13][14][15], i.e.
The local number density of superfluid pairs is defined as n p (r) = |Ψ(r)| 2 = n(r)/2 where n(r) is the total density (the density of atoms). Notice that Ψ is not the condensate wave function, because the modulus square of it gives the number of particles, not of condensate particles [18,21]. At unitarity (a = ±∞) one finds the condensate density of pairs is n 0 ≃ 0.7 n p , where n p is the total density of pairs [27,28]. In the BEC side (a > 0) of the BCS-BEC crossover n 0 reaches n p quite rapidly by reducing the positive scattering length a, while in the BCS side (a < 0) of the crossover n 0 decreases exponentially to zero by reducing the absolute value of the negative scattering length [27,28]. However, the phase of Ψ(r) is exactly the phase of the condensate wave function and its gradient gives the superfluid velocity [18,21]. Under an external potential U(r) acting on individual atoms, the properties of the Fermi superfluid in the BCS-BEC crossover can be described, in the spirit of the density functional theory [29][30][31], by the energy functional where U p (r) = 2U(r) is the external potential acting on a pair and U(r) the external potential acting on a single atom [16][17][18]. Here E[n p ] is bulk energy density as a function of density of pairs, i.e. the energy density of the uniform system in the BCS-BEC crossover [16][17][18]. The last two terms in (2.2) correspond to the local density approximation (LDA), equivalent to hydrostatics, and the first term involving the gradient corresponds to a correction to LDA [22,31]. Including only the LDA term in energy density leads to the Thomas-Fermi (TF) approximation. In the gradient term m p = 2m is the mass of a pair with m the mass of one atom. This term takes into account corrections to the kinetic energy due to spatial variations in the density of the system. For normal fermions various authors have proposed different gradient terms [32][33][34][35][36]. For superfluid fermions we are using the familiar von Weizsäcker [32] term in the case of pairs.
We notice that the present gradient term is exactly the one that emerges [37] from the Bogoliubov-de Gennes equations in the BEC regime of the BCS-BEC crossover. In this case the energy functional (2.2) is nothing but the Gross-Pitaevskii energy functional of the Bose-condensed molecules [37,38].
In addition, we stress that the gradient term [ 2 /(2m p )] ∇ n p (r) 2 of (2.2) is very close to that recently obtained in [22] for the same gradient term in the case of a Fermi gas at unitarity, in a rigorous calculation using an ǫ expansion of energy density around d = 4 − ǫ spatial dimensions [22,23]. Their final result for this term 0.032 2 [∇n(r)] 2 /[mn(r)], quoted in their (53), can be rewritten as , if we recall m p = 2m and n(r) = 2|Ψ(r)| 2 = 2n p (r) − very close to the present gradient term in (2.2).
In the BCS-BEC crossover, the bulk equation of state of energy density E(n) as a function of density of atoms n of a dilute superfluid depends on the s-wave scattering length a of the inter-atomic potential [16][17][18]. In the unitarity limit, when a → ±∞, the bulk energy density is independent of a and, for simple dimensional reasons [11,[39][40][41], must be of the form where χ is a universal coefficient. Thus, in the unitarity limit the bulk chemical potential is proportional to that of a non-interacting Fermi gas. For fermions, results of fixed-node Monte Carlo calculation give χ = (3π 2 ) 2/3 ξ/2 with ξ = 0.44 [10]. Taking into account (2.3), the energy functional of the Fermi superfluid reads By minimizing (2.4) with constraint (2.1) we obtain the following nonlinear Schrödinger equation where the chemical potential µ 0 is fixed by normalization. This is the zero-temperature 3D NLS equation of superfluid Fermi gas at unitarity. In the case of a large number of fermions, in (2.5) the nonlinear term vastly dominates over the kinetic energy term containing the Laplacian ∇ 2 , and the numerical solution of this equation is thus insensitive to the Laplacian term. In phenomenological treatment of Fermi superfluid, this term is often neglected and the resulting model is called the TF model (LDA approach), which, however, leads to a nonanalytic solution in space variable. The inclusion of the Laplacian term leads to a solution of (2.5) analytic in space variable, although lying close to the solution of the TF model for a large number of particles. The Laplacian term, however, plays an important role in the case of a small number of particles, as we shall show studying the dimensional crossover of the system. Now we compare the energies of our density functional (2.4) with those of Monte Carlo calculations [25,26] in the case of a harmonically trapped system with small number N of spin 1/2 fermions at unitarity. We solve (2.5) numerically, using the Crank-Nicholson method detailed in Sec. 5, and calculate the total energy E of the system using (2.4) for different N in a spherically-symmetric harmonic trap with ω the trap frequency. In figure 1 we plot the total energy E in units of ω vs. N. For comparison, we also plot the TF energy. The energy of the TF approximation (LDA approach) is analytically known to be E(N) = (3N) 4/3 √ ξ/4 [42]. From the data of Figure 1 we find that the average percentage deviation of the present result from the FNDMC data is 3.30%, whereas the same of the LDA result from the FNDMC data is 9.47%. Very recently it has been shown [43] that the best fit to the fixed-node diffusion Monte-Carlo data at unitarity [25] is obtained, in the case of even number of particles, by using λ = 0.18 in the gradient term λ 2 |∇Ψ| 2 /m of the energy functional (2.4).
The usefulness of the present model equation (2.5) over the commonly-used TF approximation (LDA approach) is well appreciated. In fact, the TF approximation is valid only for very large N values when the gradient term is negligible compared to the bulk chemical potential. The present density-functional model for Fermi superfluid at unitarity has been extended over the full BCS-unitarity crossover and the extended model has been found to yield [24] energy of a Fermi superfluid under spherical harmonic confinement in good agreement with Monte Carlo data [44] in the crossover domain.

External confinement: axially-symmetric harmonic potential
If we take r ≡ (ρ,z) so thatρ andz are the radial and axial variables and if we consider the axially-symmetric harmonic trap where λ 2 ω ⊥ and λ 1 ω ⊥ are the trap frequencies in radial and axial directions with λ 1 and λ 2 convenient trap parameters. In the following we shall consider for simplicity a real wave function Ψ(r). Equation (2.5) can be reduced to the following dimensionless form by scalingz = za ⊥ ,ρ = ρa ⊥ , Ψ = ψ N/2/a 3/2 where a ⊥ = /mω ⊥ is the characteristic harmonic oscillator length in the radial direction. Equation For a cigar-shaped superfluid it is convenient to define a linear density through Similarly, for a disk-shaped superfluid it is convenient to consider a radial density through 4. Cigar-shaped Fermi superfluid: 3D-1D crossover

Harmonic Confinement
Now let us suppose that the external trapping potential U(r) is given by a harmonic confinement of frequency ω ⊥ in the cylindrical radial directionρ and by a generic potential V (z) in the cylindrical axial directionz: We introduce the variational field into the fermionic energy functional (2.4) and integrate overx andỹ coordinates. After neglecting the space derivatives ofσ(z) (adiabatic approximation) we obtain the following effective energy functional which depends on two fields: the transverse widthσ(z) and the axial wave functioñ f (z). Note that the variational approach we are using here has been successfully applied in the dimensional reduction from 3D to 1D of the 3D Gross-Pitaevskii equation, which describes trapped Bose-Einstein condensates. In that case the variational approach ends up with the 1D nonpolynomial Schrödinger equation (1D NPSE) for the axial wave function [45]. The 1D NPSE has been extended to investigate the Tonks-Girardeau regime [46,47], the two-component BEC [48], transverse spatial modulations [49] and also axial vorticity [50], Minimizing E 1 with respect tof(z) one finds This is a 1D Schrödinger equation andμ 1 is fixed by the normalization Instead, minimizing E 1 with respect toσ one gets We call (4.4), equipped with (4.6), the 1D nonpolynomial Schrödinger (1D NPS) equation.
The 1D NPS equation can be conveniently written in dimensionless form by scaling z = za ⊥ ,σ(z) = σ(z)a ⊥ ,f (z) = N/2f (z)/ √ a ⊥ , andμ 1 = ω ⊥ (µ 1 + 1/2) as follows: Here we have included a constant term [=1/2, corresponding to the energy of the system in the transverse trap whose effect has been integrated out in (4.7) and (4.8)] in the scaled chemical potential, so that in the N = 0 limit the chemical potential coincides with that of the axial harmonic trap. In the N = 0 limit, from (4.8) we find that σ 2 (z) = 1/2 and the 1/(4σ 2 )+σ 2 terms cancel the constant term in (4.7). In deriving (4.7) we assumed an harmonic confinement in the z direction: V (z) = mλ 2 1 ω 2 ⊥z 2 /2 with λ 1 the anisotropic parameter. By using the 1D NPS equation (4.7) with (4.8) we can study the dimensional crossover from 3D to 1D of the superfluid Fermi gas at unitarity. In general (4.8) must be solved numerically for σ(z) in terms of f (z) and the result when substituted in (4.7) leads to the principal result of this section. No closed-form analytic expression for this equation can be given as in the bosonic case [45], except under special limiting conditions.
The first interesting limit of the formalism is obtained for small number of atoms when the nonlinear term in (4.7) is small, (which corresponds to the weak-coupling limit in the bosonic case,) so that the last term in (4.8) can be neglected and the transverse width is z independent and is given by σ(z) ≡ σ = 1/ √ 2 under the condition N p f 2 (z) ≪ 125 √ 2π/(128(3χ) 3/2 ) = 0.273293... (obviously satisfied for a small number of fermions). The cigar-shaped system is then quasi-1D and governed by − 1 4 In the opposite extreme, for a large number of fermions, N p f 2 (z) ≫ 125 √ 2π/(128(3χ) 3/2 ), the cigar-shaped system is effectively 3D. Under this condition (4.8) can also be solved for σ(z) to yield Consequently, the 1/(4σ 2 (z)) term in (4.7) can be neglected and the remaining terms combined to yield (4.10) The power of the nonlinear term has changed from 7/3 to 9/5 as we pass from quasi-1D regime governed by (4.9) to the effectively 3D regime governed by (4.10). In the quasi-1D regime the nonlinear power 7/3 is the same as that in the original three-dimensional equation (2.5), whereas in the 3D regime it has acquired a different power. Although the quasi-1D equation (4.9) has been used in different studies on a Fermi superfluid [51] and on a degenerate Fermi gas [52], (4.10) for Fermi superfluid is new. We shall see that (4.10) is already valid for a moderate number of Fermi atoms (N > 100) and is of interest in phenomenological applications.
In the 3D regime it is a good approximation to neglect the kinetic energy term in (4.10) and the following analytic expression for density is obtained in the so called TF approximation N p f 2 (z) = 125π where Θ(x) is the Heaviside step function. As we are in the 3D regime it is interesting to compare this result with the following TF approximation made on the full threedimensional equation (2.5) after integrating over the transverse variables The two TF results have the same functional dependence on the variables as well as very similar numerical coefficients, in spite of (4.10) and (2.5) having different powers of density in the nonlinear terms. The quasi-1D equation (4.9) has the same power of density in the nonlinear term as (2.5). Nevertheless, a TF approximation made in (4.9) will generate a density with an entirely different dependence on z.
It is now interesting to study the bulk chemical potential µ 1 (n 1 ) implicit in (4.7) and (4.8) given by where density n 1 (z) = N p f 2 (z). Many physical observables are determined by the bulk chemical potential so obtained. For example, assuming a power-law dependence µ 1 ∼ n s 1 for the bulk chemical potential on density n 1 (polytropic equation of state), the frequency of the lowest axial compressional mode Ω 1 is given by [46] Here we introduce an effective polytropic index s as the logarithmic derivative of the bulk chemical potential µ 1 , that is From (4.13) and (4.14) one finds that in the 1D regime s = 2/3 and Ω = Ω 1 /(λ 1 ω ⊥ ) = 8/3, while in the 3D regime s = 2/5 and Ω = Ω 1 /(λ 1 ω ⊥ ) = 12/5. Note that the result 12/5 is exactly the same one obtains setting γ = 2/3 in the formula (3γ + 2)/(γ + 1) derived by Cozzini and Stringari [53] for the axial breathing mode of 3D cigar-shaped superfluids with a 3D bulk chemical potential that scales as n γ p . In figure 2 (a) we show the power s of the dependence µ 1 ∼ n s 1 as a function of n 1 obtained by numerically solving nonlinear (4.13) and (4.14). In figure 2 (b) we plot the square of collective frequency Ω 2 in the dimensional crossover as a function of the axial density n 1 and As the number of atoms N increases, the frequency of axial compressional mode slightly decreases with a decrease of the polytropic power in µ ∼ n s 1 relation.

Uniform Density
Now we consider the axially uniform case, where V (z) = 0. In this case all the variables attain a constant value independent of z and we remove the z dependence on all variables. Setting n 1 = N p f 2 , from (4.7) and (4.8) we find These equations can be used to derive the axial sound velocity c 1 of the fermionic system, that is obtained with the formula [14] c 1 = n 1 ∂µ 1 ∂n 1 .

Disk-shaped Fermi superfluid: 3D-2D crossover
Let us suppose that the external trapping potential U(r) is given by a generic potential W (ρ) in the cylindrical radial directionρ and by a harmonic confinement of frequency ω z in the cylindrical axial directionz: We introduce the variational field into the fermionic energy functional (2.4) and integrate over thez coordinate. After neglecting the space derivatives ofη(ρ) we obtain the following effective energy functional which depends on two fields: the axial widthη(ρ) and the transverse wave functioñ φ(ρ). Also in this case we observe that the variational approach we are using has been successfully applied in the dimensional reduction from 3D to 2D of the 3D Gross-Pitaevskii equation. The resulting effective equation has been called 2D nonpolynomial Schrödinger equation (2D NPSE) [45].
The 2D NPS equation can be conveniently written in dimensionless form by scaling ρ = ρa z ,η(ρ) = η(ρ)a z ,φ(ρ) = √ N p φ(ρ)/a z , andμ 2 = ω z (µ 2 + 1/4) as follows: with normalization 2π ∞ 0 φ 2 (ρ)ρdρ = 1. In deriving (5.7) and (5.8) we assumed an harmonic confinement in the radialρ direction: W (ρ) = mω 2 zλ 2 2 ρ 2 /2. Again in defining the reduced chemical potential µ 2 we have removed the zero-point energy corresponding to the energy of the axial trap, so that in the N = 0 limit, (5.7) and (5.8) coincides with the corresponding linear harmonic oscillator problem. By using (5.7) with (5.8) we can study the dimensional crossover from 3D to 2D of the superfluid Fermi gas at unitarity. First (5.8) is to be solved numerically for η(ρ) in terms of φ(ρ) and the result when substituted in (5.7) gives the desired result for studying a crossover from 3D to 2D. Closed-form analytic result for these equations is possible only under limiting conditions.
In the 3D regime if we neglect the kinetic energy term in (5.10) and the following analytic expression for density is obtained in the TF approximation after a neglect of the kinetic energy term As we are in the 3D regime it is interesting to compare this result with the following TF approximation made on the full three-dimensional equation (2.5) after integrating over the longitudinal variable The two TF results have the same functional dependence on the variables as well as very similar numerical coefficients, in spite of (5.10) and (2.5) having different powers of density in the nonlinear terms. The quasi-2D equation (5.9) has the same power of density in the nonlinear term as (2.5). Nevertheless, a TF approximation made in (5.9) will generate a density with an entirely different dependence on ρ.
We now consider the bulk chemical potential implicit in (5.7) and (5.8) where density n 2 (ρ) = N p φ 2 (ρ). Again the bulk chemical potential can be considered to possess a power-law dependence on density: µ ∼ n s 2 , where the numerical coefficient s can be extracted from a numerical solution of nonlinear (5.13) and (5.14) using the relation The coefficient s is of interest in the study of physical observables of interest. In the quasi-2D regime s = 2/3 whereas in the 3D regime s = 1/2. However, in the quasi-2D to 3D crossover the µ(n 2 ) dependence is to be calculated numerically using (5.13) and (5.14) and then the polytropic indexs can be calculated as a function of n 2 for the dimensional crossover. In figure 3 (a) we plot the polytropic power s vs. density.

Uniform Density
Now we consider the radially uniform case, where W (ρ) = 0. In this case all the variables attain a constant value independent of ρ and we remove the ρ dependence on all variables. Setting n 2 = N p φ 2 , from (5.7) and (5.8) we find These equations can be used to derive the radial sound velocity c 2 of the fermionic system, that is obtained with the formula [14] c 2 = n 2 ∂µ 2 ∂n 2 . (5.18) In the quasi-2D regime (for small n 2 ), from (5.16) and (5.17) we obtain η 2 = 1/2, for large n 2 in the effective 3D description given by (5.10). The asymptotic result is indistinguishable from the exact result except near very small n 2 . This shows that the effective 3D description of the system is very good. The sound velocity increases with matter density as it should.

Numerical Result
Next we study the effectiveness of the dimensional reduction of the 3D GL equation (2.5) to 1D and 2D forms with a variation of the number of fermions for cigar-and disk-shaped configurations. We numerically solve the full 3D equation as well as various 1D and 2D reduced equations by discretizing them by the semi-implicit Crank-Nicholson algorithm with imaginary time propagation [55][56][57]. For numerical convenience we transform the chemical potential term µψ in the nonlinear equations to a time-dependent term idψ/dt. The space and time steps used in discretization were typically 0.05 and 0.001 respectively.  2) denoted 3D, the full one-dimensional equations (4.7) and (4.8) denoted NPS, the quasi-1D equation (4.9) denoted Q1, and the effectively 3D equation (4.10) denoted E3 in the unitarity limit using parameters ξ = 0.44, λ 1 = 0.1, and λ 2 = 1.

3D-1D crossover
For a cigar-shaped superfluid we solve four sets of equations in the unitarity limit: (a) the 3D equation (3.2), (b) the complete reduced 1D equations (4.7) and (4.8), (c) the approximate quasi-1D equation (4.9), and (d) the approximate effectively 3D equation (4.10). We use the parameters ξ = 0.44 [10], and a highly cigar-shaped trap with λ 1 = 0.1 and λ 2 = 1. In figure 4 we illustrate the results of our calculation by plotting the linear density profile f 2 (z) vs. z of the four sets of calculations for fermion number N = 2, 10 and 100.
From figure 4 we find that the 1D approximate calculations are good approximations to the solution of the 3D equation (3.2). However, some features of the different approximate models are worth commenting. All approximate 1D models lead to a density smaller than that obtained from the 3D equation (3.2). The density obtained from (4.7) and (4.8) provide the best approximation to the exact density for all N. For small N, the quasi-1D model (4.9) provides a better approximation to the exact result than the effectively 3D model (4.10). The opposite happens for large values of N. For an intermediate value of N, the quasi-1D model (4.9) and the effectively 3D model (4.10) could produce similar results.

3D-2D crossover
For a disk-shaped superfluid we again solve for sets of equations: (a) the 3D equation (3.2), (b) the complete reduced 2D equations (5.7) and (5.8), (c) the quasi-2D equation (5.9), and (d) the effectively 3D equation (5.10). We use the parameters ξ = 0.44 in the unitarity limit and λ 1 = 1, and λ 2 = 0.1 for a disk-shaped trap. In figure 5 we present the results of our calculation by plotting the radial density profile φ 2 (ρ) vs. ρ of the four sets of calculations for fermion numbers N = 2, 10, and 1000.
From figure 5 we find that the three 2D approximate models could be good approximations to the solution of the 3D equation (3.2). Again, as in the cigar-shaped superfluid, the full 3D model (3.2) produces the largest density profile with the complete 2D equations (5.7) and (5.8) providing the best approximation to it for all values of N. For small N, the quasi-2D model (5.9) produces better approximation to the exact result than the effectively 3D model (5.10). The opposite happens for large N. For an intermediate N these two latter approximations could produce similar results.
An interesting result from our calculations with cigar and disk-shaped superfluid is that, for N > 100, the quasi-1D model (4.9) and the quasi-2D model (5.9) are poorer approximations than the effectively 3D models (4.10) and (5.10), respectively. For experimental purpose N = 100 represent a small number atoms. Hence for phenomenological applications the effectively 3D models (4.10) and (5.10) with nonlinear terms with power 9/5 and 2 should be used. Note that these powers are different from the power 7/3 in the original 3D equation (3.2). The quasi-1D model (4.9) and quasi-2D model (5.9) for cigar-and disk-shaped superfluids with nonlinear terms of power 7/3 effective for a small number of fermions are useful for academic interest. 2) denoted 3D, the full two-dimensional equations (5.7) and (5.8) denoted NPS, the quasi-2D equation (5.9) denoted Q2, and the effectively 3D equation (5.10) denoted E3 in the unitarity limit using parameters ξ = 0.44, λ 1 = 1, and λ 2 = 0.1.

Conclusion
We have suggested a time-independent Schrödinger equation for a Fermi superfluid at unitarity [(3.2)] by minimizing its energy functional. This equation can also be derived as an Euler-Lagrange equation of an appropriate Lagrangian density. In a cigar-shaped superfluid assuming a Gaussian form for the order parameter, and integrating over the transverse variables we have derived an effective nonlinear nonpolynomial 1D equation for the Fermi superfluid at unitarity [(4.7) and (4.8)]. This complex equation is simplified in the limit of small and large atom numbers when it reduces to a nonlinear equation with power-law nonlinearity. The equation for small atom number N has the same nonlinear structure as the original 3D equation and is called quasi-1D model [(4.9)], whereas the equation for large N has a distinct nonlinearity and is called effective-3D model [(4.10)]. For phenomenological application the effective-3D model seems quite attractive.
In a disk-shaped superfluid assuming a Gaussian form for the order parameter, and integrating over the axial variable we also derived an effective nonlinear nonpolynomial 2D equation for the Fermi superfluid at unitarity [(5.7) and (5.8)]. This complex equation is simplified in the limit of small and large atom numbers when it reduces to a nonlinear equation with power-law nonlinearity. The equation for small atom number N has the same nonlinear structure as the original 3D equation and is called quasi-2D model [(5.9)], whereas the equation for large N has a distinct nonlinearity and is called effective-3D model [(5.10)]. The quasi-2D model has the same nonlinearity as the original 3D equation, whereas the effective-3D model produces a different nonlinearity. Again the effective-3D model is attractive for phenomenological application producing very good results. All the above models have been studied by a numerical solution of the model equations.