Time-Dependent Density Functional Theory for Fermionic Superfluids: from Cold Atomic Gases, to Nuclei and Neutron Stars Crust
pphysica status solidi
Time-Dependent Density FunctionalTheory for Fermionic Superfluids:from Cold Atomic Gases, to Nucleiand Neutron Stars Crust
Aurel Bulgac
Department of Physics, University of Washington, Seattle, WA 98195-1560, USA
Key words:
Density Functional Theory, Superfluids, Cold Gases, Nuclei, Neutron Star Crust ∗ Corresponding author: e-mail [email protected]
In cold atoms and in the crust of neutron stars the pairing gap can reach values comparable with the Fermienergy. While in nuclei the neutron gap is smaller, it is still of the order of a few percent of the Fermi energy.The pairing mechanism in these systems is due to short range attractive interactions between fermions and thesize of the Cooper pair is either comparable to the inter-particle separation or it can be as big as a nucleus,which is still relatively small in size. Such a strong pairing gap is the result of the superposition of a very largenumber of particle-particle configurations, which contribute to the formation of the Copper pairs. These sys-tems have been shown to be the host of a large number of remarkable phenomena, in which the large magnitudeof the pairing gap plays an essential role: quantum shock waves, quantum turbulence, Anderson-Higgs mode,vortex rings, domain walls, soliton vortices, vortex pinning in neutron star crust, unexpected dynamics of frag-mented condensates and role of pairing correlations in collisions on heavy-ions, Larkin-Ovchinnikov phase asan example of a Fermi supersolid, role pairing correlations control the dynamics of fissioning nuclei, self-boundsuperfluid fermion droplets of extremely low densities.
Copyright line will be provided by the publisher
There is a large class of fermionicsuperfluid systems in which the pairing correlations arevery strong and their description and especially the de-scription of their dynamics and interaction with typicallystrong external probes require and extension of the DensityFunctional Theory (DFT) [1] following the `a la Kohn andSham formulation [2], which does not involve non-localpotentials, as in the first extension suggested by Oliveira etal. [3]. If the pairing gap is large, the number of particle-particle configurations defining the anomalous density ismuch larger than the number N of fermions in the system.In the Kohn-Sham version of the DFT the energy densityfunctional (EDF) depends on mainly two types of densi-ties, the number density and the kinetic energy density,which are expressed through N single particle wave func-tions (`a la Hartree-Fock approximation). Since in the lan-guage of the number density alone one cannot distinguishbetween normal and superfluid phases, there is an obvi-ous need to introduce the anomalous density as well, which defines the order parameter, non-vanishing only in the su-perfluid phase. However, the description of fermionic su-perfluids becomes even more demanding, since in practiceone has to study very often superfluids in interaction withstrong time-dependent probes, e.g. when one is stirring afermionic superfluid, when studying non-equilibrium phe-nomena such as quantum turbulence, and when one mayeven observe the evolution of superfluid into a normalphase. Naturally, under such circumstances one needs an-other “order parameter,” capable to disentangle parts of thesystem which evolve in time at various rates, and in thiscase the appearance of current densities in the EDF be-comes unavoidable. The physical systems of interests runthe gamut from cold atom system to nuclear systems andone has to consider multi-component systems for whichthe structure the EDFbecomes quite complex. In nuclei onehas to consider at the same time both the (charged) protonand neutron miscible superfluids and in neutron star crustin addition also include the electron background as well. Copyright line will be provided by the publisher a r X i v : . [ nu c l - t h ] A p r : Figure 1
A rather exotic excitation of a superfluid ex-ist, the Anderson-Bogoliubov-Higgs mode, which corre-sponds to the amplitude oscillation of the order parame-ter [4]. Surprisingly the evolution of the magnitude of thismode in time is unlike the motion of a ball rolling.In the case of cold atoms one is interested lately in misci-ble mixtures of either Fermi-Bose, Fermi-Fermi or Bose-Bose superfluids. In neutron stars various mesons, whichare bosons, are expected to appear at relatively large den-sities, close to the core of the star.
Two prevailing the-oretical models are used to describe the dynamics of su-perfluids. The oldest is the Landau two-fluid hydrodynam-ics [5,6], which at zero-temperature reduces to the hydro-dynamics of a single perfect classical fluid [7,8], namelyof the superfluid component alone. Naturally in such a for-malism Planck’s constant is absent and the two-fluid hy-drodynamics is unable to describe the formation of quan-tized vortices [9], their dynamics, their crossing and re-combination, which is at the heart of the venerable field ofquantum turbulence [10] conjectured by Feyman in 1955.Classical turbulence is due to viscosity, which is absent insuperfluids at zero-temperature. On the same note, thereis no mechanism within the two-fluid hydrodynamics todescribe either the conversion of the superfluid into thenormal component, when a superfluid is stirred vigorously.The quantizations of vortices in the two-fluid hydrodynam-ics has to be enforced by hand [6], in a manner similarto the Bohr 1913 quantization of the hydrogen atom. It isthus impossible to describe the evolution of a superfluid atrest and brought into rotation, when quantized vortices areformed, and later on when they might cross and reconnectas well.An extremely attractive alternative approach was de-veloped by Gross [11] and Pitaevskii [12] to describea weakly interacting Bose at zero-temperature, the cel-ebrated Gross-Pitaevskii equation. This is a non-linearSchr¨odinger equation in which both the density and thesuperfluid order parameter are described by the same com- plex field ψ ( r , t ) : i ¯ h ∂ψ ( r , t ) ∂t = − ¯ h m ∇ ψ ( r , t ) (1) + g | ψ ( r , t ) | ψ ( r , t ) + V ext ( r , t ) ψ ( r , t ) , where g > is the strength of the weak interparticle repul-sion and V ext ( r , t ) is some external potential. While thisfully quantum mechanical formalism is adequate to de-scribe a large range of properties of a weakly repulsiveBose gas, in particular the formation of quantized vortices,it is clearly inadequate to describe the properties of super-fluid liquid He, which is a strongly interacting system.Under the influence of a strong time-dependent externalfield a weakly interacting Bose gas can become normal,at least in some regions of space, but the Gross-Pitaevskiiequation is unable to disentangle between the normal andthe superfluid components.Zaremba, Nikuni, and Griffin[13] presented a nice so-lution to many of these issues in the case of weakly inter-acting bosons. They coupled the quantum Gross-Pitaevskiiequation describing the condensate, with a classical kineticequation for the bosons in the normal component, whichincludes collision between non-condensed bosons andalso collision between the condensed and non-condensedbosons. Landau two fluid hydrodynamics emerges whenthe frequency of the mode satisfies the condition ωτ (cid:28) ,where τ is the relaxation time associated with the collisionsbetween the condensed and non-condensed bosons. Thisnew framework is valid only for weakly interacting bosons.Later on other more sophisticated approaches have beensuggested for the weakly interacting Bose systems [14].Apart from these difficulties discussed above, none ofthese models of superfluids allow for the existence of theAnderson-Bogoliubov-Higgs mode. This mode was no-ticed by P. W. Anderson a long time ago [17,18,19]. Thepotential energy of a system with a complex order param-eter has a shape similar to a Mexican hat, see Fig. 1 takenfrom Ref. [4]). Typically this potential depends only on themagnitude of the complex field but not on its phase, e.g. V ( φ ) = a | φ | + b | φ | , where a < and b > and a mini-mum value at | φ | = − b/ a . The mode | φ | is known inhigh-energy physics as the Higgs boson and it became anessential element of the Standard Model. The existence ofthe Higgs boson leads to masses of quarks, gluons, elec-trons, Z , and W ± bosons, and other elementary particlesand its existence has been determined experimentally [4].If in a homogeneous unitary Fermi Gas (see next sec-tion for its characterization) one would bring very slowlythe pairing gap out of the equilibrium position ∆ in theground state one would be able to observe at least twokinds of rather unexpected excitation modes. The first kindof excitation correspond to small amplitude oscillationsaround the equilibrium, with an unexpected slow algebraicdamping [20,15], when the magnitude of the pairing gap Copyright line will be provided by the publisher ss header will be provided by the publisher 3 −150 −100 −50 0 50 100 150 200 250 30000.20.40.60.81 t ε F | ∆ ( t ) | / ∆ t ε F | ∆ ( t ) | / ∆ (t ε F ) −1/2ab c Figure 2
While exciting the Anderson-Bogoliubov-Higgsmode in an uniform Fermi superfluid the number densityremains constant, while the magnitude of the pairing gapand the phase change in time, but not in space. Unlike a ballrolling down the hill in Fig. 1, the magnitude of the orderparameter will oscillate only between zero and the equi-librium value, < | ∆ ( t ) | < | ∆ | [15]. There is no suchcollective mode in either Landau two-fluid hydrodynamicsor Gross-Pitaevskii equation. The unitary Fermi gas is leftto evolve freely in time with the TDDFT extended to su-perfluid system, after the magnitude pairing gap is broughtvery slowly out of equilibrium to either a relatively smallvalue at t = 0 , panels a and b , or only slightly changedfrom its equilibrium value, panel c . (Here ε F is the Fermienergy and ¯ h = 1 .) In cold Fermi gases one can excitethese modes by means of the Feshbach resonance [16], bypreparing the system with a value of the scattering lengthclose to the BCS limit and suddenly changing it to unitarityat t = 0 . In a uniform Fermi gas by controlling the magni-tude of the scattering length one can control the magnitudeof the pairing gap while maintaining the number densityconstant.behaves as a function of time as | ∆ ( t ) = | ∆ ∞ | + A (cid:112) | ∆ ∞ | t sin( | ∆ ∞ | t + θ ) , (2)where | ∆ ∞ | < | ∆ . The asymptotic state corresponds toa partially fermionic paired state plus quasi-particle ex-citations. If instead the pairing gap is slowly brought toa very small value and after that the system is left toevolve freely, the magnitude of the pairing gap oscillateswith a maximum value smaller than the equilibrium value, < | ∆ ( t ) | < | ∆ | , see Fig. 2. The same results illustratedin Fig. 2 can be obtained by preparing initially the systemwith an interaction strength corresponding to an equilib-rium pairing gap smaller than | ∆ | and suddenly changingthe coupling strengths to a value corresponding to an equi-librium value of the pairing gap | ∆ | . We will illustrate the extension of the Kohn-Sham local density approximation (LDA) of the DFT to super-fluid fermionic systems at first with the case of an unitaryFermi gas, which is both methodologically clear and ofgreat practical value. In a subsection section we will laterbriefly discuss the structure of the EDF for nuclei and neu-tron stars.As in the case of normal electron systems, one of themost frustrating aspects of DFT is the construction of theEDF, which we know it exits according to the Hohenberg-Kohn theorem [21], but for which however we have no welldefined procedure to construct with enough high accuracy.The unitary Fermi gas [22,16], which is a system of spin-up and down fermions, interacting with a zero-range po-tential, characterized by an infinite scattering length and azero range effective become an object of extremely inten-sive study both experimentally and theoretically in the lasttwo decades. George Bertsch noticed [23] that the neutronmatter in the crust of neutron stars is very close to such anidealized system.Neutron interaction in the s -wave is characterized by avery large scattering length a and a relatively small effec-tive range r , with an low-energy s -wave scattering ampli-tude f = 1 − ik − a + r k + . . . , (3)where k is the relative wave vector of two scatteringfermions. The wave function outside the potential is ψ ( r ) = exp( r · k ) + fr exp( ikr ) ≈ − ar + O ( kr ) . (4)The total scattering cross section of two low-energyfermions σ = 4 π | f | → πk reaches the maximum possi-ble value allowed by unitarity if r → and | a | → ∞ . Ifthe scattering length a is positive then a bound state existwith the wave function ψ ( r ) = 1 √ ar exp (cid:16) − ra (cid:17) + O (cid:16) r a (cid:17) . (5)If in a many-fermion system, in which the average inter-particle separation is ∝ n − / , the conditions, r (cid:28) n − / (cid:28) | a | (6)are satisfied such a system is called a unitary Fermi gas(UFG) [16].In 1999 the dilute neutron matter (for which the condi-tion (6) is weakly satisfied) was he closest physical systemto an UFG one could envision, and the calculation of thevalue of the dimensionless Bertsch parameter ξ became atheoretical challenge. If ξ < the system would collapseinto a high density liquid or solid with an average inter-particle separation likely of the order of the range of theinteraction and with properties determined by the particu-lar features of the interaction between fermions. However,if < ξ ≤ the ground state would be that of a gas, and Copyright line will be provided by the publisher : a very unusual gas at that, a superfluid with a pairing gapof the order of the Fermi energy, the largest pairing gapin any physical system, and universal properties largely in-dependent of the details of the interaction. One can eas-ily establish using dimensional arguments alone that theground state energy of a uniform UFG should be given bya function depending only on the volume V , particle num-ber N , ¯ h , and the fermion mass M , in the unitary limitwhen r → and | a | → ∞ : E gs ( N, V, ¯ h, m, a, r ) → ξN ε F , (7)where ε F = ¯ h k F m is the Fermi energy of a uniform non-interacting Fermi gas with the Fermi wave-vector k F = (cid:0) π NV (cid:1) / , and ξ is the dimensionless Bertsch parameter. ξ ≤ , as the scattering length a can become infinite onlyin the case of attractive interactions. Apart from the univer-sal value of the Bertsch parameter ξ , the EDF for the UFGshould have the same structure as for a non-interactingFermi gas and the only parameter specifying the natureof the fermions is their mass. Both theoretical QuantumMonte Carlo (QMC) ξ = 0 . [24] and experimental ξ = 0 . [25] values of the Bertsch parameter are nowin very good agreement with each other.A fermionic superfluid under the influence of a time-dependent external field might become normal in somespatial regions but not in others, and the number densityalone cannot discriminate between different phases. Thecase of the UFG is particularly attractive from the point ofview of a DFT aficionado, as only rather general require-ments suffice to narrow down the structure of the EDF.Dimensional arguments, rotation, translational, and par-ity symmetries, gauge symmetry related to the transforma-tions of the complex order parameter, Galilean covariance,and renormalization and regularization of the EDF com-bined with the so called adiabatic local density approxima-tion (ALDA) [26] extended to the case of superfluid sys-tems restrict the form of the (unregulated) EDF for UFG toa rather simple form [27,28,29], namely: E ( r , t ) = ¯ h m (cid:20) ατ ( r , t ) + β
35 (3 π ) / n / ( r , t ) (8) + γ | ν ( r , t ) | n / ( r , t ) + (1 − α ) j ( r , t ) n ( r , t ) (cid:21) + V ext ( r , t ) n ( r , t ) , where α , β , and γ are dimensionless constants. n ( r , t ) , τ ( r , t ) , ν ( r , t ) , and j ( r , t ) are the unregulated number, ki-netic, anomalous densities and current densities of a fullyunpolarized UFG (equal number of spin-up and spin-downparticles) and expressed through the Bogoliubov quasi-particle amplitudes u n ( r , t ) and v n ( r , t ) [30]. V ext ( r t ) isan arbitrary external field with which one might probeor excite the system. We refer to this form of DFT asthe (Time-Dependent) Superfluid Local Density Approx-imation ((TD)SLDA), which is a natural extension of theLDA for normal systems of Kohn and Sham formulation of the DFT [2] to superfluid systems. The emerging TD-SLDA equations have the expected form, identical to theBogoliubov-de Gennes equations, i ¯ h ∂∂t (cid:32) u n v n (cid:33) = (cid:32) h − µ ∆∆ ∗ − h ∗ + µ (cid:33) (cid:32) u n v n (cid:33) , (9)where h = ∂ E ∂n , ∆ = ∂ E ∂ν ∗ , and µ are the single-particleHamiltonian, the pairing field, and the chemical potential.Both kinetic energy and anomalous densities diverge as afunction of the ultraviolet cutoff in a very similar mannerin the case of a zero-range interaction [28]. In particularthe anomalous density matrix has the behavior ν ( r , r , t ) ∝ | r − r | , if | r − r | → . (10)This divergence is the same as the divergence of the scat-tering wave in Eqs. (4) or of the bound state wave function(5) when r → r → . This divergence is “real” and nota deficiency of the theoretical formalism. The divergenceof the anomalous density matrix reflects nothing else butthe increase of the wave function of the Cooper pair whenthe separation between the fermions approaches the radiusof the interaction. One can relate the anomalous densitymatrix ν ( r , r , t ) with the Cooper pair wave function.The mathematical difficulty in extending DFT in ´ala Kohn and Sham manner to superfluid system ariseswhen one attempts to reach the limit n / r → , whenthe average interparticle separation is much larger thanthe radius of the interaction and avoid the appearanceof infinities in the calculations of various densities andof the pairing gap. In meanfield approximation the pair-ing field ∆ ( r , r , t ) = − V ( r − r ) ν ( r , r , t ) , where V ( r − r ) is the fermion-fermion interaction responsiblefor the pairing correlations, is formally non-local. In thelimit n / r → densities should be calculated with acutoff [31,28]: n ( r , t ) = 2 (cid:88) In the ground state the single-particle Hamiltonian of anunpolarized UFG has the structure h = − α ¯ h m ∆ + U ( r ) . (15)One needs to introduce the momentum dependent wave-vectors and the renormalized coupling constant α ¯ h k ( r )2 m + U ( r ) − µ = 0 , (16) α ¯ h k c ( r )2 m + U ( r ) − µ = E cut , (17) g eff ( r ) = n / ( r ) γ (18) − mk c ( r )2 π ¯ h α (cid:20) − k ( r )2 k c ( r ) ln k c ( r ) + k ( r ) k c ( r ) − k ( r ) (cid:21) , in roder to derive the renormalized form of the pairing gap ∆ ( r ) = − g eff ( r ) ν ( r ) . (19)One can then show that the combination α ¯ h m τ ( r ) − ∆ ( r ) ν ∗ ( r ) (20)does not diverge when E cut → ∞ . While the wave vector k ( r ) can become imaginary in the classically forbiddenregions of space, the effective coupling consnat g eff ( r ) re-mains real. If k c ( r ) becomes imaginary in any spatial re-gion the recipe is that the last term on the rhs side of Eq.(18) should be dropped [31,28].The EDF for the UFG Eq. (8) depends on three di-mensionless parameters α, β and γ , for both superfluid andnormal phases, which can be extracted from values of theBertsch parameter ξ , the pairing gap, and the momentumdependence of the quasi-particle excitations obtained in theQMC for the uniform UFG [32,33,24,34,35,36], whichagree well with extracted experimental values [37,38,25].In the case of polarized UFG one needs two numberdensities n a,b ( r , t ) and from dimensional arguments aloneit follows that the energy density of a uniform polarizedUFG is: E ( n a , n b ) = 35 ¯ h m (6 π ) / (cid:20) n a g (cid:18) n b n a (cid:19)(cid:21) / (21)with g (1) = (2 ξ ) / , see Fig. 3. At this point one is inthe position to validate the accuracy of the suggested formof the EDF for the UFG (8) by placing various number offermions with spin-up and spin-down in external fields. AUFG in a harmonic trap is particularly interesting, sincefrom theory we know that its properties are determined bythe only energy scale in the system ¯ hω [43,44,45], theirproperties can be calculated numerically with controlled . . . . . . . . . x = n b /n a g ( x ) Figure 3 The function g ( x ) can be extracted from QMCand other theoretical calculations [36,39,40], blue pointswith error bars for the normal state of UFG and the blackpoint with error bars at x = 1 for the unpolarized super-fluid UFG) and also directly from experiment [41] (greencrosses). The black dashed line represents the mixed phasediscussed in Refs. [36,39,40]. The red line is for theLarkin-Ovchinnikov (LO) phase: the solid line is for thepure state and the dashed line is for the mixed state withthe superfluid unpolarized state[42]. At the left end of thered solid line there is a second order phase transition be-tween the normal and the LO phase, while at the right endof the red solid line there is first order transition betweenthe pure LO phase and the mixed phase. Unlike in the caseof weakly BCS superfluids, the LO phase exits at unitarityin a very wide window of polarization x = n b /n a .accuracy within QMC and unlike the the homogeneoussystem density gradients are important. The EDF for theUFG has so far been constructed for uniform systems andit is not obvious that such an EDF would perform well forinhomogeneous systems, where density gradients are sig-nificant. Fortunately there exist a sufficiently large num-ber of QMC calculations of unpolarized and polarized fi-nite systems of fermions at unitarity in a external harmonictrap, both in the superfluid and normal phases. The EDFfor the polarized UFG has a little more complex structure,as in this case there are two independent number densities,for spin-up and spin-down fermions. Similar arguments asthose used for the derivation of the EDF for an unpolar-ized UFG can be used for its derivation [46,42,47,28].In case of inhomogeneous systems a new gradient termmight need to be included in Eq. (8), namely one of theform ∝ (cid:12)(cid:12)(cid:12) ∇ (cid:112) n ( r , t ) (cid:12)(cid:12)(cid:12) , but comparison with QMC calcu-lations of inhomogeneous systems indicate that it can beneglected [27]. Polarized fermions at unitarity can be in ei-ther superfluid or normal phase, depending on the degreeof polarization N ↓ /N ↑ [46,42,47,28] and the same EDFcan be used to describe with surprisingly high accuracyall such systems. When comparing the results of the QMCcalculations performed for trapped unitary fermions [48,49] with the predictions of the SLDA one finds almost per-fect agreement [27,28]. The differences between the QMCcalculations and the SLDA predictions are almost alwayswithin the statistical errors of the QMC calculations, withinat most 1-3%. There are two exceptions. The QMC and Copyright line will be provided by the publisher : SLDA calculations for the ( N ↑ , N ↓ ) = (15 , have anerror of 9.5%, attributed to the inaccuracies in the QMCresults for this particular system. The largest disagreementis for the two-fermion system ( N ↑ , N ↓ ) = (1 , , about15%. The rest of the differences between the QMC re-sults and the SLDA predictions are at the level of . . . %, which is also the level of accuracy of the QMC calcu-lations. Surprisingly, the odd-even staggering, namely theenergy differences between the ground state energies ofeven N ↑ = N ↓ and odd systems N ↑ − N ↓ = 1 are withinstatistical errors as well. One should keep in mind theground state energies of harmonically unpolarized trappedsystems of unitary fermions scale as √ ξE NI ∝ ¯ hωN / , where E NI is the energy of N non-interacting fermions ina harmonic trap. The odd-even energy differences dependrelatively weakly on the total particle number [48,49] andthe relatively small odd-even energies differences are re-produced within SLDA with surprising accuracy as well.There are only a few exact solutions of the time-dependent Schr¨odinger equation for interacting many-fermion systems. The linear response theory predictsdamped harmonic oscillations with a frequency | ∆ | [17,20], while Eq. (2) emerges when nonlinearities are takeninto account. Equally unexpected is the behavior of thepairing gap when the initial disturbance is large, panels a and b in Fig. 2, when one would naively expect thatthe pairing gap would oscillate somehow around the equi-librium value | ∆ | . While these modes have been put inevidence in a meanfield-like framework and one mightsuspect that they are artifacts of possible approximations,the same behavior was demonstrated to appear in an ex-actly soluble time-dependent realistic many-body model ofsuperconductivity, the Richardson or Gaudin model [50].We will illustrate below however the power of TDSLDAin confronting actual non-equilibrium phenomena in realexperiments. Nuclear systems are signifi-cantly more complex that the UFG. While we know fordecades quite a bit about the nuclear EDF (NEDF), the ex-act from and accuracy of mostly empirical NEDF is stillinsufficient for many applications, like predicting the ori-gin of elements in Universe. In case of nuclear systemsthere two type of fermions, protons and neutrons, and thespin-orbit interaction is very strong, and has the oppositesig when compared to atomic systems. Pairing correlationsare relatively strong as well, though not as strong as in thecase of UFG. However, the dilute neutron matter, foundin the crust of neutron stars, has quite a lot of similaritieswith the UFG [23,22]. QMC calculations for nuclear sys-tems a significantly more complex than for electronic orcold atom systems, and so far only the pure neutron mat-ter equation of state as function of density is know withreasonable accuracy from it ab initio calculations. The in-teraction between nucleons (neutrons and protons) is verycomplex and in nuclei not only two-body interactions areimportant, but also three-body (and even four-body) inter- actions play a crucial role. As a result the form of NEDFis typically obtained with significant phenomenological in-put [51,52]. The accuracy of phenomenological NEDFs inpredicting the binding energies of about 2300 known nucleiis nowadays at the sub-percent level, but still not accurateenough for many applications.In the resulting evolution equations we suppressed forthe sake of simplicity the space and time coordinates ( r , t ) .The ensuing equations represent an infinite set of couplednonlinear time-dependent 3D partial differential equationsfor the quasi-particle wave functions, i ¯ h ˙ u k ↑ ˙ u k ↓ ˙ v k ↑ ˙ v k ↓ = h ↑↑ h ↑↓ ∆h ↓↑ h ↓↓ − ∆ − ∆ ∗ − h ∗↑↑ − h ∗↑↓ ∆ ∗ − h ∗↓↑ − h ∗↓↓ u k ↑ u k ↓ v k ↑ v k ↓ . (22)Here both the local mean field h σσ (cid:48) and pairing field ∆ depend on the various single-particle densities. The index k labels each quasi-particle wave function, and is bothdiscrete and continuous. This index must also run overisospin, so that there are similar sets of equations for bothprotons and neutrons, which naturally are coupled. Wehave explicitly included the spin indices ( σ, σ (cid:48) ) ∈ {↑ , ↓} ,allowing for mixing between the spin-up and spin-downstates by the spin-orbit interaction, thus capturing the ef-fects of proton-proton and neutron-neutron pairing. The TDSLDAequations are discretized in space and time. The systemof interest is placed on a spatial lattice with N x N y N z lattice points, for a chosen lattice constant l , which deter-mines the momentum cutoff ¯ hπ/l [53]. When these sim-ulation box parameters are chosen appropriately one canensure that with further discretization the corrections are(exponentially) small [53]. The equations are propagatedin time using a 6th order multi-step predictor-corrector-modifier, which typically requires an application of thequasi-particle Hamiltonian only twice per time-step. No-tice that in the case of the popular Runge-Kutta 4th ordermethod one would need to apply the Hamiltonian fourtimes per time step. Unitarity and accuracy of various con-served quantities during the time evolution are satisfiedwith high accuracy for up to millions of time steps. Thenumber of coupled complex 3D time-dependent nonlin-ear partial differential equations which has to be evolvedin time is up to O (10 ) , see supplemental material inRefs. [28,54,55]. The numerical solution of these equa-tions requires the use of the leading edge supercomputersand it ranks among some of the largest direct numericalsimulations attempted so far. This section will briefly describe quali-tative aspects of static and dynamic aspects of fermionicsuperfluids in strongly interacting cold atomic gases, nu-clei and neutron star crust. Copyright line will be provided by the publisher ss header will be provided by the publisher 7 One of the first ap-plications of the SLDA for UFG was to determine thestructure of a quantized vortex [56,57]. It was shown thatthe actual number density profile of a quantized vortex issomewhere between that in a BCS superfluid and in a BECsuperfluid. While in a BCS superfluid the density in thecore is practically the same as the density far away fromthe vortex core and only the order parameter vanishes atthe vortex core, in a BEC quantized vortex both the orderparameter and the number density vanish at the core, inthe case of dilute superfluids. In the case of the UFG whilethe order parameter vanishes at the core, the vortex coreis only partially filled with fermions, with a density onlyabout half the value of the asymptotic value. This densitydepletion of the vortex core was used to visualize in ex-periments the Abrikosov vortex lattice formed in a rotatingUFG [58], which was the deciding experimental argumentin demonstrating that UFG is indeed a superfluid. Once the EDFof the UFG was established and validated, it was used toin order to establish if a UFG can sustain an inhomoge-neous state of the order parameter [59,60]. In Ref. [42] theselfconsistent SLDA equations for a system in which theorder parameter was allowed to oscillate in the z -direction,while in the x and y directions the properties of the systemremains homogeneous, were solved Unlike in the case ofa weakly coupled BCS superfluid, where the LO phase ex-ist only for a very narrow window of spin polarization, theUFG was shown to sustain such a phase in a surprisinglywide range, see Fig. 3. Only the amplitude of the oscilla-tions of the number density of the minority component aresignificant. This LO phase has not been observed yet incold fermi atom systems. The most fascinating applications of theSLDA is to non-equilibrium phenomena. In the first simu-lation of a fermionic superfluid in real time [61] we placeda UFG in a container resembling a soda can, homogeneousand with periodic boundary conditions along the longitu-dinal direction, along with other geometries as well. Weinserted a “straw” and started stirring the fluid with vari-ous constant angular velocities and linear velocities smallerand larger than the Landau critical velocity. A UFG, apartfrom being characterized by a very large pairing gap, whichin appropriate units is even bigger than in high T c super-conductor [62]. An UFG one has perhaps the largest Lan-dau critical velocity (in appropriate units) of any super-fluid [63,64,65]. Many of the results obtained are avail-able in the form of videos online, for various geometries,various ways to stir the superfluid, and a range of stirringvelocities ranging from very slow to well above the Landaucritical velocity [66]. If one introduces an object and movesit relatively slowly through the superfluid and eventuallyextract that object slowly also the UFG returns practicallyto its initial state, as one would have naturally expected Figure 4 The spatial profile of the pairing gap ∆ ( z ) andof the number densities n a,b ( z ) of the majority (dotted)and minority species (solid) in the region where a pure LOphase exist, see solid red line in Fig. 3. For polarizationclose to the left end of the solid red curve in Fig. 3, thespatial shape of the order parameter is very similar to asine-function and the amplitude of the oscillation of theorder parameter is small (blue curves). Close to the rightend of the solid red curve in Fig. 3 the spatial shape of theorder parameter starts resembling a domain wall of finitewidth. For each polarization the optimal period of the LOphase determined.for an adiabatic evolution. We have also noticed that some-times we can bring the atomic cloud into rotation even ifthe linear speed of the “straw” exceeds Landau critical ve-locity and the cloud remained superfluid. This is possiblesince UFG is gas, during rotation accumulates along thewall, the density and therefore the local Landau critical ve-locity increases. In the same work [61] we have demon-strated in a real-time treatment for the first time that in afermionic superfluid quantized vortices can cross and re-combine, exactly as Feynman [9] envisioned and suggestedthat quantum turbulence emerges, see also subsection 4.6. Thomas and collabora-tors demonstrated that quantum shock waves can be ex-cited in a cold atomic fermionic cloud [ ? ]. Thee shockwave front was directly visualized as a clear number den-sity discontinuity. In a TDSLDA simulation of this exper-iment [67] is was shown that the character of these shockwaves is controlled by the dispersive effects and not dis-sipative effects as was assumed in the analysis of the ex-periment [68]. At the front the shock wave both the num-ber density and the velocity field are discontinuous, whencoarse grained appropriately. At the front of a shock wavethe matter flows in opposite directions. But in addition withquantum shock waves also domain walls are formed, whichpropagate thought the cloud, see Ref. [67] and Fig. 6. A do-main wall is the region where the phase of the order param-eter changes from − π to π and is topological in character,similar to domain walls in magnetism. At a domain wall Copyright line will be provided by the publisher : Figure 5 The first two rows show the magnitude and thephase of the pairing gap ∆ ( x, y ) at various times duringstirring. The position of the “straw” (which parallel at alltimes to the axis of the “can”) is visible as small notchin | ∆ ( x, y ) | at the edge of the “can” and its linear speedis always smaller than Landau critical velocity. The vor-tex Abrikosov lattice emerges quite rapidly and the UFGis also partially depleted along the rotation axis. If the lin-ear speed of the “straw” exceeds Landu critical velocity,the order parameter vanishes quite rapidly in time, see thethird row.number density has a significant depletion, as in a vortexcore. In a surpris-ing experiment performed in 2013 [69] it was announcedthat new kind of many-body excitation was observed inan experiment performed on an UFG. Half of a very elon-gated superfluid cloud was illuminated with a laser, whichresulted in a difference in the phase of the order param-eter between the two halves. Theoretically is was knownfor quite some time that in this case such a planar solitonpropagates with a know speed, and in 3D can be unstabledue to the snake-instability. What the MIT experiment es-tablished was this the excitations they created was movingwith a speed two orders of magnitude slower than the the-ory predicted and the dubbed this mode a “heavy soliton.”In a very careful TDSLDA analysis [70] we establishedhowever, that under the conditions described in the originalpaper [69] most likely the authors observed a vortex ring. Figure 6 Three sequential frames of the number density.The speed of the front of the shock wave in experimentand simulations agree, though not perfectly. The quan-tum shock wave form after two clouds collide and startexpanding, but in their wake one can see several densitydepletions, which are due to the excitation of the domainwalls [67].A planar soliton in an inhomogeneous cloud quite rapidlyevolves into a vortex ring in simulations and it starts prop-agating back and forth in an almost perfect harmonic mo-tion. However, while moving in one direction the vortexring is large, and while is moving in the opposite directionthe vortex ring shrinks, and then the motion is repeated,see Fig.7. In subsequent more detailed experiments [71,72] the MIT group confirmed the transformation of the pla-nar soliton, into a vortex ring. However, since their traplack azimuthal symmetry, which was not adequately estab-lished in the original paper, a vortex ring in an axially non-symmetric trap rather quickly touches the walls and turnseither into a single or two line vortices, dubbed howeveror solitonic vortices, if vortex line is perpendicular to thepropagation direction. This is fully in agreement with ear-lier simulations of an UFG [67], see supplemental mate-rial [66], and subsequent analysis of the new MIT experi-ment [73,74]. At the hearty of quantumturbulence lie the formation of a tangle many quantizedvortices. In experiments with liquid helium such tangleshave been created in laboratories for decades [10]. In thecase of cold fermionic gases vortices can be easily gen-erated with rotating laser beams [58] but also by simplyilluminating part of a cloud with a laser [69]. By com-bining these two methods, which clearly are not the onlypossibility, one can generate a tangle of quantized vorticesin both unpolarized [70,74] and Fig. 8, and in polarizedsystems [75]. The great advantages of cold atom systemover liquid helium are multiple: i) one can control easilymany parameters of the system, including the interactionstrength, and create a wide array of external probes [16]; Copyright line will be provided by the publisher ss header will be provided by the publisher 9 Figure 7 After a planar soliton is created by illuminat-ing with a laser half of an elongated fermionic superfluida vortex ring is formed. This vortex ring propagates in onedirection with a larger radius and in the opposite directionwith a smaller radius [70].ii) one can describe with very good accuracy theoreticallyboth static properties and particularly the non-equilibriumdynamics of such systems with very good control of thetheoretical ingredients. being able to confront theory andexperiment in great detail is a great advantage of cold atomsystems over studies performed in liquid helium, whereonly phenomenological models exist. Our main goal in developingSLDA was to describe non-equilibrium dynamics of nu-clear systems, but in order to verify and validate the the-oretical framework and the complex numerical implemen-tation we made quite a long detour studying cold atomssystems, for which theoretical tools are both simpler andin better control, and also there is a great wealth of ex-perimental data which can be confronted with theory pre-dictions and postdictions ! One of the oldest problems ofstrongly interacting quantum many-body systems i s nu-clear fission. Unlike superconductivity for example, whichrequired less than five decades form the initial experimen-tal observation [76] in 1911, until a microscopic theorywas put forward in 1957 [77], nuclear fission observed in1939 [78] will turn 80 years old in 2019 and likely will stillhave no adequate microscopic description. There are manyreasons why this is the case and here are some of them:i) the nuclear interactions are extremely complex and theyare not yet accurately known, unlike Coulomb interactionbetween electrons and nuclei; ii) nuclei are finite systemsand at the same time they have too many particles, cor-relations are very stroong; iii) when a heavy nucleus fis-sions the number of final channels is in the hundreds, corre-sponding to various possible splittings of protons and neu-trons between the fission fragments, the fission fragmentsemerge excited with various quantum numbers; iv) apart v/v F l og P D F ( v ) v ⊥ v k Figure 8 Three consecutive frames illustrating the cross-ing and recombination of quantized vortices in a coldatomic Fermi superfluid are shown in the left column.In the right column the longitudinal (red) and transversal(blue) momentum distributions compared with a thermalmomentum distribution in logarithmic scale. density (fm −3 ) pairing gap (MeV) −π −π/2 0 π/2 π pairing phase Pu fission with the normal pairing gap Figure 9 The three column show the evolution of the neu-tron/proton number density, magnitude and phase of thepairing respectively (upper/lower half of each frame), fromthe top of the outer fission barrier, past scission, until fis-sion fragments are separated.from fission fragments typically quite a number of neu-trons are emitted, during fission, immediately after fission,and much later, along with gamma rays and beta decays. Copyright line will be provided by the publisher One cannot declare “victory” until theory is able to predictwith reasonable confidence most of these fission fragmentproperties, which are critical for many applications as well,and also for clarifying many fundamental questions, suchas the origin of elements in the Universe, and the structureand evolution of stars.Using TDSLDA and an NEDF of reasonable quality in2016 we were able for the first time for describe the evo-lution of a fissioning nucleus from the outer fission barrierto scission, until fission fragments were separated [79] andFig.9 and recently also in Ref. [55]. The fission fragmentsemerge with properties similar to those determined exper-imentally, while the fission dynamics appears to be quitecomplex, with various shape and pairing modes being ex-cited during the evolution. The time scales of the evolutionare found to be much slower than previously expected andthe role of the collective inertia in the dynamics is found tobe negligible. Even though in this first study of its kind wedid not obtain a perfect agreement with experiment, our re-sults clearly demonstrate that rather complex calculationsof the real-time fission dynamics without any restrictionsare feasible and further improvements in the quality of theNEDF, and especially in its dynamics properties, can leadto a theoretical microscopic framework with great predic-tive power. TDSLDA will offer insights into nuclear pro-cesses which are either very difficult or even impossible toobtain in the laboratory. Figure 10 A nucleus moving at constant velocity repelsthe quantized vortex line. The little green arrows displaythe relative magnitude and direction of the repulsive forcebetween the nucleus and the small linear element of thevortex line, while the long green line is the force exertedby the vortex on the nucleus.. A ratherold puzzle in nuclear astrophysics is to explain why neu-tron stars/pulsars experience glitches. Neutron stars are the most compact in the Universe (not counting black hole,which are not made of typical matter) with a density reach-ing ≈ g/cm and they rotate extremely fast, with pe-riods as short as milliseconds. As they evolve they radi-ate energy and their rotation period increases very slowly.However, at practically random times they experience sud-den accelerations and their periods increase very rapidly,practically instanteneously. Anderson and Itoh [80] sug-gested that since neutron stars are mainly superfluid neu-tron matter, there are many quantized vortices, which likelyare pinned in the so called neutron star crust. Neutron starsare not composed entirely of neutrons, which are unstableagainst beta-decay, but also of a relatively small fraction ofprotons and electrons. Electrons are relativistic, very delo-calized, and likely in a normal phase. But the neutrons andprotons lead to the formation of the so called pasta phase,a series of nuclei ranging in shape from small spheres(meat balls), long sticks (spaghetti) , plates (lasagna), and“swiss cheese-like” matter, with holes like bubbles andtubes. These “nuclei” from a Coulomb lattice in a super-fluid neutron superfluid. The neutron superfluid under ro-tation is threaded with quantized vortices, which might at-tach to nuclei or interstitially. During the neutron star rota-tion these vortices can unpin and find a new position, whichleads to a change in the star moment of inertia and there-fore its period, or glitches, in analogy with flux creeping inhard superconducting magnets according to Anderson andItoh [80].For decades nuclear physicists and nuclear astrophysi-cists attempted to figure out whether a quantized vortex ispinned on a nucleus or whether it is repelled by a nucleus inthe neutron star crust, and the results have been all over themap, from attraction to repulsion, with no clear indicationeven on the magnitude of this interaction. In 2013 we sug-gest a new approach to this problem, which could simplyindicate whether the interaction is either repulsive or at-tractive. A nucleus scattering on a quantized vortex wouldbe deflected either to right or to the left depending on thecharacter of the force [81]. The size and sign of the nucleusvortex interaction arises in any theoretical approach as adifference between two very large quantiles, which can-not be computed with enough accuracy. In 2016 we haveimplemented this recipe [82], see also Fig. 10 and we wereable to extract even more finer details of the nucleus-vortexinteraction, namely how each small linear element of thevortex line interacts with the nucleus. During their inter-action the vortex line is deformed and stretched. We haveestablished not only that the interaction between a nucleusand a vortex is repulsive, for all possible densities on theneutron matter, but we have extracted with quite high ac-curacy the magnitude of this interaction. I will address only a single aspect of this rather widefield, which likely is of interest outside nuclear physics, thephase locking between two colliding superfluid droplets.Since the gauge symmetry is spontaneously broken in su- Copyright line will be provided by the publisher ss header will be provided by the publisher 11 Figure 11 The evolution of the phase of the pairingfield (time runs top to bottom) in the head-on collision of Sn+ Sn. The right and the left columns correspondto a realistic or artificially increased pairing field strengthrespectively. The upper and lower half of each frame cor-responds to an initial phase difference between the two ini-tial pairing condensates of 0 and π respectively. The phaselocking of the pairing field is clearly manifest after fusionin the left column, but absent in the right column.perfluids, it is reasonable to wonder under what conditionsthe relative phase of two superfluids is physically relevant.The Josephson effect [83,84], experiments with cold Boseor Fermi atoms [85,69,71,72,86,73], and the superfluidfragments emerging from nuclear fission [87], are just afew examples where that is the case.Recently Magierski, et al. [88] reported on a surprisingand strong dependence of the properties of the emergingfragments on the relative phase of the pairing condensatesin the initial colliding nuclei. Related theoretical issues insuperfluid helium have been discussed in literature by An-derson [89]. In Ref. [90] we have demonstrated that this istypical behavior of superfluids in the weak coupling limit.With increasing interaction strength however, the initiallyindependent phases of the two order parameters in the col-liding partners quickly become phase locked, as the strongcoupling favors an overall phase rigidity of the entire con-densate, and upon their separation the emerging superfluidfragments become entangled, see Fig. 11. One of themain difficulties with TDDFT is that in a time-dependentframework is likely to always lead to the same final state,even if one were to consider a range of different initialconditions[55]. As I discussed briefly above, at least inthe case of nuclear fission and nuclear reactions in generalthat is definitely not the case. In nuclear fission one has awide distribution of fission fragments masses and charges.Can one find within TDDFT a sensible solution to thisproblem?The prevalent approach to perform real time evolutionof many-nucleon systems was developed in the 1970’s us-ing path-integral techniques; see [91] for a review. Startingwith a system described by a many-body Hamiltonian ˆ H ,one performs a Hubbard-Stratonovich transformation onthe many-nucleon evolution operator by introducing aux-iliary one-body fields σ ( t ) , formally, ˆ U ( t i , t f ) = exp (cid:20) − i ¯ h ˆ H ( t f − t i ) (cid:21) = (cid:90) D [ σ ( t )] W [ σ ( t )] exp (cid:18) − i ¯ h (cid:90) t f t i ˆ h [ σ ( t )] (cid:19) , where D [ σ ( t )] is an appropriate measure depending on allauxiliary fields, W [ σ ( t )] is a Gaussian weight and ˆ h [ σ ( t ] isa one-body Hamiltonian built with the auxiliary one-bodyfields σ ( t ) .Using the stationary phase approximation, one selects asingle mean field trajectory ¯ σ ( t ) , which one may simulatewith the TDDFT trajectory. If the initial state is a (gener-alized) Slater determinant, the final state is also a (gener-alized) Slater determinant under the evolution of this sta-tionary phase mean field trajectory. After a trivial changeof integration variables σ ( t ) = ¯ σ ( t ) + η ( t ) , the true many-body wave function can be put into the form Ψ ( t ) = (23) (cid:90) D [ η ( t )] ¯ W [ η ( t )] exp (cid:18) − i ¯ h ˆ h [¯ σ ( t ) + η ( t )] (cid:19) Ψ (0) , where Ψ (0) is the initial wave function, and η ( t ) are fluc-tuations around the stationary phase trajectory ¯ σ ( t ) at time t . In Eq. (23) the weight functions are Gaussian-shaped.Thus, the true many-nucleon wave function is now a time-dependent linear superposition of many time-dependent(generalized) Slater determinants. In this respect the truemany-nucleon wave function has a similar mathematicalstructure as the wave function in the time-dependent gener-ator coordinate method (TDGCM) introduced by Wheeler et al. [92,93]. One cannot but see the analogy in treat-ing fluctuations around the mean field trajectory with theclassical Langevin description of nuclear collective motionas well [94]. The representation (23) (which is an exactone) of the many-body wave function has the great advan-tage that each trajectory is independent of all the others.One particular aspect of this general structure of the many-nucleon wave function is the nature of the initial wave Copyright line will be provided by the publisher function. One choice is a single (generalized) Slater deter-minant and another is a superposition of many such (gen-eralized) Slater determinants.In our fission studies we observed that one-body dissi-pation (i.e. the collisions of fermions with the moving sur-face of the nucleus) or Landau type of dissipation leads to avery quick dissipation of the f=collective flow energy intointrinsic/thermal of the fissioning nucleus [55]. The intrin-sic motion of the descending nucleus from the outer saddletowards the scission configuration is similar to the down-ward motion of a heavy railway car on a very steep hill withits wheels blocked. The wheels dot not rotate but slip andbecome extremely hot, since almost the entire gravitationalpotential energy of the railway car at the top of the hill isconverted into heat and very little of it is converted intocollective kinetic energy. In this case, the railway car ve-locity is equal to the terminal velocity. An object attains aterminal velocity when the conservative force is balancedby the friction force, the acceleration of the system van-ishes and the inertia plays no role in its dynamics. The mo-tion of the railway car is strongly non-adiabatic. Similarlyfor a nucleus the flow of energy from the collective/shapedegrees of freedom is controlled by the entropy of the in-trinsic degrees of freedom. This energy transfer is to someextend stochastic, as is brownian motion, and can be sim-ulated with TDDFT evolution equations augmented to in-corporate dissipation and fluctuations we introduce havethe form i ¯ h ˙ ψ k ( r , t ) = h [ n ] ψ k ( r , t ) + γ [ n ] ˙ n ( r , t ) ψ k ( r , t ) (24) − 12 [ u ( r , t ) · ˆ p + ˆ p · u ( r , t )] ψ k ( r , t ) + u ( r , t ) ψ k ( r , t ) , where ˆ p = − i ¯ h ∇ (not to be confused with p ( r , t ) ), theindex k runs over the neutron and proton quasi-particlestates and where ψ k ( r , t ) are 4-component quasi-particlewave functions and h [ n ] is a × partial differential oper-ator [95,87]. The fields u ( r , t ) and u ( r , t ) generate bothrotational and irrotational dynamics and the term propor-tional to γ is a quantum frictional term. In the presence ofthis additional term alone ˙ E tot ≤ , as in the case of thepresence of a classical friction term. This is only a phe-nomenological solution, which likely will lead to reason-able results and might serve as inspiration. I presented an extension of the DFTto superfluid systems using a local pairing field and a fur-ther extension of this framework to time-dependent phe-nomena using the adiabatic approximation. This static ex-tension SLDA can be correlated with ab initio calculations,it is strongly constrained by a number of theoretical ar-guments, and at least in the case of cold atoms it appearsto have a quite good accuracy. The applications to nuclearphenomena is based to a large extent to a phenomenologi-cal energy density functional, which has a sub-percent ac-curacy for a large number of nuclei and their static proper-ties. The time-dependent extension, which in addition isrequired to satisfy local Galilean covariance, appears to provide a correct description of many non-equilibrium pro-cesses in nuclei and neutron stars as well. Acknowledgements I am grateful to my former studentsand many collaborators, whose names are in the papers we havepublished together and referenced here, and to many of my col-leagues with whom I have discussed many of the issues cov-ered in this review. The results presented here have been obtainedwith support by US DOE Grant No. DE-FG02-97ER-41014 andin part by NNSA cooperative agreement de-na0003841. Calcu-lations have been performed at the OLCF Titan and Piz Daintand for generating initial configurations for direct input into theTDDFT code at OLCF Titan and NERSC Edison. This researchused resources of the Oak Ridge Leadership Computing Facil-ity, which is a U.S. DOE Office of Science User Facility sup-ported under Contract No. DE- AC05-00OR22725 and of theNational Energy Research Scientific computing Center, which issupported by the Office of Science of the U.S. Department ofEnergy under Contract No. DE-AC02-05CH11231. We acknowl-edge PRACE for awarding us access to resource Piz Daint basedat the Swiss National Supercomputing Centre (CSCS), decisionNo. 2016153479. This work is supported by ”High PerformanceComputing Infrastructure” in Japan, Project ID: hp180048. Somesimulations were carried out on the Tsubame 3.0 supercomputerat Tokyo Institute of Technology. References [1] R. M. Dreizler and E. K. U. Gross, Density Functional The-ory: An Approach to the Quantum Many–Body Problem(Springer-Verlag, Berlin, 1990).[2] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965).[3] L. N. Oliveira, E. K. U. Gross, and W. Kohn, Phys. Rev.Lett. , 356 (1941).[6] I. M. Khalatnikov, An introduction to the theory of su-perfluidity (Advanced Book Program, Perseus Pub., Cam-bridge, Massachusetts, 2000).[7] H. Lamb, Hydrodynamics, 6 edition (Dover, New York,1945).[8] L. D. Landau and E. M. Lifshitz, Fluid Mechanics(Butterworth-Heinemann, Oxford, 1966).[9] R. Feynman, Prog. Low Temp. Phys. , 17 (1955).[10] W. F. Vinen and J. J. Niemela, J. Low Temp. Phys. , 167(2002).[11] E. P. Gross, Nuovo Cimento , 454 (1961).[12] L. P. Pitaevskii, Soviet JETP , 451 (1961).[13] E. Zaremba, T. Nikuni, and A. Griffin, J. Low Temp. Phys. , 277 (1999).[14] N. P. Proukakis and B. Jackson, Journal of Physics B:Atomic, Molecular and Optical Physics (20), 203002(2008).[15] A. Bulgac and S. Yoon, Phys. Rev. Lett. , 085302(2009).[16] W. Zwerger (ed.), The BCS–BEC Crossover and theUnitary Fermi Gas, Lecture Notes in Physics, Vol. 836(Springer-Verlag, Berlin Heidelberg, 2012).[17] P. W. Anderson, Phys. Rev. , 1900 (1958).[18] P. W. Anderson, Nature Phys. , 93 (2015). Copyright line will be provided by the publisher ss header will be provided by the publisher 13 [19] P. B. Littlewood and C. M. Varma, Phys. Rev. B (Nov),4883–4893 (1982).[20] A. F. Volkov and S. M. Kogan, Sov. Phys. JETP , 1018(1974).[21] P. Hohenberg and W. Kohn, Phys. Rev. , B864 (1964).[22] G. A. Baker Jr., Phys. Rev. C , 054311 (1999).[23] The Many-Body Challenge Problem (MBX) formulated byG. F. Bertsch in 1999. [22,16].[24] J. Carlson, S. Gandolfi, K. E. Schmidt, and S. Zhang, Phys.Rev. A , 061602 (2011).[25] M. J. H. Ku, A. T. Sommer, L. W. Cheuk, and M. W. Zwier-lein, Science , 563 (2012).[26] M. A. L. Marques, N. T. Maitra, F. M. S. Nogueira, E. K. U.Gross, and A. Rubio (eds.), Fundamentals of Time-Dependent Density Functional Theory, Lecture Notes inPhysics, Vol. 837 (Springer, Heidelberg, 2012).[27] A. Bulgac, Phys. Rev. A , 040502 (2007).[28] A. Bulgac, M. M. Forbes, and P. Magierski, The UnitaryFermi Gas: From Monte Carlo to Density Functionals,Vol. 836 of Zwerger [16], 2012, chap. 9, pp. 305 – 373.[29] A. Bulgac, Ann. Rev. Nucl. Part. Sci. , 97 (2013).[30] P. G. de Gennes, Superconductivity of metals and alloys(Benjamin, New York, 1966).[31] A. Bulgac and Y. Yu, Phys. Rev. Lett., 042504 (2002).[32] J. Carlson, S. Y. Chang, V. R. Pandharipande, and K. E.Schmidt, Phys. Rev. Lett. , 050401 (2003).[33] J. Carlson and S. Reddy, Phys. Rev. Lett. , 060401(2005).[34] S. Y. Chang, V. R. Pandharipande, J. Carlson, and K. E.Schmidt, Phys. Rev. A , 043602 (2004).[35] G. E. Astrakharchik, J. Boronat, J. Casulleras, andS. Giorgini, Phys. Rev. Lett. , 200404 (2004).[36] C. Lobo, A. Recati, S. Giorgini, and S. Stringari, Phys. Rev.Lett. , 200403 (2006).[37] J. Carlson and S. Reddy, Phys. Rev. Lett. , 150403(2008).[38] A. Schirotzek, Y. Shin, C. H. Schunck, and W. Ketterle,Phys. Rev. Lett. , 140403 (2008).[39] S. Pilati and S. Giorgini, Phys. Rev. Lett. (Jan), 030401(2008).[40] R. Combescot, A. Recati, C. Lobo, and F. Chevy, Phys.Rev. Lett. (May), 180402 (2007).[41] Y. i. Shin, Phys. Rev. A (Apr), 041603 (2008).[42] A. Bulgac and M. M. Forbes, Phys. Rev. Lett. , 215301(2008).[43] J. E. Thomas, J. Kinast, and A. Turlapov, Phys. Rev. Lett. , 120402 (2005).[44] F. Werner and Y. Castin, Phys. Rev. A (Nov), 053604(2006).[45] D. T. Son, Three comments on the Fermi gas at unitarity ina harmonic trap, arXiv:0707.1851 (2007).[46] A. Bulgac and M. M. Forbes, Phys. Rev. A , 031605(R)(2007).[47] A. Bulgac and M. M. Forbes, The Asymmetric Su-perfluid Local Density Approximation (ASLDA),arXiv:0808.1436, (2008).[48] D. Blume, J. von Stecher, and C. H. Greene, Phys. Rev.Lett. (23), 233201 (2007).[49] J. von Stecher, C. H. Greene, and D. Blume, Phys. Rev. A , 043619 (2008). [50] E. A. Yuzbashyan, O. Tsyplyatyev, and B. L. Altshuler,Phys. Rev. Lett. (Mar), 097005 (2006).[51] M. Bender, P. H. Heenen, and P. G. Reinhard (1), 121–180 (2003).[52] A. Bulgac, M. M. Forbes, S. Jin, R. N. Perez, andN. Schunck, Phys. Rev. C (Apr), 044313 (2018).[53] A. Bulgac and M. M. Forbes, Phys. Rev. C (5),051301(R) (2013).[54] A. Bulgac, M. M. Forbes, M. M. Kelley, K. J. Roche, andG. Wlazłowski, Phys. Rev. Lett. (Jan), 025301 (2014).[55] A. Bulgac, S. Jin, K. J. Roche, N. Schunck, and I. Stetcu,Fission Dynamics, 2018.[56] A. Bulgac and Y. Yu, Phys. Rev. Lett. , 190404 (2003).[57] Y. Yu and A. Bulgac, Phys. Rev. Lett. , 161101 (2003).[58] M. W. Zwierlein, J. R. Abo-Shaeer, A. Schirotzek, C. H.Schunck, and W. Ketterle, Nature , 1047–1051 (2005).[59] A. I. Larkin and Y. N. Ovchinnikov, Zh. Eksp. Teor. Fiz. , 1136 (1964).[60] P. Fulde and R. A. Ferrell, Phys. Rev. , A550 (1964).[61] A. Bulgac, Y. L. Luo, P. Magierski, K. J. Roche, and Y. Yu,Science , 1288 (2011).[62] P. Magierski, G. Wlazłowski, and A. Bulgac, Phys. Rev.Lett. (Sep), 145304 (2011).[63] R. Combescot, M. Y. Kagan, and S. Stringari, Phys. Rev. A (4), 042717 (2006).[64] R. Sensarma, M. Randeria, and T. L. Ho, Phys. Rev. Lett. (9), 090403 (2006).[65] S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod.Phys. , 1215 (2008).[66] Real-Time Dynamics of Quantized Vortices in a UnitaryFermi Gas, http://faculty.washington.edu/bulgac//UFG/.[67] A. Bulgac, Y. L. Luo, and K. J. Roche (15), 150401(2012).[68] J. A. Joseph, J. E. Thomas, M. Kulkarni, and A. G.Abanov (15), 150401 (2011).[69] T. Yefsah, A. T. Sommer, M. J. H. Ku, L. W. Cheuk, W. Ji,W. S. Bakr, and M. W. Zwierlein, Nature , 426 (2013).[70] A. Bulgac, M. M. Forbes, M. M. Kelley, K. J. Roche, andG. Wlazłowski, Phys. Rev. Lett. , 025301 (2014).[71] M. J. H. Ku, W. Ji, B. Mukherjee, E. Guardado-Sanchez,L. W. Cheuk, T. Yefsah, and M. W. Zwierlein, Phys. Rev.Lett. , 065301 (2014).[72] M. J. H. Ku, B. Mukherjee, T. Yefsah, and M. W. Zwierlein,Phys. Rev. Lett. , 045304 (2016).[73] G. Wlazłowski, A. Bulgac, M. M. Forbes, and K. J. Roche,Phys. Rev. A , 031602 (2015).[74] A. Bulgac, M. M. Forbes, and G. Wlazłowski, Journal ofPhysics B: Atomic, Molecular and Optical Physics ,014001 (2017).[75] G. Wlazłowski, K. Sekizawa, M. Marchwiany, andP. Magierski, Phys. Rev. Lett. , 253002 (2018).[76] H. Kamerlingh Onnes, KNAW, Proceedings 13 II , 1093–1113 (1910-1911).[77] J. Bardeen, L. N. Cooper, and J. R. Schrieffer (5), 1175–1204 (1957).[78] O. Hahn and F. Stassmann, Naturwissenschaften , 11(1939).[79] A. Bulgac, P. Magierski, K. J. Roche, and I. Stetcu, Phys.Rev. Lett. , 122504 (2016). Copyright line will be provided by the publisher [80] P. W. Anderson and N. Itoh, Nature (5512), 25–27(1975).[81] A. Bulgac, M. M. Forbes, and R. Sharma, Phys. Rev. Lett. , 241102 (2013).[82] G. Wlazlowski, K. Sekizawa, P. Magierski, and M. M. Bul-gacA.and Forbes., Phys. Rev. Lett. , 232701 (2016).[83] B. D. Josephson, Rev. Mod. Phys. , 216 (1964).[84] B. D. Josephson, Rev. Mod. Phys. , 251 (1974).[85] Y. Shin, M. Saba, T. A. Pasquini, W. Ketterle, D. E.Pritchard, and A. E. Leanhardt, Phys. Rev. Lett. , 050405(2004).[86] G. Wlazłowski, A. Bulgac, M. M. Forbes, and K. J. Roche,Phys. Rev. A , 031602(R) (2015).[87] A. Bulgac, M. M. Forbes, S. Jin, R. N. Perez, andN. Schunck, Phys. Rev. C , 044313 (2018).[88] P. Magierski, K. Sekizawa, and G. Wlazłowski, Phys. Rev.Lett. , 042501 (2017).[89] P. W. Anderson, Measurement in Quantum Theory and theProblem of Complex Systems, in: The Lesson of Quan-tum Theory, edited by J. d. Boer, E. Dal, and O. Ulfbeck,(North-Holland, Elsevier, Amsterdam, 1986), pp. 23–34.[90] A. Bulgac and S. Jin, Phys. Rev. Lett. , 052501 (2017).[91] J. W. Negele and H. Orland, Quantum Many-Particle Sys-tems (Westview Press, Reading, MA, November 1998).[92] D. L. Hill and J. A. Wheeler, Phys. Rev. (Mar), 1102–1145 (1953).[93] J. J. Griffin and J. A. Wheeler, Phys. Rev. (Oct), 311–327 (1957).[94] P. Fr¨obrich and I. Gontchar, Physics Reports (3), 131 –237 (1998).[95] A. Bulgac, M. M. Forbes, and G. Wlazłowski, J. Phys. B , 014001 (2017). Graphical Table of Contents