Expansion of the strongly interacting superfluid Fermi gas: symmetries and self-similar regimes

TL;DR

This paper analyzes the expansion dynamics of a strongly interacting superfluid Fermi gas at unitarity, revealing symmetries, self-similar regimes, and the connection to classical hydrodynamics and the virial theorem.

Contribution

It introduces a self-similar solution framework for the gas expansion, linking quantum symmetries with classical hydrodynamics and the virial theorem in the unitary regime.

Findings

Asymptotic expansion matches the quasi-classical self-similar solution.

Virial theorem governs the expansion dynamics.

Self-similar solutions describe shape deformations during expansion.

Abstract

We consider an expansion of the strongly interacting superfluid Fermi gas in a vacuum, assuming absence of the trapping potential, in the so-called unitary regime (see, for instance, \cite{pitaevskii2008superfluid}) when the chemical potential $\mu \propto \hbar^2n^{2/3}/m$ where $n$ is the density of the Bose-Einstein condensate of Cooper pairs of fermionic atoms. In low temperatures, $T\to 0$, such expansion can be described in the framework of the Gross-Pitaevskii equation (GPE). Because of the chemical potential dependence on the density, $\sim n^{2/3}$, the GPE has additional symmetries, resulting in the existence of the virial theorem \cite% {vlasov1971averaged}, connecting the mean size of the gas cloud and its Hamiltonian. It leads asymptotically at $t\to\infty$ to the gas cloud expansion, linearly growing in time. We study such asymptotics, and reveal the perfect match between the quasi-classical self-similar solution and the asymptotic expansion of the non-interacting gas. This match is governed by the virial theorem, derived through utilizing the Talanov transformation \cite{talanov1970focusing}, which was first obtained for the stationary self-focusing of light in media with a cubic nonlinearity due to the Kerr effect. In the quasi-classical limit, the equations of motion coincide with 3D hydrodynamics for the perfect monoatomic gas with $\gamma=5/3$. Their self-similar solution describes, on the background of the gas expansion, the angular deformities of the gas shape in the framework of the Ermakov--Ray--Reid type system.

Equations

$$m\frac{d^2}{d{t}^2}\int r^2|\psi {|}^{2}d\mathbf{r}=4H,$$

$$H=\int \left[\frac{{\hslash }^2}{2m}|\nabla \psi {|}^2+g|\psi {|}^4\right]d\mathbf{r}.$$

$$i\hslash \frac{\partial \psi }{\partial t}=-\frac{{\hslash }^2}{2(2m)}\Delta \psi +\mu (n)\psi ,$$

$$\mu (n)=2(1+\beta ){\epsilon }_F,$$

$${\epsilon }_F=\frac{{\hslash }^2}{2m}{\left(6{\pi }^{2}n\right)}^{2/3}.$$

$$i\frac{\partial \psi }{\partial t}=-\frac{1}{2}\Delta \psi +\mu (n)\psi$$

$$\mu (n)=2(1+\beta ){\left(6{\pi }^{2}n\right)}^{2/3}.$$

$$\frac{\partial n}{\partial t}+\left(\nabla \cdot n\nabla \phi \right)$$

$$\frac{\partial \phi }{\partial t}+\left[\mu (n)+\frac{{\left(\nabla \phi \right)}^2}{2}+T_{QP}\right]$$

$$T_{QP}=-\frac{\Delta \sqrt{n}}{2\sqrt{n}}.$$

$$2(1+\beta ){\left(6{\pi }^{2}n\right)}^{2/3}\to \frac{5}{3}{n}^{2/3}$$

$$i\frac{\partial \psi }{\partial t}+\frac{1}{2}\Delta \psi -(\upsilon +1)|\psi {|}^{2\upsilon }\psi =0,$$

$$i\frac{\partial \psi }{\partial t}=\frac{\delta H}{\delta {\psi }^{\ast }},$$

$$H=\int \left[\frac{1}{2}|\nabla \psi {|}^2+|\psi {|}^{2(\upsilon +1)}\right]d\mathbf{r},$$

$$\frac{\partial n}{\partial t}=\frac{\delta H}{\delta \phi },\frac{\partial \phi }{\partial t}=-\frac{\delta H}{\delta n}$$

$$H=\int \left[\frac{n{\left(\nabla \phi \right)}^2}{2}+\frac{{\left(\nabla \sqrt{n}\right)}^2}{2}+n^{\upsilon +1}\right]dr.$$

$$\int r^2|\psi {|}^{2}d\mathbf{r}=2H{t}^2+C_{1}t+C_2.$$

$$\int r^2|\psi {|}^{2}d\mathbf{r}\to 2H{t}^2.$$

$${⟨r^2⟩}^{1/2}\propto t\sqrt{2H/N}.$$

$$\frac{\partial \phi }{\partial t}+\frac{{\left(\nabla \phi \right)}^2}{2}+\frac{5}{3}{n}^{2/3}=0.$$

$$n=\frac{1}{a_x{a}_y{a}_z}f\left(\frac{x}{a_x},\frac{y}{a_y},\frac{z}{a_z}\right),$$

$$\frac{\partial n}{\partial t}=-\frac{1}{a_x{a}_y{a}_z}\sum _i\frac{\dot{a_i}}{a_i}\frac{\partial }{\partial {\xi }_i}(f{\xi }_i),$$

$${\xi }_x=\frac{x}{a_x},{\xi }_y=\frac{y}{a_y},{\xi }_z=\frac{z}{a_z}.$$

$$\phi ={\phi }_0(t)+\sum _l\frac{\dot{a_l}{a}_l}{2}{\xi }_l^2,$$

$$\frac{\partial \phi }{\partial t}+\frac{{\left(\nabla \phi \right)}^2}{2}=\frac{d{\phi }_0(t)}{dt}+\sum _l\frac{\ddot{a_l}{a}_l}{2}{\xi }_l^2.$$

$$\frac{5}{3}{n}^{2/3}=\frac{5}{3}\frac{1}{{\left(a_x{a}_y{a}_z\right)}^{2/3}}{\left[f\left({\xi }_x,{\xi }_y,{\xi }_z\right)\right]}^{2/3}.$$

$$\frac{d{\phi }_0(t)}{dt}=-\frac{5}{3}\frac{1}{{\left(a_x{a}_y{a}_z\right)}^{2/3}}f(0{)}^{2/3},$$

$$\sum _l\frac{\ddot{a_l}{a}_l}{2}{\xi }_l^2=\frac{5}{3{\left(a_x{a}_y{a}_z\right)}^{2/3}}[f{\left(0\right)}^{2/3}-f{\left(\xi \right)}^{2/3}].$$

$$\ddot{a_x}{a}_x=\ddot{a_y}{a}_y=\ddot{a_z}{a}_z=\frac{\lambda }{{\left(a_x{a}_y{a}_z\right)}^{2/3}}$$

$$f\left(\xi \right)={\left[f{\left(0\right)}^{2/3}-\frac{3\lambda }{10}{\xi }^2\right]}^{3/2},$$

$$n=\frac{1}{a_x{a}_y{a}_z}{\left[1-\frac{3\lambda }{10}{\xi }^2\right]}^{3/2}$$

Full text

I Introduction

Since the discovery of Bose-Einstein condensation (BEC) in alkali gases of bosonic isotopes ${}^{7}Li,^{23}Na,^{87}Rb$ [4, 5, 6] the time of flight experiments connected with the expansion of the Bose condensate cloud from the trap, when the trapping potential is switched off, served as one of the important proves of the superfluid transition in the gas. Note that in these experiments the velocity distribution in the gas expansion has a typical bimodal shape which corresponds to two components - a normal and a superfluid ones. The velocity distribution of the normal component has a thermal (Maxwell-type) shape while the velocity of the superfluid component is governed only by the interaction parameter of the weakly non-ideal (in the Gross-Pitaevskii meaning) Bose gas and the total number of particles.

The first scaling time-dependent solutions both for the condensate and thermal gas in the hydrodynamic regime for the anisotropic trap was obtained by Yu.Kagan, Surkov and Shlyapnikov [7]. In particular, in [7] the spectrum of breathing modes of the oscillating type in the trapping potential was also determined. Later on the self-similar regimes were observed in the experiments of the Thomas’ group [8] for the anisotropic expansion of strongly interacting degenerate Fermi gas of atoms ${}^{6}Li$ from the optical trap. The measurements of this group were performed in the regime of Feshbach resonance [9, 10] which corresponds to the BCS-BEC crossover [11, 12, 13] between extended Cooper pairs of fermionic atoms (BCS) and BEC of tightly bound fermion pairs (molecules or dimers ${}^{6}Li_{2}$).

Note that while exploiting Feshbach resonance it is possible to reach large absolute values of the $s$-wave scattering length $|a_{s}|$ and thus to increase sharply the critical temperature of the superfluid transition for the fixed particle density $n$ in the trap reaching experimentally accessible values of critical temperature $T_{c}$. Remind that BCS phase of extended pairs on the phase diagram of the BCS-BEC crossover corresponds to the positive values of the chemical potential $\mu>0$ and negative values of the scattering length $a_{s}<0$. In the same time the BEC phase of local pairs vice versa corresponds to a positive scattering length $a_{s}>0$ and negative chemical potential $\mu<0$. For small positive values of the gas parameter $0<a_{s}k_{F}\ll 1$ (where $p_{F}=\hbar k_{F}$ is Fermi momentum) we are in the dilute BEC domain which describes a Bogoliubov gas of weakly repulsive composed bosons (local fermion pairs). In the same time, small negative values of the gas parameter $-1\ll(a_{s}k_{F})^{-1}<0$ corresponds to the dilute BCS regime describing weakly attractive degenerate Fermi gas of atoms.

In the regime of Feshbach resonance it is convenient to consider the phase diagram of the BCS-BEC crossover in terms of the dimensionless temperature $T/\varepsilon_{F}$ (the vertical axis in Fig. 1) and the inverse gas parameter $(a_{s}k_{F})^{-1}$ which scales linearly with the shift of external magnetic field $\Delta B=B-B_{0}$ from the resonance value $B_{0}$ at the horizontal axis (see [13]). Then on the phase diagram between dilute BCS and BEC domains there appears an intermediate region of strong correlations where the inverse gas parameter varies in the interval $-1<(ak_{F})^{-1}<1$ . In this region it is difficult to develop rigorous diagrammatic expansions and only some reasonable (conserving) approximations like the self-consistent $T$-matrix approximations [13] are available.

However, we have one special point in this regime which lies in the middle of strongly interacting region and corresponds to the unitarian limit $(a_{s}k_{F})^{-1}\rightarrow 0$ where exact (and universal) results can be obtained (see e.g. [1] for the review). The reason for the universality of the unitarian limit is connected with the fact that there is no other energy scales besides Fermi energy $\varepsilon_{F}$ in this point. Thus, both the chemical potential $\mu$ and the critical temperature $T_{c}$ in the unitarian limit scale linearly with $\varepsilon_{F}$.

Note also that at the middle of 1990-es in the field of ultra-cold quantum gases only anisotropic traps were experimentally available. Anisotropy degree for the gas clouds distributions in traps changed from 3D (with weak anisotropy) up to quasi-1D (with strong anisotropy for cigar-shaped traps). For instance, the cigar-shaped traps were considered in the Thomas’ group experiments [8]. Almost spherical trap, to our knowledge, appeared only in 2015 (see [15]). For the spherical traps all three trapping frequencies are equal, $\omega_{x}=\omega_{y}=\omega_{z}$. For the cigar-shaped traps we have the following hierarchy of the trapping frequencies: $\omega_{z}\ll\{\omega_{x},\omega_{y}\}$.

Quite recently also the disk-shaped traps with quasi-2D cloud distributions become experimentally available (see [16] and references therein). Parameters of these traps correspond to inverse hierarchy of the trapping frequencies with respect to cigar-shaped traps, namely $\omega_{z}\gg\{\omega_{x},\omega_{y}\}$. In this paper, we are discussing also expansion of the superfluid Fermi gas into vacuum from disk-shaped traps [17]. It is interesting to note that the simplest ballistic picture of the Fermi gas expansion for the disk-shaped superfluid in the unitarian regime considered in [17] coincides with the exact quasi-classical self-similar solution first found by Anisimov and Lysikov for the hydrodynamic expansion of the classical gas with adiabatic constant $\gamma=5/3$ [18] (in the absence of vorticity).

In this paper, we consider an expansion of the strongly interacting superfluid Fermi gas in a vacuum in the unitary regime when the chemical potential $\mu\propto\hbar^{2}n^{2/3}/m$ where $n$ is the density of the Bose-Einstein condensate of Cooper pairs of fermionic atoms assuming temperature $T\to 0$. Such expansion can be described in the framework of the Gross-Pitaevskii equation (GPE). Because of the chemical potential dependence on the density $\sim n^{2/3}$ the GPE has additional symmetries resulting in existence of the virial theorem [2] connected the mean size of the gas cloud and its Hamiltonian. It leads asymptotically at $t\rightarrow\infty$ to the linear in time expansion of the gas. We carefully study such asymptotics and reveal a perfect matching between the quasi-classical self-similar solution and the asymptotic expansion of the non-interacting gas.

The paper is organized as follows. In the next Section we discuss the problem concerning symmetries of the GPE in the unitarian limit and how these symmetries are connected with those found for the nonlinear Schrodinger equation (NLSE) in the critical case when the virial theorem can be applied for description of the stationary self-focusing of light in media with the Kerr nonlinearity [2] and symmetries in the case of ideal monoatomic gas. Section 3 mainly deals with self-similar solution of the anisotropic type of the quasi-classical GPE in the unitarian limit for expansion of the Fermi superfluid gas. This solution describes the angular deformations of the gas shape on the background of the gas expansion. In Section 4, we discuss in which extent the analytical results obtained in the previous sections are related with experimental data. Conclusion summarizes all results of this paper.

II Symmetries and integrals of motion

First of all, we would like to remind that the topic of gas expansion was very popular in the hydrodynamic content in 60-s of the XXth century. The first classical works were performed by L.V. Ovsyannikov (1956) [19] and F.J. Dyson (1968) [20]. These studies had a lot of applications not only in hydrodynamics but also in astrophysics (see, e.g. the original paper by Ya. B. Zel’dovich [21]).

In 1970 S.I. Anisimov and Yu.I. Lysikov [18] discovered very interesting phenomenon connected with the nonlinear angular deformation of the gas cloud while its expansion. Such behavior directly follows from their remarkable solution for a gas with specific heat ratio $5/3$ (see e.g. reviews [22, 23] and references therein). This result, as was pointed out by I.E. Dzyaloshinskii (private communication, 1970), represents a consequence of the symmetry which is well known in quantum mechanics for motion of a non-relativistic particle in the potential $V(r)=\beta/r^{2}$. This symmetry, independent on the sign of $\beta$, is dilatations of both spatial coordinates and time for which $\mathbf{r}\rightarrow\alpha\mathbf{r}$ and $t\rightarrow\alpha^{2}t$ where $\alpha$ is a scaling parameter. Indeed, such symmetry first time was exploited by V.P. Ermakov in 1880 [24] to construct solutions for some mechanical systems including motion of a particle in the potential which is a combination of the oscillator potential and $V(r)=\beta/r^{2}$. In seventies of the XXth century, Ray and Reid [25] rediscovered the Ermakov results. Now all such equations are accepted to call the Ermakov-Ray-Reid systems (see, e.g. [26] and references therein). As we will show in this paper, this additional symmetry for the GPE takes place for both attractive and repulsive interactions (in optical content, corresponding to focusing and defocusing nonlinearities). Note that, in quantum mechanics (see [27]), for the attractive potential, $V<0$, with constant $|\beta|$ larger some critical value ($=\hbar^{2}/\left(8m\right)$) the quantum falling of a particle with mass $m$ into the center is possible which can be understood as collapse. Moreover, this falling becomes more quasiclassical while approach the center. In the case of the Gross-Pitaevskii equation (GPE) [28, 29] which can be applied for description of nonlinear dynamics of the Bose condensate for diluted gases, the kinetic energy has the same scaling as in the usual quantum mechanics, i.e. $\propto\alpha^{-2}$. The nonlinear interaction term in the GPE, due to the $s$-scattering, has a scaling $\propto\alpha^{-d}$ which appears from the conservation of the total number of particles $N=\int|\psi|^{2}d\mathbf{r}$ with $d$ being the space dimension and $\psi$ the wave function of the Bose condensate. Thus, at $d=2$ only we have the situation analogous to that in the quantum mechanics for potentials $V(r)=\beta/r^{2}$. This case, as it was first time demonstrated by Vlasov, Petrishchev and Talanov [2] for the 2D nonlinear Schrodinger equation, is very special for which the so-called virial relation is valid:

[TABLE]

where the Hamiltonian $H$ in the case of the GPE for the Bose condensate has the form

[TABLE]

Here the coupling coefficient $g=4\pi\hbar^{2}a_{s}/m$ with $a_{s}$ being scattering length and $m$ particle mass. It is necessary to emphasize that the virial relation (1) is valid for any sign of $g$. The only restriction follows from the requirement of convergence of the integrals standing in (1). It is worth noting that in classical mechanics the virial theorem establishes the ratio between mean values of the total kinetic and potential energies. The simplest way to derive this theorem is calculation of the second time derivative of a moment of inertia (this results in the virial relation like Eq. (1)) and then averaging it in time. Further we will call relation (1) as the virial theorem.

In this paper we consider another example of the same symmetry, when the generalized Gross-Pitaevskii equation [1] can be applied for description of the strongly interacting Fermi gas in the superfluid phase (at $T=0$):

[TABLE]

where $\psi$ is the wave function of the Bose condensate of fermion pairs, $m$ is a fermion mass ($2m$ is a mass of a fermion pair), $\mu$ is the chemical potential. In the unitarian limit (when $(k_{F}a)^{-1}\rightarrow 0$) the chemical potential reads (see e.g. [1]):

[TABLE]

where universal interaction parameter $\beta=-0.63$, in accordance with [30], [31, 32], [33], and local Fermi energy

[TABLE]

Here $n=|\psi|^{2}$ is concentration of fermionic pairs. Below we will normalize density $n$ by its initial maximum value $n_{0}$, inverse time $t^{-1}$ by $\frac{\hbar}{2m}n_{0}^{2/3}$ and coordinate $r$ by $n_{0}^{-1/3}$. In these new (dimensionless) units equation (2) reads as

[TABLE]

where

[TABLE]

Choosing the standard ansatz for the $\psi$-function, $\psi=\sqrt{n(r,t)}\exp\left(i\varphi(r,t)\right)$, and separating then real and imaginary parts in (2) we get the system of continuity and Euler (eiconal) equations

[TABLE]

where $\mathbf{v=}\nabla\varphi$ has a meaning of velocity. Here we used the condition of the absence of vortices $\nabla\times\mathbf{v}=0$.

Note that the term in Eq. (7) represents the quantum pressure given by

[TABLE]

Throughout the main part of the present paper we will neglect this term and will discuss its possible role in two last Sections. Neglecting quantum pressure corresponds to the quasi-classical (or eikonal ) approximation (called also the time-dependent Thomas-Fermi approximation) which assumes more rapid space and time variations of phase (larger phase gradients and time-derivatives) in comparison with the space and time variations of the modulus of the $\psi$-function in Eq.(4).

It is important to emphasize that the generalized Gross-Pitaevskii equation (4) coincides with the nonlinear Schroedinger equation (NLSE) widely used in nonlinear optics and plasma physics. It is convenient to exclude in (5) the factor $2(1+\beta)\left(6\pi^{2}\right)^{2/3}$ by simple rescaling of the density $n$,

[TABLE]

so that equation (4) takes the standard form accepted for the NLSE,

[TABLE]

with the exponent $\upsilon=2/3$. This equation can be written in the Hamiltonian form

[TABLE]

where Hamiltonian

[TABLE]

with the first term coinciding with the total kinetic energy and the second one responsible for nonlinear interaction of the repulsion type. After applying the transformation $\psi=\sqrt{n(r,t)}\exp\left(i\varphi(r,t)\right)$ equations for density $n$ and phase $\varphi$ remain the Hamiltonian form,

[TABLE]

where the Hamiltonian coincides with (10). In terms of $n$ and $\varphi$, $H$ takes the form

[TABLE]

The Hamiltonian equations of motion (11) are the same Eqs. (6) and (7) transformed under simple rescaling; $n$ and $\varphi$ in this case play the role of canonically conjugated quantities.

The second term in $H$ is responsible for the quantum pressure in Eq. (7). In the quasiclassical limit (the Thomas-Fermi approximation), this term becomes small and can be neglected so that we arrive at the hydrodynamic equations for potential flow of monoatomic gas with specific heat ratio (adiabatic index) $\gamma=5/3$ ( $\nu=2/3$). This $\gamma$ is remarkable for both NLSE and its quasiclassical limit. It turns out that the equations of motion in this case have two additional symmetries. The first symmetry forms dilatation group of the scaling type: $\mathbf{r\rightarrow}\alpha\mathbf{r}$ and $t\mathbf{\rightarrow}\alpha^{2}t$ . In usual quantum mechanics, this symmetry appears for the potential $V(r)\sim r^{-2}$ independently on both the potential sign (attraction or repulsion) and space dimension $d$. However, for the NLSE (9) such symmetry appears as a result of the conservation of the total number of particles $N=\int|\psi|^{2}d\mathbf{r}$ so that at $d=3$ only the nonlinear potential $\sim|\psi|^{4/3}$ in Eq. (9) has the same scaling as the Laplace operator $\Delta$. At $d=2$ such symmetry takes place already for the nonlinear potential $\sim|\psi|^{2}$ (in this case the NLSE describes the stationary self-focusing of light in a media with the Kerr nonlinearity). In the general case the dilatation symmetry arises at $\nu=2/d$ (see, for instance, [34, 35]). The second symmetry of the conformal type first time was found by V.I. Talanov for the cubic NLSE at $d=2$ [3] (1970) and called now the Talanov transformations. In optical content these are the lens transformations well known in a linear optics.

These symmetries are of the Noether type and generate two additional integrals of motion. They can be obtained from the virial theorem (1) (first time obtained for the 2D cubic NLSE in [2]), after twice integration in time (dimensionless variables):

[TABLE]

Hence we get asymptotically at $t\rightarrow\infty$, independently on $C_{1}$ and $C_{2}$,

[TABLE]

Therefore the mean size (indeed, r.m.s) of the gas cloud varies at large $t$ linearly in time,

[TABLE]

This result is very important also since it perfectly matches quasi-classical solution in Eq.(13) with the linearly varying solutions for the non-interacting particles (ballistic expansion).

It should be emphasized that the virial theorem (1) is the exact result, it can be applied in particular in the quasi-classical limit also when the quantum pressure term in $H$ is eliminated. The latter corresponds to the classical gas expansion with $\gamma=5/3$. In this case Anisimov and Lysikov (1970) [18] constructed exact anisotropic self-similar solution based, in fact, on existence of two integrals of motion $C_{1}$ and $C_{2}$ (see the next subsection). This solution describes the gas expansion in time in correspondence with (12) with nonlinear angular deformation of the gas shape.

III Self-similar quasi-classical solution

Note that the corresponding system of gas dynamics equations in the quasi-classical limit is described by Eqs. (6,7,9,10) by neglecting quantum pressure in (7) and with the eikonal equation for the phase which in the unitarian limit (for $\upsilon+1=5/3$) is given by

[TABLE]

Let us search for a solution of these equations in the self-similar form (see e.g. [18, 36, 37, 38]),

[TABLE]

assuming that three scaling parameters $a_{x},a_{y},a_{z}$ are functions of time. Note that the ansatz (15) conserves the total number of particles.

Calculating the (partial) time derivative from density we get

[TABLE]

where $i=(x,y,z)$. Here we introduced convenient notations for the self-similar variables

[TABLE]

Then the continuity equation admits integration resulting in relations for the phase $\varphi$

[TABLE]

where function $\varphi_{0}(t)$ can be found after substitution in the eikonal equation. Calculating first two terms in (14) yields:

[TABLE]

The third term in (14) reads

[TABLE]

Hence in order to satisfy the self-similar ansatz one needs to require that

[TABLE]

Hence we conclude that

[TABLE]

where $\lambda$ is arbitrary positive constant. For $f\left(\mathbf{\xi}\right)$ we have

[TABLE]

where further we will put $f\left(0\right)=1$. Respectively, the density is written as

[TABLE]

here the constant $\lambda$ will be found from the initial condition.

We will assume that initially the density distribution is defined from the Thomas-Fermi approximation. In the presence of harmonic trap, at the stationary state we have the equilibrium condition

[TABLE]

where $\mu(n)$ is defined by (5). Remind that because of pairing in this expression $2m$ stands instead of $m$. This gives the initial density distribution:

[TABLE]

where $\omega_{m}=\max(\omega_{i})$, $a_{i}(0)=\omega_{m}/\omega_{i}$. This profile matches precisely with self-similar solution (19) at $t=0$. Hence we have that

[TABLE]

or in terms of $N$

[TABLE]

Function $f\left(\mathbf{\xi}\right)$ (18) is spherically symmetric with respect to $\mathbf{\xi}$ . It varies from $1$ at $\xi=0$ up to zero at $\xi_{\max}=$ $\sqrt{10/3\lambda}$; above $\xi_{\max}$ the density $n$ is equal to zero (see Fig.2).

In accordance with (17) dynamics of three scaling parameters $a_{i}(t)$ ($i=1,2,3$) is described by the Newton equations for motion of a particle

[TABLE]

where potential

[TABLE]

It is worth noting that at the point $\xi=\xi_{\max}$ the obtained quasiclassical solution given by (15-20) breaks down that follows from estimation of the quantum pressure term which becomes infinitely large. In this case, $\xi=\xi_{\max}$ plays the role of a reflection point in the usual quasiclassical approximation in quantum mechanics. This means that at the vicinity $\Delta\xi$ around $\xi=\xi_{\max}$ one needs to match the constructed solution at $\xi<\xi_{\max}$ (inner region) with that at $\xi>\xi_{\max}$ (outer region) where we should neglect nonlinearity in the NLSE (free Schrodinger equation). This problem was discussed in details in [36] for the strong collapse regime in the supercritical NLSE with $d=3$ and $\nu=1$. In the given case, the matching problem can be also resolved if $\Delta\xi\ll\xi_{\max}$ when the matching solution is expressed in terms of the Painleve function. It should be noted that for the NLSE (9) with account of $T_{QP}$ the scaling remains the same as for the quasiclassical solution and by this reason $\Delta\xi$ can be considered as time-independent quantity so that the ratio between $\Delta\xi$ and $\xi_{\max}$ remains in time the same, unlike [36] where this ratio in 3D collapse regime is time-dependent and vanishes while approaching the collapse time.

We should notice also that near $\xi=\xi_{\max}$ the unitarian limit is not also applicable because $k_{F}$ becomes infinitely large. Near this point, however, the nonlinear term is small and we again return to the same matching problem like in the previous case.

III.1 The virial theorem for the scaling parameters

It is more or less evident that Newton equations (20) have to have the same symmetry properties as the original NLSE (16) what can be easily verified. Note first that for Eqs. (20) the energy integral is written in the standard form

[TABLE]

Secondly, for Eqs. (20) by direct calculation it is possible to get the virial theorem (1), written in terms of $a_{i}$. For $\sum a_{i}^{2}$ we have

[TABLE]

Then substitution of (20) into this relation gives finally

[TABLE]

that coincides with the virial identity (1). Its twice integration gives two constants $C_{1}$ and $C_{2}$ (integrals of motion):

[TABLE]

Hence

[TABLE]

In the isotropic (spherically symmetric) case when $a_{x}=a_{y}=a_{z}\equiv a$ the equations of motion transforms into one equation

[TABLE]

with the energy $E=\frac{3}{2}\left(\dot{a}^{2}+\frac{\lambda}{a^{2}}\right)$ and $3a^{2}=2Et^{2}+C_{1}t+C_{2}$. From the second relation we immediately have that gas cloud expands in radial direction asymptotically at $t\rightarrow\infty$ with constant velocity

[TABLE]

(ballistic regime). This result is in agreement with the virial theorem (13).

If we change the sign of the potential in Eq.(25) then we get the falling of the particle on the potential center which, as known in quantum mechanics, becomes more quasiclassical while approaching the center (see [27]).

For the expansion of a noninteracting gas from a harmonic potential

[TABLE]

and $v_{\infty}=const$ in agreement with our intuitive considerations and with Eq.(26) as well. (Let us remind that for quasi-2D disk-shaped traps the trapping frequency $\omega_{z}\gg\omega_{x}\simeq\omega_{y}$). Thus, we have almost perfect matching of ballistic results for non-interacting gas and quasi-classical results derived for strongly-interacting Fermi gas in the eikonal approximation.

It is worth noting that in the virial relation (22) besides total energy $E$ there enter two more constants (integrals of motion) $C_{1},C_{2}$ . In principle, if $C_{1}>0$ then the solution with $a^{2}\propto C_{1}t$ is possible for some intermediate times, but not initially. The regime with $a\propto\left(t_{0}-t\right)^{1/2}$ is typical for weak self-similar collapse (see [36, 39]).

III.2 Anisotropic self-similar solution

The simplest anisotropic case corresponds to the cylindrically symmetric expansion and is governed by the scaling parameters $a_{x}=a_{y}=a/\sqrt{2}$, $a_{z}=b$ . For $a>>b$ we have the case of an initially disc-shaped cloud while for $b\gg a$ we are effectively in the cigar-shape limit. An isotropic limit obviously corresponds to $b=a/\sqrt{2}$.

In the anisotropic cylindrically symmetric case Eqs.(17) read

[TABLE]

Note that in the initial moment of expansion when the trapping potential is switched off

[TABLE]

The effective potential in accordance with Eq. (21) is given by

[TABLE]

The corresponding Newton equations in agreement with (27) acquire the form

[TABLE]

Note that this system belongs to the so called Ermakov type of equations [24]. These equations describe the motion of two degrees of freedom and therefore to integrate this system it is enough one to have two autonomous integrals of motion which should be in involution. In our case, however, we have three integrals of motion. The first one is the total energy,

[TABLE]

The second and the third integrals are two constants $C_{1},C_{2}$ which appear in (12) while the double integration over time of the virial identity (1),

[TABLE]

The integrals (23), however, are not autonomous, they contain explicit dependence of time, and therefore can not provide a complete integration of the system. As we will see only their combination defines the needed integral of motion for the Ermakov type of equations.

Let us introduce now the polar coordinates for $a$ and $b$

[TABLE]

In these variables the virial theorem (29) acquires the evident form

[TABLE]

where the total energy $E$ in accordance with (28) is

[TABLE]

and correspondingly

[TABLE]

Multiplying now (30) by $r^{2}$ and using the relations (31) simple calculations give that the combination

[TABLE]

is a constant (the Ermakov integral). As the result, we arrive at conservation law for new ”energy”

[TABLE]

with new time $\tau$

[TABLE]

where

[TABLE]

plays a role of potential energy. It is always positive and goes to infinity for $\Phi\rightarrow 0$ and $\Phi\rightarrow\pi/2$.  The minimum of $U_{eff}(\Phi)=9\lambda/2$ corresponds to the isotropic case when $\sin\Phi_{\min}=1/\sqrt{3}$ . Graphically the effective potential $U_{eff}(\Phi)$ is shown on Fig.3.

The new time $\tau$  (33) can be easily expressed  through $t$ ,

[TABLE]

where $\chi^{2}=\widetilde{E}/E$ and $t_{0}=\frac{C_{1}}{4E}$ so that $\tau=0$ at $t=0$. If the initial velocity is equal zero (that is typical for experiment) the constant $C_{1}=0$ and

[TABLE]

In this case, asymptotically at $t\rightarrow\infty$

[TABLE]

The trajectory $\Phi(\tau)$ is defined from integration of (32)

[TABLE]

Hence the $\tau$-period of the oscillations in the potential $U_{eff}(\Phi)$ (34) is expressed through the integral

[TABLE]

where $\Phi^{(\pm)}$ are roots of equation $\widetilde{E}=U_{eff}(\Phi)$ (reflection points). This integral is expressed via elliptic integrals of the third order (see [18]). At large value of $\widetilde{E}$ oscillations are almost independent on the details of $U_{eff}(\Phi)$. Asymptotically in this case the angular velocity $\frac{d\Phi}{d\tau}\rightarrow$ $\pm\sqrt{2\widetilde{E}}$ and the $\tau$-period

[TABLE]

namely, in this limit $T$ exceeds in two times $\tau_{\infty}$ (35). Notice also that dependence of $T$ with respect to $\widetilde{E}$ is monotonic for the given potential $U_{eff}(\Phi)$ with the maximum corresponding to the potential minimum. This means that in the real experiment (what we will discuss in the next Section) in the better case it is possible one to observe only half of such oscillation, $t_{osc}$. Important, that a recurrence to the initial shape is impossible in this case. In the quasi-classical regime, the gas shape behavior will be different for cigar and disk initial conditions. For example, in the cigar case we start at fixed $\widetilde{E}$ from the left reflection point of the potential $U_{eff}(\Phi)$, in the disk case – from the right reflection point. Therefore the shape forms will coincide only for intermediate moments of time, far from the initial reflection points. We should take into account that at fixed $\widetilde{E}$ starting from any reflection point we can not reach its opposite reflection point.

It should be emphasized that the solution presented here was obtained first time by Anisimov and Lysikov [18] for expansion of ideal gas with $\gamma=5/3.$

III.3 The general anisotropic case

In the general anisotropic case, when all the scaling parameters are different $a_{x}\neq a_{y}\neq a_{z}$ it is convenient to introduce the spherical coordinates ( $r,\theta,\varphi$) where the total energy acquires the form

[TABLE]

Correspondingly introducing again Ermakov reduced energy $\widetilde{E}$ , being a sequence of the dilatation symmetry, and new time $\tau$, according the same prescriptions as in the preceding subsection, we get

[TABLE]

where the effective potential is now

[TABLE]

Thus, we arrive to the system for two degrees of freedom. As it was pointed out in the previous subsection the integral (36) is a consequence of the scaling symmetry, but for integration of the system it is not enough. As it was shown by Gaffet [40], this system indeed has one additional integral (besides $\widetilde{E}$) which follows from the Painleve test. Existence of these two integrals of motion guarantees complete integration of this system. As in the previous limit motion in potential (37) remains its nonlinear quasi-oscillation character.

IV Discussion of experimental data and

comparison with obtained results

The self-similar expansion of a strongly interacting Fermi gas from a cigar-shaped trap was observed in [8]. The images of the expanding gas are shown in Fig. 4(a). The transverse size grows rapidly, while the longitudinal size is nearly stationary, with a weak growth. In Fig. 4(b) one may see qualitative agreement between the time behavior of the gas expending shape and the self-similar-expansion model represented by Eqs. (15), (27).

The cloud images are changing on Fig. 4(a) from almost ellipsoid significantly stretched along $z$-axis (exposition $t=100$ $\mu s$), later on to the almost spherical shape (at $t=600$ $\mu s$) and, finally, from the spherical shape to the ellipsoid stretched now in the direction perpendicular to $z$. The total time of the observation was 2000 $\mu s$ which can be taken as a half period (or less) of the period of the angular shape oscillations, $t\leq t_{osc}/2$, in accordance with the results of the previous section. The frequency ratio (and thus the anisotropy ratio up to factor $\sqrt{2}$) in the experiments [8] was rather large initially (around 30) that follows from the Thomas-Fermi estimation.

Small deviation of the data from the self-similar behavior in Fig. 4(b) has been attributed [41] to the contribution of quantum pressure (8) into hydrodynamic model (6), (7). That $T_{QP}$ term is neglected to obtain self-similar solution of Eqs. (15), (27). This difference, however, may be explained even without the quantum pressure, by possible deviation of the equation of state from the $\mu\propto n^{2/3}$ dependence (3) since in experiment the interaction parameter is not tuned exactly on resonance $1/(k_{F}a_{s})=0$, with the estimate $1/(k_{F}a_{s})\simeq-0.14$ [8].

Exactly on resonance, the mean squared cloud size $\langle\mathbf{r}^{2}\rangle\equiv\langle x^{2}\rangle+\langle y^{2}\rangle+\langle z^{2}\rangle$ is found [42] to evolve as

[TABLE]

where $U(\mathbf{r})$ is the initial trapping potential. Expansion law (38) was obtained within the Thomas-Fermi approximation and coincides with the quasi-classical dependence for $\langle\mathbf{r}^{2}\rangle$ (22) in the unitarian limit. It should be emphasized that, according to (22), $\langle\mathbf{r}^{2}\rangle$ indeed depends linearly on energy $E$ that was verified in experiments [42].

In the case when the system is far from the unitarian point $(k_{F}a_{s})^{-1}=0$ experiments nevertheless give the parabolic time dependence for $\langle\mathbf{r}^{2}\rangle$. Small deviation of the data from the self-similar behavior in Fig. 4(b) has been attributed [41] to the contribution of quantum pressure (8) into hydrodynamic model (6), (14). That $T_{QP}$ term is neglected to obtain self-similar solution of Eqs. (15), (27). This difference, however, may be explained even without the quantum pressure, by reasons which include deviation of the equation of state from the $\mu\propto n^{2/3}$ dependence (3) since in experiment the interaction parameter is not tuned exactly on resonance $1/(k_{F}a)=0$, with the estimate $1/(k_{F}a)\simeq 1/(k_{\text{FI}}a)=-0.14$, where $k_{\text{FI}}$ is the Fermi wave vector of a noninteracting Fermi gas with the same atom number and in the same trap [8].

Expansion law (38) is indeed the same as for the ideal gas and coincides with the virial theorem (1) . Note that the expansion law is obtained for both zero- and finite-temperature gas with equation of state

[TABLE]

where $P$ is the pressure and $\mathcal{E}$ is the energy density, while the relation $\mu\propto n^{2/3}$ is a particular case of (39) that corresponds to the isentropic regime for $\gamma=5/3$. Equation of state (39) and expansion law (38) are consequences of resonant interaction with $1/(k_{F}a)=0$. Away from resonance, the expansion laws slightly differ from each other, which is seen in the measurements displayed in Fig. 5. Nevertheless, these laws for different parameters $1/(k_{F}a)$ have the same parabolic dependence on $t$. It should be emphasized that during the expansion the interaction parameter $1/(k_{F}a_{s})$ changes due to a drop in gas density. When the parameter values fall outside the interval $(-1.1)$, the quantum effects become less significant, and the gas expansion approaches the law for a classical monoatomic gas, which coincides, however, with (38). By these reasons we guess that expansions with $1/(k_{F}a)\simeq 1/(k_{\text{FI}}a_{s})=0.59$ and $1/(k_{F}a_{s})\simeq 1/(k_{\text{FI}}a_{s})=-0.61$ (red and blue curves in Fig. 5) correspond to normal Fermi gas.

V Conclusion

We have demonstrated that symmetry for the GPE in the unitarian limit, describing strongly interacting superfluid Fermi gas, provides existence of the virial theorem ((1)). As its consequence, independently on the ratio between quantum pressure and chemical potential while the Fermi superfluid gas expansion the size of the gas cloud scales linearly with time asymptotically so that the expansion velocity tends to the constant value, $v_{\infty}=\left(2H/N\right)^{1/2}$.

For description of the expansion of the strongly interacting superfluid Fermi gas we have applied the self-similar quasiclassical theory. For large time scales the theory matches quite well with simple ballistic ansatz and also with the initial quasi-classical distribution of trapping gas. This self-similar solution is a consequence of the scaling symmetry of the Ermakov type. At large times the size of the gas clouds scales linearly with time that is a consequence of the virial theorem. In the unitary limit, when both kinetic and potential energy scale linearly with the Fermi energy, our quasiclassical solution for superfluid quantum gas coincides with the Anisimov-Lysikov solution [18] for classical gas expansion in the isentropic regime. This anisotropic solution describes the nonlinear deformations of the cloud shape while self-similar gas expansion. For the initial condition in the cigar-shape form this solution demonstrates successively all the stages of gas expansion, starting from the distribution extended along the cigar axis, bypassing the spherically symmetrical one and ending with the distribution, turned at angle $\pi/2$ with respect to the initial cigar form. Such behavior was observed first time in experiments [8]. For the initial distribution in the form of a quasi-2D disk, all stages of expansion are inverse to those for the initial distribution in the cigar form.

In order to understand the role of the quantum pressure while the Fermi gas expansion we would like to note that the GPE (2) admits also the following self-similar substitution,

[TABLE]

where $\nu$ plays the role of the nonlinear eigenvalue for the differential equation for the function $F\left(\mathbf{\xi}\right)$ which is assumed to vanish at large $\mathbf{\xi=r/}\sqrt{t}$. For this substitution the nonlinear interaction term in (2) is of the same order as the quantum pressure one. However, this ansatz contradicts to the relation (12) following from the virial theorem. If one substitutes (40) in (12) it follows immediately that the quantity $\left\langle r^{2}\right\rangle$ grows linearly in time but $\left\langle r^{2}\right\rangle$ must grow quadratically as $t\rightarrow\infty$ because the Hamiltonian $H$ is strictly positive. This means that (40) can not be applied for description of the system for the whole space. We may hope only that this ansatz can be used in some region as, for instance it happens for the cubic NLSE describing weak collapse regime [36]. By this reason, we can suppose that account of the quantum pressure should provide a transition to the quasi-classical self-similar asymptotics because in the general case $\left\langle r^{2}\right\rangle$ behaves similarly in the quasi-classical limit, compare with (22). The verification of this assumption will be a subject of our future numerical study.

V.1 Acknowledgments

The work of E.A.K. was supported by the Russian Science Foundation (grant 19-72-30028), M.Yu.K is grateful for the support to the Russian Science Foundation under grant 18-12-00002, A. V. T. was supported by the Presidium of the Russian Academy of Sciences (Program “Newest methods of mathematical modeling of nonlinear dynamic systems”).