Time-Dependent Density Functional Theory for Fermionic Superfluids: from Cold Atomic Gases, to Nuclei and Neutron Stars Crust

TL;DR

This paper discusses the application of time-dependent density functional theory to fermionic superfluids across various systems, highlighting the role of large pairing gaps in phenomena from cold atomic gases to neutron star crusts.

Contribution

It introduces a theoretical framework for modeling superfluid fermionic systems with strong pairing correlations across different physical contexts.

Findings

Large pairing gaps influence quantum shock waves and turbulence.

Superfluid phenomena like vortex rings and domain walls are explained.

Pairing correlations affect nuclear fission dynamics and neutron star crust behavior.

Abstract

In cold atoms and in the crust of neutron stars the pairing gap can reach values comparable with the Fermi energy. 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 the size 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 large number of particle-particle configurations, which contribute to the formation of the Copper pairs. These systems have been shown to be the host of a large number of remarkable phenomena, in which the large magnitude of 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 fragmented condensates and role of pairing correlations in collisions on heavy-ions, Larkin-Ovchinnikov phase as an example of a Fermi supersolid, role pairing correlations control the dynamics of fissioning nuclei, self-bound superfluid fermion droplets of extremely low densities.

Figures

Equations

$$i\hslash \frac{\partial \psi (\mathbf{r},t)}{\partial t}$$

$$|\Delta (t)=|{\Delta }_{\infty }|+\frac{A}{\sqrt{2|{\Delta }_{\infty }|t}}\sin (2|{\Delta }_{\infty }|t+\theta ),$$

$$f=\frac{1}{-ik-\frac{1}{a}+\frac{1}{2}{r}_0{k}^2+\dots },$$

$$\psi (\mathbf{r})=\exp (\mathbf{r}\cdot \mathbf{k})+\frac{f}{r}\exp (ikr)\approx 1-\frac{a}{r}+\mathcal{O}(kr).$$

$$\psi (\mathbf{r})=\frac{1}{\sqrt{2a}r}\exp \left(-\frac{r}{a}\right)+\mathcal{O}\left(\frac{r_0}{a}\right).$$

$$r_0\ll n^{-1/3}\ll |a|$$

$$E_{\text{gs}}(N,V,\hslash ,m,a,r_0)\to \xi N\frac{3}{5}{\epsilon }_F,$$

$$\mathcal{E}(\mathbf{r},t)=\frac{{\hslash }^2}{2m}[\alpha \tau (\mathbf{r},t)+\beta \frac{3}{5}(3{\pi }^2{)}^{3/2}{n}^{5/3}(\mathbf{r},t)$$

$$+\gamma \frac{|\nu (\mathbf{r},t){|}^2}{n^{1/3}(\mathbf{r},t)}+(1-\alpha )\frac{{\mathbf{j}}^2(\mathbf{r},t)}{n(\mathbf{r},t)}]+V_{\text{ext}}(\mathbf{r},t)n(\mathbf{r},t),$$

$$i\hbar\frac{\partial}{\partial t}\left(\begin{array}[]{c}u_{n}\\ v_{n}\end{array}\right)=\left(\begin{array}[]{cc}h-\mu&\Delta\\ \Delta^{*}&-h^{*}+\mu\end{array}\right)\left(\begin{array}[]{c}u_{n}\\ v_{n}\end{array}\right),$$

$$\nu ({\mathbf{r}}_1,{\mathbf{r}}_2,t)\propto \frac{1}{|{\mathbf{r}}_1-{\mathbf{r}}_2|},\text{if}|{\mathbf{r}}_1-{\mathbf{r}}_2|\to 0.$$

$$n(\mathbf{r},t)$$

$$\tau (\mathbf{r},t)$$

$$\nu (\mathbf{r},t)$$

$$\mathbf{j}(\mathbf{r},t)$$

$$h=-\alpha \frac{{\hslash }^2}{2m}\Delta +U(\mathbf{r}).$$

$$\alpha \frac{{\hslash }^2{k}_0(\mathbf{r})}{2m}+U(\mathbf{r})-\mu =0,$$

$$\alpha \frac{{\hslash }^2{k}_c(\mathbf{r})}{2m}+U(\mathbf{r})-\mu =E_{\text{cut}},$$

$$\frac{1}{g_{\text{eff}}(\mathbf{r})}=\frac{n^{1/3}(\mathbf{r})}{\gamma }$$

$$-\frac{m{k}_c(\mathbf{r})}{2{\pi }^2{\hslash }^2\alpha }\left[1-\frac{k_0(\mathbf{r})}{2k_c(\mathbf{r})}\ln \frac{k_c(\mathbf{r})+k_0(\mathbf{r})}{k_c(\mathbf{r})-k_0(\mathbf{r})}\right],$$

$$\Delta (\mathbf{r})=-g_{\text{eff}}(\mathbf{r})\nu (\mathbf{r}).$$

$$\alpha \frac{{\hslash }^2}{2m}\tau (\mathbf{r})-\Delta (\mathbf{r}){\nu }^{\ast }(\mathbf{r})$$

$$\mathcal{E}(n_a,n_b)=\frac{3}{5}\frac{{\hslash }^2}{2m}(6{\pi }^2{)}^{2/3}{\left[n_{a}g\left(\frac{n_b}{n_a}\right)\right]}^{5/3}$$

$$\displaystyle i\hbar\left(\begin{array}[]{c}\dot{u}_{k\uparrow}\\ \dot{u}_{k\downarrow}\\ \dot{v}_{k\uparrow}\\ \dot{v}_{k\downarrow}\end{array}\right)=\left(\begin{array}[]{cccc}h_{\uparrow\uparrow}&h_{\uparrow\downarrow}&0&\Delta\\ h_{\downarrow\uparrow}&h_{\downarrow\downarrow}&-\Delta&0\\ 0&-\Delta^{*}&-h^{*}_{\uparrow\uparrow}&-h^{*}_{\uparrow\downarrow}\\ \Delta^{*}&0&-h^{*}_{\downarrow\uparrow}&-h^{*}_{\downarrow\downarrow}\end{array}\right)\left(\begin{array}[]{c}u_{k\uparrow}\\ u_{k\downarrow}\\ v_{k\uparrow}\\ v_{k\downarrow}\end{array}\right).$$

$$\hat{U}(t_i,t_f)$$

$$\Psi (t)=$$

$$\int \mathcal{D}[\eta (t)]\bar{W}[\eta (t)]\exp \left(-\frac{i}{\hslash }\hat{h}[\bar{\sigma }(t)+\eta (t)]\right)\Psi (0),$$

$$i\hslash {\dot{\psi }}_k(\mathbf{r},t)=h[n]{\psi }_k(\mathbf{r},t)+\gamma [n]\dot{n}(\mathbf{r},t){\psi }_k(\mathbf{r},t)$$

$$-\frac{1}{2}\left[\mathbf{u}(\mathbf{r},t)\cdot \hat{\mathbf{p}}+\hat{\mathbf{p}}\cdot \mathbf{u}(\mathbf{r},t)\right]{\psi }_k(\mathbf{r},t)+u_0(\mathbf{r},t){\psi }_k(\mathbf{r},t),$$

Full text

1 Introduction.

There is a large class of fermionic superfluid systems in which the pairing correlations are very strong and their description and especially the description of their dynamics and interaction with typically strong external probes require and extension of the Density Functional Theory (DFT) [2] following the à la Kohn and Sham formulation [3], which does not involve non-local potentials, as in the first extension suggested by Oliveira et al. [4]. If the pairing gap is large, the number of particle-particle configurations defining the anomalous density is much larger than the number $N$ of fermions in the system. In the Kohn-Sham version of the DFT the energy density functional (EDF) depends on mainly two types of densities, the number density and the kinetic energy density, which are expressed through $N$ single particle wave functions (à la Hartree-Fock approximation). Since in the language of the number density alone one cannot distinguish between normal and superfluid phases, there is an obvious need to introduce the anomalous density as well, which defines the order parameter, non-vanishing only in the superfluid phase. However, the description of fermionic superfluids becomes even more demanding, since in practice one has to study very often superfluids in interaction with strong time-dependent probes, e.g. when one is stirring a fermionic superfluid, when studying non-equilibrium phenomena such as quantum turbulence, and when one may even observe the evolution of superfluid into a normal phase. Naturally, under such circumstances one needs another “order parameter,” capable to disentangle parts of the system which evolve in time at various rates, and in this case the appearance of current densities in the EDF becomes unavoidable. The physical systems of interests run the gamut from cold atom system to nuclear systems and one has to consider multi-component systems for which the structure the EDFbecomes quite complex. In nuclei one has to consider at the same time both the (charged) proton and neutron miscible superfluids and in neutron star crust in addition also include the electron background as well. In the case of cold atoms one is interested lately in miscible mixtures of either Fermi-Bose, Fermi-Fermi or Bose-Bose superfluids. In neutron stars various mesons, which are bosons, are expected to appear at relatively large densities, close to the core of the star.

2 Why is there a need for an extension of DFT to time-dependent and superfluid?

Two prevailing theoretical models are used to describe the dynamics of superfluids. The oldest is the Landau two-fluid hydrodynamics [6, 7], which at zero-temperature reduces to the hydrodynamics of a single perfect classical fluid [8, 9], namely of the superfluid component alone. Naturally in such a formalism Planck’s constant is absent and the two-fluid hydrodynamics is unable to describe the formation of quantized vortices [10], their dynamics, their crossing and recombination, which is at the heart of the venerable field of quantum turbulence [11] conjectured by Feyman in 1955. Classical turbulence is due to viscosity, which is absent in superfluids at zero-temperature. On the same note, there is no mechanism within the two-fluid hydrodynamics to describe either the conversion of the superfluid into the normal component, when a superfluid is stirred vigorously. The quantizations of vortices in the two-fluid hydrodynamics has to be enforced by hand [7], in a manner similar to the Bohr 1913 quantization of the hydrogen atom. It is thus impossible to describe the evolution of a superfluid at rest and brought into rotation, when quantized vortices are formed, and later on when they might cross and reconnect as well.

An extremely attractive alternative approach was developed by Gross [12] and Pitaevskii [13] to describe a weakly interacting Bose at zero-temperature, the celebrated Gross-Pitaevskii equation. This is a non-linear Schrödinger equation in which both the density and the superfluid order parameter are described by the same complex field $\psi({\bf r},t)$:

[TABLE]

where $g>0$ is the strength of the weak interparticle repulsion and $V_{\text{ext}}({\bf r},t)$ is some external potential. While this fully quantum mechanical formalism is adequate to describe a large range of properties of a weakly repulsive Bose gas, in particular the formation of quantized vortices, it is clearly inadequate to describe the properties of superfluid liquid 4He, which is a strongly interacting system. Under the influence of a strong time-dependent external field a weakly interacting Bose gas can become normal, at least in some regions of space, but the Gross-Pitaevskii equation is unable to disentangle between the normal and the superfluid components.

Zaremba, Nikuni, and Griffin[14] presented a nice solution to many of these issues in the case of weakly interacting bosons. They coupled the quantum Gross-Pitaevskii equation describing the condensate, with a classical kinetic equation for the bosons in the normal component, which includes collision between non-condensed bosons and also collision between the condensed and non-condensed bosons. Landau two fluid hydrodynamics emerges when the frequency of the mode satisfies the condition $\omega\tau\ll 1$, where $\tau$ is the relaxation time associated with the collisions between the condensed and non-condensed bosons. This new framework is valid only for weakly interacting bosons. Later on other more sophisticated approaches have been suggested for the weakly interacting Bose systems [15].

Apart from these difficulties discussed above, none of these models of superfluids allow for the existence of the Anderson-Bogoliubov-Higgs mode. This mode was noticed by P. W. Anderson a long time ago [18, 19, 20]. The potential energy of a system with a complex order parameter has a shape similar to a Mexican hat, see Fig. 1 taken from Ref. [5]). Typically this potential depends only on the magnitude of the complex field but not on its phase, e.g. $V(\phi)=a|\phi|^{2}+b|\phi|^{4}$, where $a<0$ and $b>0$ and a minimum value at $|\phi_{0}|^{2}=-b/2a$. The mode $|\phi_{0}|$ is known in high-energy physics as the Higgs boson and it became an essential element of the Standard Model. The existence of the Higgs boson leads to masses of quarks, gluons, electrons, Z0, and W± bosons, and other elementary particles and its existence has been determined experimentally [5].

If in a homogeneous unitary Fermi Gas (see next section for its characterization) one would bring very slowly the pairing gap out of the equilibrium position $\Delta_{0}$ in the ground state one would be able to observe at least two kinds of rather unexpected excitation modes. The first kind of excitation correspond to small amplitude oscillations around the equilibrium, with an unexpected slow algebraic damping [21, 16], when the magnitude of the pairing gap behaves as a function of time as

[TABLE]

where $|\Delta_{\infty}|<|\Delta_{0}$. The asymptotic state corresponds to a partially fermionic paired state plus quasi-particle excitations. If instead the pairing gap is slowly brought to a very small value and after that the system is left to evolve freely, the magnitude of the pairing gap oscillates with a maximum value smaller than the equilibrium value, $0<|\Delta(t)|<|\Delta_{0}|$, see Fig. 2. The same results illustrated in Fig. 2 can be obtained by preparing initially the system with an interaction strength corresponding to an equilibrium pairing gap smaller than $|\Delta_{0}|$ and suddenly changing the coupling strengths to a value corresponding to an equilibrium value of the pairing gap $|\Delta_{0}|$.

3 Energy density functional for the unitary Fermi gas.

We will illustrate the extension of the Kohn-Sham local density approximation (LDA) of the DFT to superfluid fermionic systems at first with the case of an unitary Fermi gas, which is both methodologically clear and of great practical value. In a subsection section we will later briefly discuss the structure of the EDF for nuclei and neutron stars.

As in the case of normal electron systems, one of the most frustrating aspects of DFT is the construction of the EDF, which we know it exits according to the Hohenberg-Kohn theorem [22], but for which however we have no well defined procedure to construct with enough high accuracy. The unitary Fermi gas [23, 17], which is a system of spin-up and down fermions, interacting with a zero-range potential, characterized by an infinite scattering length and a zero range effective become an object of extremely intensive study both experimentally and theoretically in the last two decades. George Bertsch noticed [24] that the neutron matter in the crust of neutron stars is very close to such an idealized system.

Neutron interaction in the $s$-wave is characterized by a very large scattering length $a$ and a relatively small effective range $r_{0}$, with an low-energy $s$-wave scattering amplitude

[TABLE]

where $k$ is the relative wave vector of two scattering fermions. The wave function outside the potential is

[TABLE]

The total scattering cross section of two low-energy fermions $\sigma=4\pi|f|^{2}\rightarrow\tfrac{4\pi}{k^{2}}$ reaches the maximum possible value allowed by unitarity if $r_{0}\rightarrow 0$ and $|a|\rightarrow\infty$. If the scattering length $a$ is positive then a bound state exist with the wave function

[TABLE]

If in a many-fermion system, in which the average interparticle separation is $\propto n^{-1/3}$, the conditions,

[TABLE]

are satisfied such a system is called a unitary Fermi gas (UFG) [17].

In 1999 the dilute neutron matter (for which the condition (6) is weakly satisfied) was he closest physical system to an UFG one could envision, and the calculation of the value of the dimensionless Bertsch parameter $\xi$ became a theoretical challenge. If $\xi<0$ the system would collapse into a high density liquid or solid with an average interparticle separation likely of the order of the range of the interaction and with properties determined by the particular features of the interaction between fermions. However, if $0<\xi\leq 1$ the ground state would be that of a gas, and a very unusual gas at that, a superfluid with a pairing gap of the order of the Fermi energy, the largest pairing gap in any physical system, and universal properties largely independent of the details of the interaction. One can easily establish using dimensional arguments alone that the ground state energy of a uniform UFG should be given by a function depending only on the volume $V$, particle number $N$, $\hbar$, and the fermion mass $M$, in the unitary limit when $r_{0}\rightarrow 0$ and $|a|\rightarrow\infty$:

[TABLE]

where $\varepsilon_{F}=\tfrac{\hbar^{2}k_{F}^{2}}{2m}$ is the Fermi energy of a uniform non-interacting Fermi gas with the Fermi wave-vector $k_{F}=\left(3\pi^{2}\tfrac{N}{V}\right)^{1/3}$, and $\xi$ is the dimensionless Bertsch parameter. $\xi\leq 1$, as the scattering length $a$ can become infinite only in the case of attractive interactions. Apart from the universal value of the Bertsch parameter $\xi$, the EDF for the UFG should have the same structure as for a non-interacting Fermi gas and the only parameter specifying the nature of the fermions is their mass. Both theoretical Quantum Monte Carlo (QMC) $\xi=0.372(5)$ [25] and experimental $\xi=0.376(5)$ [26] values of the Bertsch parameter are now in very good agreement with each other.

A fermionic superfluid under the influence of a time-dependent external field might become normal in some spatial regions but not in others, and the number density alone cannot discriminate between different phases. The case of the UFG is particularly attractive from the point of view of a DFT aficionado, as only rather general requirements suffice to narrow down the structure of the EDF. Dimensional arguments, rotation, translational, and parity symmetries, gauge symmetry related to the transformations of the complex order parameter, Galilean covariance, and renormalization and regularization of the EDF combined with the so called adiabatic local density approximation (ALDA) [27] extended to the case of superfluid systems restrict the form of the (unregulated) EDF for UFG to a rather simple form [28, 29, 30], namely:

[TABLE]

where $\alpha$, $\beta$, and $\gamma$ are dimensionless constants. $n({\bm{r}},t)$, $\tau({\bm{r}},t)$, $\nu({\bm{r}},t)$, and ${\bm{j}}({\bm{r}},t)$ are the unregulated number, kinetic, anomalous densities and current densities of a fully unpolarized UFG (equal number of spin-up and spin-down particles) and expressed through the Bogoliubov quasiparticle amplitudes $u_{n}({\bf r},t)$ and $v_{n}({\bf r},t)$ [31]. $V_{\text{ext}}({\bm{r}}t)$ is an arbitrary external field with which one might probe or excite the system. We refer to this form of DFT as the (Time-Dependent) Superfluid Local Density Approximation ((TD)SLDA), which is a natural extension of the LDA for normal systems of Kohn and Sham formulation of the DFT [3] to superfluid systems. The emerging TDSLDA equations have the expected form, identical to the Bogoliubov-de Gennes equations,

[TABLE]

where $h=\tfrac{\partial{\cal E}}{\partial n}$, $\Delta=\tfrac{\partial{\cal E}}{\partial\nu^{*}}$, and $\mu$ are the single-particle Hamiltonian, the pairing field, and the chemical potential. Both kinetic energy and anomalous densities diverge as a function of the ultraviolet cutoff in a very similar manner in the case of a zero-range interaction [29]. In particular the anomalous density matrix has the behavior

[TABLE]

This divergence is the same as the divergence of the scattering wave in Eqs. (4) or of the bound state wave function (5) when $r\rightarrow r_{0}\rightarrow 0$. This divergence is “real” and not a deficiency of the theoretical formalism. The divergence of the anomalous density matrix reflects nothing else but the increase of the wave function of the Cooper pair when the separation between the fermions approaches the radius of the interaction. One can relate the anomalous density matrix $\nu({\bm{r}}_{1},{\bm{r}}_{2},t)$ with the Cooper pair wave function.

The mathematical difficulty in extending DFT in á la Kohn and Sham manner to superfluid system arises when one attempts to reach the limit $n^{1/3}r_{0}\rightarrow 0$, when the average interparticle separation is much larger than the radius of the interaction and avoid the appearance of infinities in the calculations of various densities and of the pairing gap. In meanfield approximation the pairing field $\Delta({\bm{r}}_{1},{\bm{r}}_{2},t)=-V({\bm{r}}_{1}-{\bm{r}}_{2})\nu({\bm{r}}_{1},{\bm{r}}_{2},t)$, where $V({\bm{r}}_{1}-{\bm{r}}_{2})$ is the fermion-fermion interaction responsible for the pairing correlations, is formally non-local. In the limit $n^{1/3}r_{0}\rightarrow 0$ densities should be calculated with a cutoff [32, 29]:

[TABLE]

in which the pre-factor 2 is for the spin multiplicity. In time-dependent simulations one usually starts with the system prepared in its ground state and the initial quasi-particle amplitudes $u_{n}({\bf r},t)=u_{n}({\bf r})\exp\left(-i\tfrac{E_{n}t}{\hbar}\right)$ and $v_{n}({\bf r},t)=v_{n}({\bf r})\exp\left(-i\tfrac{E_{n}t}{\hbar}\right)$. $E_{n}$ are the eigenvalues of the initial SLDA stationary equations and $E_{\text{cut}}$ is an ultraviolet cut-off energy, which when chosen large enough does not affect the values of any physical observables [32]. In the ground state the single-particle Hamiltonian of an unpolarized UFG has the structure

[TABLE]

One needs to introduce the momentum dependent wave-vectors and the renormalized coupling constant

[TABLE]

in roder to derive the renormalized form of the pairing gap

[TABLE]

One can then show that the combination

[TABLE]

does not diverge when $E_{\text{cut}}\rightarrow\infty$. While the wave vector $k_{0}({\bm{r}})$ can become imaginary in the classically forbidden regions of space, the effective coupling consnat $g_{\text{eff}}({\bm{r}})$ remains real. If $k_{c}({\bm{r}})$ becomes imaginary in any spatial region the recipe is that the last term on the rhs side of Eq. (18) should be dropped [32, 29].

The EDF for the UFG Eq. (8) depends on three dimensionless parameters $\alpha,\beta$ and $\gamma$, for both superfluid and normal phases, which can be extracted from values of the Bertsch parameter $\xi$, the pairing gap, and the momentum dependence of the quasi-particle excitations obtained in the QMC for the uniform UFG [33, 34, 25, 35, 36, 37], which agree well with extracted experimental values [38, 39, 26].

In the case of polarized UFG one needs two number densities $n_{a,b}({\bm{r}},t)$ and from dimensional arguments alone it follows that the energy density of a uniform polarized UFG is:

[TABLE]

with $g(1)=(2\xi)^{3/5}$, see Fig. 3.

3.1 Validation of (TD)SLDA.

At this point one is in the position to validate the accuracy of the suggested form of the EDF for the UFG (8) by placing various number of fermions with spin-up and spin-down in external fields. A UFG in a harmonic trap is particularly interesting, since from theory we know that its properties are determined by the only energy scale in the system $\hbar\omega$ [44, 45, 46], their properties can be calculated numerically with controlled accuracy within QMC and unlike the the homogeneous system density gradients are important. The EDF for the UFG has so far been constructed for uniform systems and it is not obvious that such an EDF would perform well for inhomogeneous systems, where density gradients are significant. Fortunately there exist a sufficiently large number of QMC calculations of unpolarized and polarized finite systems of fermions at unitarity in a external harmonic trap, both in the superfluid and normal phases. The EDF for 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 as those used for the derivation of the EDF for an unpolarized UFG can be used for its derivation [47, 43, 48, 29]. In case of inhomogeneous systems a new gradient term might need to be included in Eq. (8), namely one of the form $\propto\left|{\bm{\nabla}}\sqrt{n({\bf r},t)}\right|^{2}$, but comparison with QMC calculations of inhomogeneous systems indicate that it can be neglected [28]. Polarized fermions at unitarity can be in either superfluid or normal phase, depending on the degree of polarization $N_{\downarrow}/N_{\uparrow}$ [47, 43, 48, 29] and the same EDF can be used to describe with surprisingly high accuracy all such systems. When comparing the results of the QMC calculations performed for trapped unitary fermions [49, 50] with the predictions of the SLDA one finds almost perfect agreement [28, 29]. The differences between the QMC calculations and the SLDA predictions are almost always within the statistical errors of the QMC calculations, within at most 1-3%. There are two exceptions. The QMC and SLDA calculations for the $(N_{\uparrow},N_{\downarrow})=(15,14)$ have an error of 9.5%, attributed to the inaccuracies in the QMC results for this particular system. The largest disagreement is for the two-fermion system $(N_{\uparrow},N_{\downarrow})=(1,1)$, about 15%. The rest of the differences between the QMC results and the SLDA predictions are at the level of $1\dots 3$ %, which is also the level of accuracy of the QMC calculations. Surprisingly, the odd-even staggering, namely the energy differences between the ground state energies of even $N_{\uparrow}=N_{\downarrow}$ and odd systems $N_{\uparrow}-N_{\downarrow}=1$ are within statistical errors as well. One should keep in mind the ground state energies of harmonically unpolarized trapped systems of unitary fermions scale as $\sqrt{\xi}E_{NI}\propto\hbar\omega N^{4/3},$ where $E_{NI}$ is the energy of $N$ non-interacting fermions in a harmonic trap. The odd-even energy differences depend relatively weakly on the total particle number [49, 50] and the relatively small odd-even energies differences are reproduced within SLDA with surprising accuracy as well.

There are only a few exact solutions of the time-dependent Schrödinger equation for interacting many-fermion systems. The linear response theory predicts damped harmonic oscillations with a frequency $2|\Delta_{0}|$ [18, 21], while Eq. (2) emerges when nonlinearities are taken into account. Equally unexpected is the behavior of the pairing gap when the initial disturbance is large, panels $a$ and $b$ in Fig. 2, when one would naively expect that the pairing gap would oscillate somehow around the equilibrium value $|\Delta_{0}|$. While these modes have been put in evidence in a meanfield-like framework and one might suspect that they are artifacts of possible approximations, the same behavior was demonstrated to appear in an exactly soluble time-dependent realistic many-body model of superconductivity, the Richardson or Gaudin model [51]. We will illustrate below however the power of TDSLDA in confronting actual non-equilibrium phenomena in real experiments.

3.2 Nuclear systems.

Nuclear systems are significantly more complex that the UFG. While we know for decades quite a bit about the nuclear EDF (NEDF), the exact from and accuracy of mostly empirical NEDF is still insufficient for many applications, like predicting the origin of elements in Universe. In case of nuclear systems there two type of fermions, protons and neutrons, and the spin-orbit interaction is very strong, and has the opposite sig when compared to atomic systems. Pairing correlations are relatively strong as well, though not as strong as in the case of UFG. However, the dilute neutron matter, found in the crust of neutron stars, has quite a lot of similarities with the UFG [24, 23]. QMC calculations for nuclear systems a significantly more complex than for electronic or cold atom systems, and so far only the pure neutron matter equation of state as function of density is know with reasonable accuracy from  it ab initio calculations. The interaction between nucleons (neutrons and protons) is very complex and in nuclei not only two-body interactions are important, but also three-body (and even four-body) interactions play a crucial role. As a result the form of NEDF is typically obtained with significant phenomenological input [52, 53]. The accuracy of phenomenological NEDFs in predicting the binding energies of about 2300 known nuclei is nowadays at the sub-percent level, but still not accurate enough for many applications.

In the resulting evolution equations we suppressed for the sake of simplicity the space and time coordinates $({\bm{r}},t)$. The ensuing equations represent an infinite set of coupled nonlinear time-dependent 3D partial differential equations for the quasi-particle wave functions,

[TABLE]

Here both the local mean field $h_{\sigma\sigma^{\prime}}$ and pairing field $\Delta$ depend on the various single-particle densities. The index $k$ labels each quasi-particle wave function, and is both discrete and continuous. This index must also run over isospin, so that there are similar sets of equations for both protons and neutrons, which naturally are coupled. We have explicitly included the spin indices $(\sigma,\sigma^{\prime})\in\{\uparrow,\downarrow\}$, allowing for mixing between the spin-up and spin-down states by the spin-orbit interaction, thus capturing the effects of proton-proton and neutron-neutron pairing.

3.3 Numerical implementation.

The TDSLDA equations are discretized in space and time. The system of interest is placed on a spatial lattice with $N_{x}N_{y}N_{z}$ lattice points, for a chosen lattice constant $l$, which determines the momentum cutoff $\hbar\pi/l$ [54]. When these simulation box parameters are chosen appropriately one can ensure that with further discretization the corrections are (exponentially) small [54]. The equations are propagated in time using a 6th order multi-step predictor-corrector-modifier, which typically requires an application of the quasi-particle Hamiltonian only twice per time-step. Notice that in the case of the popular Runge-Kutta 4th order method one would need to apply the Hamiltonian four times per time step. Unitarity and accuracy of various conserved quantities during the time evolution are satisfied with high accuracy for up to millions of time steps. The number of coupled complex 3D time-dependent nonlinear partial differential equations which has to be evolved in time is up to ${\cal O}(10^{6})$, see supplemental material in Refs. [29, 55, 56]. The numerical solution of these equations requires the use of the leading edge supercomputers and it ranks among some of the largest direct numerical simulations attempted so far.

4 SLDA for cold Atomic gases, nuclei, and neutron star crust.

This section will briefly describe qualitative aspects of static and dynamic aspects of fermionic superfluids in strongly interacting cold atomic gases, nuclei and neutron star crust.

4.1 Quantized vortex in cold superfluid Fermi gases and dilute neutron matter.

One of the first applications of the SLDA for UFG was to determine the structure of a quantized vortex [57, 58]. It was shown that the actual number density profile of a quantized vortex is somewhere between that in a BCS superfluid and in a BEC superfluid. While in a BCS superfluid the density in the core is practically the same as the density far away from the vortex core and only the order parameter vanishes at the vortex core, in a BEC quantized vortex both the order parameter and the number density vanish at the core, in the case of dilute superfluids. In the case of the UFG while the order parameter vanishes at the core, the vortex core is only partially filled with fermions, with a density only about half the value of the asymptotic value. This density depletion of the vortex core was used to visualize in experiments the Abrikosov vortex lattice formed in a rotating UFG [59], which was the deciding experimental argument in demonstrating that UFG is indeed a superfluid.

4.2 The Larkin-Ovchinnikov phase.

Once the EDF of the UFG was established and validated, it was used to in order to establish if a UFG can sustain an inhomogeneous state of the order parameter [60, 61]. In Ref. [43] the selfconsistent SLDA equations for a system in which the order parameter was allowed to oscillate in the $z$-direction, while in the $x$ and $y$ directions the properties of the system remains homogeneous, were solved Unlike in the case of a weakly coupled BCS superfluid, where the LO phase exist only for a very narrow window of spin polarization, the UFG was shown to sustain such a phase in a surprisingly wide range, see Fig. 3. Only the amplitude of the oscillations of the number density of the minority component are significant. This LO phase has not been observed yet in cold fermi atom systems.

4.3 Real-time generation and dynamics of quantized vortices.

The most fascinating applications of the SLDA is to non-equilibrium phenomena. In the first simulation of a fermionic superfluid in real time [62] we placed a UFG in a container resembling a soda can, homogeneous and with periodic boundary conditions along the longitudinal direction, along with other geometries as well. We inserted a “straw” and started stirring the fluid with various constant angular velocities and linear velocities smaller and larger than the Landau critical velocity. A UFG, apart from being characterized by a very large pairing gap, which in appropriate units is even bigger than in high $T_{c}$ superconductor [63]. An UFG one has perhaps the largest Landau critical velocity (in appropriate units) of any superfluid [64, 65, 66]. Many of the results obtained are available in the form of videos online, for various geometries, various ways to stir the superfluid, and a range of stirring velocities ranging from very slow to well above the Landau critical velocity [67]. If one introduces an object and moves it relatively slowly through the superfluid and eventually extract that object slowly also the UFG returns practically to its initial state, as one would have naturally expected for an adiabatic evolution. We have also noticed that sometimes we can bring the atomic cloud into rotation even if the linear speed of the “straw” exceeds Landau critical velocity and the cloud remained superfluid. This is possible since UFG is gas, during rotation accumulates along the wall, the density and therefore the local Landau critical velocity increases. In the same work [62] we have demonstrated in a real-time treatment for the first time that in a fermionic superfluid quantized vortices can cross and recombine, exactly as Feynman [10] envisioned and suggested that quantum turbulence emerges, see also subsection 4.6.

4.4 Quantum shock waves.

Thomas and collaborators demonstrated that quantum shock waves can be excited in a cold atomic fermionic cloud [Jpseph:2011]. Thee shock wave front was directly visualized as a clear number density discontinuity. In a TDSLDA simulation of this experiment [68] is was shown that the character of these shock waves is controlled by the dispersive effects and not dissipative effects as was assumed in the analysis of the experiment [69]. At the front the shock wave both the number density and the velocity field are discontinuous, when coarse grained appropriately. At the front of a shock wave the matter flows in opposite directions. But in addition with quantum shock waves also domain walls are formed, which propagate thought the cloud, see Ref. [68] and Fig. 6. A domain wall is the region where the phase of the order parameter changes from $-\pi$ to $\pi$ and is topological in character, similar to domain walls in magnetism. At a domain wall number density has a significant depletion, as in a vortex core.

4.5 Vortex ring versus heavy soliton.

In a surprising experiment performed in 2013 [70] it was announced that new kind of many-body excitation was observed in an experiment performed on an UFG. Half of a very elongated superfluid cloud was illuminated with a laser, which resulted in a difference in the phase of the order parameter between the two halves. Theoretically is was known for quite some time that in this case such a planar soliton propagates with a know speed, and in 3D can be unstable due to the snake-instability. What the MIT experiment established was this the excitations they created was moving with a speed two orders of magnitude slower than the theory predicted and the dubbed this mode a “heavy soliton.” In a very careful TDSLDA analysis [71] we established however, that under the conditions described in the original paper [70] most likely the authors observed a vortex ring. A planar soliton in an inhomogeneous cloud quite rapidly evolves into a vortex ring in simulations and it starts propagating back and forth in an almost perfect harmonic motion. However, while moving in one direction the vortex ring is large, and while is moving in the opposite direction the vortex ring shrinks, and then the motion is repeated, see Fig.7. In subsequent more detailed experiments [72, 73] the MIT group confirmed the transformation of the planar soliton, into a vortex ring. However, since their trap lack azimuthal symmetry, which was not adequately established in the original paper, a vortex ring in an axially non-symmetric trap rather quickly touches the walls and turns either into a single or two line vortices, dubbed however or solitonic vortices, if vortex line is perpendicular to the propagation direction. This is fully in agreement with earlier simulations of an UFG [68], see supplemental material [67], and subsequent analysis of the new MIT experiment [74, 75].

4.6 Quantum turbulence.

At the hearty of quantum turbulence lie the formation of a tangle many quantized vortices. In experiments with liquid helium such tangles have been created in laboratories for decades [11]. In the case of cold fermionic gases vortices can be easily generated with rotating laser beams [59] but also by simply illuminating part of a cloud with a laser [70]. By combining these two methods, which clearly are not the only possibility, one can generate a tangle of quantized vortices in both unpolarized [71, 75] and Fig. 8, and in polarized systems [76]. The great advantages of cold atom system over liquid helium are multiple: i) one can control easily many parameters of the system, including the interaction strength, and create a wide array of external probes [17]; ii) one can describe with very good accuracy theoretically both static properties and particularly the non-equilibrium dynamics of such systems with very good control of the theoretical ingredients. being able to confront theory and experiment in great detail is a great advantage of cold atom systems over studies performed in liquid helium, where only phenomenological models exist.

4.7 Nuclear fission.

Our main goal in developing SLDA was to describe non-equilibrium dynamics of nuclear systems, but in order to verify and validate the theoretical framework and the complex numerical implementation we made quite a long detour studying cold atoms systems, for which theoretical tools are both simpler and in better control, and also there is a great wealth of experimental data which can be confronted with theory predictions and postdictions! One of the oldest problems of strongly interacting quantum many-body systems i s nuclear fission. Unlike superconductivity for example, which required less than five decades form the initial experimental observation [77] in 1911, until a microscopic theory was put forward in 1957 [78], nuclear fission observed in 1939 [79] will turn 80 years old in 2019 and likely will still have no adequate microscopic description. There are many reasons why this is the case and here are some of them: i) the nuclear interactions are extremely complex and they are not yet accurately known, unlike Coulomb interaction between electrons and nuclei; ii) nuclei are finite systems and at the same time they have too many particles, correlations are very stroong; iii) when a heavy nucleus fissions the number of final channels is in the hundreds, corresponding to various possible splittings of protons and neutrons between the fission fragments, the fission fragments emerge excited with various quantum numbers; iv) apart from fission fragments typically quite a number of neutrons are emitted, during fission, immediately after fission, and much later, along with gamma rays and beta decays. One cannot declare “victory” until theory is able to predict with reasonable confidence most of these fission fragment properties, which are critical for many applications as well, and also for clarifying many fundamental questions, such as the origin of elements in the Universe, and the structure and evolution of stars.

Using TDSLDA and an NEDF of reasonable quality in 2016 we were able for the first time for describe the evolution of a fissioning nucleus from the outer fission barrier to scission, until fission fragments were separated [80] and Fig.9 and recently also in Ref. [56]. The fission fragments emerge with properties similar to those determined experimentally, while the fission dynamics appears to be quite complex, with various shape and pairing modes being excited during the evolution. The time scales of the evolution are found to be much slower than previously expected and the role of the collective inertia in the dynamics is found to be negligible. Even though in this first study of its kind we did not obtain a perfect agreement with experiment, our results clearly demonstrate that rather complex calculations of the real-time fission dynamics without any restrictions are feasible and further improvements in the quality of the NEDF, and especially in its dynamics properties, can lead to a theoretical microscopic framework with great predictive power. TDSLDA will offer insights into nuclear processes which are either very difficult or even impossible to obtain in the laboratory.

4.8 Vortex pinning in neutron star crust.

A rather old puzzle in nuclear astrophysics is to explain why neutron 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 reaching $\approx 10^{14}g/cm^{3}$ and they rotate extremely fast, with periods as short as milliseconds. As they evolve they radiate energy and their rotation period increases very slowly. However, at practically random times they experience sudden accelerations and their periods increase very rapidly, practically instanteneously. Anderson and Itoh [81] suggested that since neutron stars are mainly superfluid neutron matter, there are many quantized vortices, which likely are pinned in the so called neutron star crust. Neutron stars are not composed entirely of neutrons, which are unstable against beta-decay, but also of a relatively small fraction of protons and electrons. Electrons are relativistic, very delocalized, and likely in a normal phase. But the neutrons and protons 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 and tubes. These “nuclei” from a Coulomb lattice in a superfluid neutron superfluid. The neutron superfluid under rotation is threaded with quantized vortices, which might attach to nuclei or interstitially. During the neutron star rotation these vortices can unpin and find a new position, which leads to a change in the star moment of inertia and therefore its period, or glitches, in analogy with flux creeping in hard superconducting magnets according to Anderson and Itoh [81].

For decades nuclear physicists and nuclear astrophysicists attempted to figure out whether a quantized vortex is pinned on a nucleus or whether it is repelled by a nucleus in the neutron star crust, and the results have been all over the map, from attraction to repulsion, with no clear indication even on the magnitude of this interaction. In 2013 we suggest a new approach to this problem, which could simply indicate whether the interaction is either repulsive or attractive. A nucleus scattering on a quantized vortex would be deflected either to right or to the left depending on the character of the force [82]. The size and sign of the nucleus vortex interaction arises in any theoretical approach as a difference between two very large quantiles, which cannot be computed with enough accuracy. In 2016 we have implemented this recipe [83], see also Fig. 10 and we were able to extract even more finer details of the nucleus-vortex interaction, namely how each small linear element of the vortex line interacts with the nucleus. During their interaction the vortex line is deformed and stretched. We have established not only that the interaction between a nucleus and a vortex is repulsive, for all possible densities on the neutron matter, but we have extracted with quite high accuracy the magnitude of this interaction.

4.9 Collisions of heavy nuclei and phase locking.

I will address only a single aspect of this rather wide field, which likely is of interest outside nuclear physics, the phase locking between two colliding superfluid droplets. Since the gauge symmetry is spontaneously broken in superfluids, it is reasonable to wonder under what conditions the relative phase of two superfluids is physically relevant. The Josephson effect [84, 85], experiments with cold Bose or Fermi atoms [86, 70, 72, 73, 87, 74], and the superfluid fragments emerging from nuclear fission [88], are just a few examples where that is the case.

Recently Magierski, et al. [89] reported on a surprising and strong dependence of the properties of the emerging fragments on the relative phase of the pairing condensates in the initial colliding nuclei. Related theoretical issues in superfluid helium have been discussed in literature by Anderson [90]. In Ref. [91] we have demonstrated that this is typical behavior of superfluids in the weak coupling limit. With increasing interaction strength however, the initially independent phases of the two order parameters in the colliding partners quickly become phase locked, as the strong coupling favors an overall phase rigidity of the entire condensate, and upon their separation the emerging superfluid fragments become entangled, see Fig. 11.

4.10 Fluctuations and dissipation.

One of the main difficulties with TDDFT is that in a time-dependent framework is likely to always lead to the same final state, even if one were to consider a range of different initial conditions[56]. As I discussed briefly above, at least in the case of nuclear fission and nuclear reactions in general that is definitely not the case. In nuclear fission one has a wide distribution of fission fragments masses and charges. Can one find within TDDFT a sensible solution to this problem?

The prevalent approach to perform real time evolution of many-nucleon systems was developed in the 1970’s using path-integral techniques; see [92] for a review. Starting with a system described by a many-body Hamiltonian $\hat{H}$, one performs a Hubbard-Stratonovich transformation on the many-nucleon evolution operator by introducing auxiliary one-body fields $\sigma(t)$, formally,

[TABLE]

where $\mathcal{D}[\sigma(t)]$ is an appropriate measure depending on all auxiliary fields, $W[\sigma(t)]$ is a Gaussian weight and $\hat{h}[\sigma(t]$ is a one-body Hamiltonian built with the auxiliary one-body fields $\sigma(t)$.

Using the stationary phase approximation, one selects a single mean field trajectory $\bar{\sigma}(t)$, which one may simulate with the TDDFT trajectory. If the initial state is a (generalized) Slater determinant, the final state is also a (generalized) Slater determinant under the evolution of this stationary phase mean field trajectory. After a trivial change of integration variables $\sigma(t)=\bar{\sigma}(t)+\eta(t)$, the true many-body wave function can be put into the form

[TABLE]

where $\Psi(0)$ is the initial wave function, and $\eta(t)$ are fluctuations around the stationary phase trajectory $\bar{\sigma}(t)$ at time $t$. In Eq. (35) 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 true many-nucleon wave function has a similar mathematical structure as the wave function in the time-dependent generator coordinate method (TDGCM) introduced by Wheeler et al. [93, 94]. One cannot but see the analogy in treating fluctuations around the mean field trajectory with the classical Langevin description of nuclear collective motion as well [95]. The representation (35) (which is an exact one) of the many-body wave function has the great advantage 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 function. One choice is a single (generalized) Slater determinant and another is a superposition of many such (generalized) Slater determinants.

In our fission studies we observed that one-body dissipation (i.e. the collisions of fermions with the moving surface of the nucleus) or Landau type of dissipation leads to a very quick dissipation of the f=collective flow energy into intrinsic/thermal of the fissioning nucleus [56]. The intrinsic motion of the descending nucleus from the outer saddle towards the scission configuration is similar to the downward motion of a heavy railway car on a very steep hill with its wheels blocked. The wheels dot not rotate but slip and become extremely hot, since almost the entire gravitational potential energy of the railway car at the top of the hill is converted into heat and very little of it is converted into collective kinetic energy. In this case, the railway car velocity is equal to the terminal velocity. An object attains a terminal velocity when the conservative force is balanced by the friction force, the acceleration of the system vanishes and the inertia plays no role in its dynamics. The motion of the railway car is strongly non-adiabatic. Similarly for a nucleus the flow of energy from the collective/shape degrees of freedom is controlled by the entropy of the intrinsic degrees of freedom. This energy transfer is to some extend stochastic, as is brownian motion, and can be simulated with TDDFT evolution equations augmented to incorporate dissipation and fluctuations we introduce have the form

[TABLE]

where $\hat{\bm{p}}=-i\hbar\bm{\nabla}$ (not to be confused with ${\bf p}({\bf r},t)$), the index $k$ runs over the neutron and proton quasi-particle states and where $\psi_{k}({\bf r},t)$ are 4-component quasi-particle wave functions and $h[n]$ is a $4\times 4$ partial differential operator [96, 88]. The fields ${\bf u}({\bf r},t)$ and $u_{0}({\bf r},t)$ generate both rotational and irrotational dynamics and the term proportional to $\gamma$ is a quantum frictional term. In the presence of this additional term alone $\dot{E}_{\text{tot}}\leq 0$, as in the case of the presence of a classical friction term. This is only a phenomenological solution, which likely will lead to reasonable results and might serve as inspiration.

5 Conclusions.

I presented an extension of the DFT to superfluid systems using a local pairing field and a further extension of this framework to time-dependent phenomena using the adiabatic approximation. This static extension SLDA can be correlated with ab initio calculations, it is strongly constrained by a number of theoretical arguments, and at least in the case of cold atoms it appears to have a quite good accuracy. The applications to nuclear phenomena is based to a large extent to a phenomenological energy density functional, which has a sub-percent accuracy for a large number of nuclei and their static properties. The time-dependent extension, which in addition is required to satisfy local Galilean covariance, appears to provide a correct description of many non-equilibrium processes in nuclei and neutron stars as well.

{acknowledgement}

I am grateful to my former students and many collaborators, whose names are in the papers we have published together and referenced here, and to many of my colleagues with whom I have discussed many of the issues covered in this review. The results presented here have been obtained with support by US DOE Grant No. DE-FG02-97ER-41014 and in part by NNSA cooperative agreement DE-NA0003841. Calculations have been performed at the OLCF Titan and Piz Daint and for generating initial configurations for direct input into the TDDFT code at OLCF Titan and NERSC Edison. This research used resources of the Oak Ridge Leadership Computing Facility, which is a U.S. DOE Office of Science User Facility supported under Contract No. DE- AC05-00OR22725 and of the National Energy Research Scientific computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We acknowledge PRACE for awarding us access to resource Piz Daint based at the Swiss National Supercomputing Centre (CSCS), decision No. 2016153479. This work is supported by ”High Performance Computing Infrastructure” in Japan, Project ID: hp180048. Some simulations were carried out on the Tsubame 3.0 supercomputer at Tokyo Institute of Technology.

Graphical Table of Contents

In normal fluids turbulence appears only when viscosity is sufficiently strong. While viscosity vanishes in superfluids at zero temperature turbulence is possible if quantum vortices exist. These vortices while moving thorough the fluid will cross and recombine, in a manner similar to DNA crossing and recombination, leading to a chaotic motion dubbed by Feyman as quantum turbulence.

GTOC image: