An Effective Field Description for Fermionic Superfluids

Abstract

In this chapter, we present the details of the derivation of an effective field theory (EFT) for a Fermi gas of neutral dilute atoms and apply it to study the structure of both vortices and solitons in superfluid Fermi gases throughout the BEC-BCS crossover. One of the merits of the effective field theory is that, for both applications, it can provide some form of analytical results. For one-dimensional solitons, the entire structure can be determined analytically, allowing for an easy analysis of soliton properties and dynamics across the BEC-BCS interaction domain. For vortices on the other hand, a variational model has to be proposed. The variational parameter can be determined analytically using the EFT, allowing to also study the vortex structure (variationally) throughout the BEC-BCS crossover.

1. Introduction

When cooling down a dilute cloud of fermionic atoms to ultralow temperatures, particles of different spin type can form Cooper pairs and condense into a superfluid state. The properties and features of these superfluid Fermi gases have been the subject of a considerable amount of theoretical and experimental research [1, 2]. The opportunity to investigate a whole continuum of inter-particle interaction regimes and the possibility to create a population imbalance result in an even richer physics than that of superfluid Bose gases. In this chapter, we present an effective field theory (EFT) suitable for the description of ultracold Fermi gases across the BEC-BCS interaction regime in a wide range of temperatures. The merits of this formalism mainly lie in the fact that it is computationally much less requiring than the Bogoliubov-de Gennes method, and that, in some cases, it can provide exact analytical solutions for the problem at hand. In Section 2, we give a short overview of the path integral theory that forms the basis for the EFT. In Section 3, we study the associated mean field theory for the description of homogeneous superfluids. In Section 4, we go beyond the mean-field approximation and describe the framework of the EFT. Sections 5 and 6 are dedicated to the application of the EFT to two important topological excitations: dark solitons and vortices.

2. Path integral theory and bosonification

The effective field theory for fermionic superfluids presented in this chapter is based on the path integral formalism of quantum field theory. The advantage of this formalism lies in the fact that the operators are replaced by fields, which can yield a more intuitive interpretation for the physics of the system. Moreover, the fact that there are no operators make working with functions of the quantum fields a lot easier.

In this section, the path integral description for ultracold Fermi gases will be briefly introduced. Using the Hubbard-Stratonovich identity, the fermionic degrees of freedom can be integrated out, resulting in an effective bosonic action. This effective bosonic action is the object of interest of this chapter and will lie at the basis of the effective field theory. An extended discussion of this section and the mean-field theory of the next section are given in an earlier publication [3]. Comprehensive introductions to the path integral method include [4] (Quantum Field Theory with Path Integrals), [5, 6] (The “classical” Path Integral), and [7] (General review book on the Path Integrals and most of its applications).

2.1. A brief introduction to the path integral formalism

The partition function of a system described by the quantum field action functional $S\left[\varphi \left(\mathbf{x},t\right),\overline{\varphi }\left(\mathbf{x},t\right)\right]$ can be expressed as a path integral [7]:

$$\mathcal{Z}=\int \mathcal{D}{\overline{\varphi }}_{\mathbf{x},\tau }\mathcal{D}{\varphi }_{\mathbf{x},\tau }\text{exp}\left(-S_E\left[{\varphi }_{\mathbf{x},\tau },{\overline{\varphi }}_{\mathbf{x},\tau }\right]\right).$$E1

Here, $\mathcal{D}{\varphi }_{\mathbf{x},\tau }$ represents a sum over all possible space-time configurations of the field $\varphi \left(\mathbf{x},\tau \right)$, and $\tau =\mathit{it}$ indicates imaginary times running from $\tau =0$ to $\tau =ħ\beta$ with $\beta =1/\left(k_{B}T\right)$. The Euclidian action $S_E\left(\beta \right)$ of the system is found from the real-time action functional $S\left(t_b,t_a\right)$ through the substitution

$$t\to -\mathit{i\tau }\Rightarrow S\left(t_b,t_a\right)\to {\mathit{\text{iS}}}_E\left(\beta \right).$$E2

For systems with an Euclidean action which is at most quadratic in the fields, the path integral (1) can be calculated analytically. In particular, two distinct cases can be considered:

Bosonic path integral: The path integral sums over a bosonic (scalar, complex valued) field $\mathrm{\Psi }\left(\mathbf{x},\tau \right)$:

$${\mathcal{Z}}_B=\int \mathcal{D}\overline{\mathrm{\Psi }}\mathcal{D}\mathrm{\Psi }\text{exp}\left(-\int \mathrm{d}\tau \int \mathrm{d}\mathbf{x}\int \mathrm{d}{\tau }^{\prime }\int \mathrm{d}{\mathbf{x}}^{\text{'}}\left[\overline{\mathrm{\Psi }}\left(\mathbf{x},\tau \right)\mathbb{A}(\mathbf{x},\tau ;{\mathbf{x}}^{\text{'}},{\tau }^{\prime }\left)\mathrm{\Psi }\right({\mathbf{x}}^{\text{'}},{\tau }^{\prime })\right]\right)=\frac{1}{\text{det}\left(\mathbb{A}\right)},$$E3

For the case of a quadratic bosonic path integral, the integration over the complex field $\mathrm{\Psi }$ reduces to a convolution of Gaussian integrals, which reduces to the inverse of the determinant of the matrix $\mathbb{A}$ containing the coefficients of the quadratic form.

Fermionic path integral: The path integral sums over a fermionic (Grassmann, complex valued) field $\psi \left(\mathbf{x},\tau \right)$:

$${\mathcal{Z}}_F=\int \mathcal{D}\overline{\psi }\int \mathcal{D}\psi \text{exp}\left(-\int \mathrm{d}\tau \int \mathrm{d}\mathbf{x}\int \mathrm{d}{\tau }^{\prime }\int \mathrm{d}{\mathbf{x}}^{\text{'}}\left[\overline{\psi }\left(\mathbf{x},\tau \right)\mathbb{A}(\mathbf{x},\tau ;{\mathbf{x}}^{\text{'}},{\tau }^{\prime }\left)\psi \right({\mathbf{x}}^{\text{'}},{\tau }^{\prime })\right]\right)=\text{det}\left(\mathbb{A}\right),$$E4

In the case of spin-dependent fermionic fields, the matrix $\mathbb{A}$ becomes slightly more complex since the spinor fields have multiple components¹ to account for the spin degree of freedom. The spinors $\psi$ are described by anti-commuting Grassmann numbers [4, 8], thus satisfying ${\psi }^2=0$. For the quadratic case, the fermionic path integral simply returns the determinant of the matrix $\mathbb{A}$.

Using the trace-log formula, these results can also be rewritten as:

$${\mathcal{Z}}_B=\text{exp}\left(-\mathrm{Tr}\left[\text{ln}\left(\mathbb{A}\right)\right]\right),$$E5

$${\mathcal{Z}}_F=\text{exp}\left(+\mathrm{Tr}\left[\text{ln}\left(\mathbb{A}\right)\right]\right).$$E6

Partition functions with quadratic action functionals form the basis of the path integral formalism. The usual approach for solving path integrals with higher order action functionals is to reduce them to the quadratic forms given above by the means of transformations and/or approximations.

In this chapter, the system of interest is an ultracold Fermi gas in which fermionic particles of opposite pseudo-spin interact via an $s$-wave contact potential. The Euclidian action functional for this system is given by

$$\begin{array}{ll}{S}_{\mathrm{E}}=& \underset{0}{\overset{\mathrm{\hslash }\beta }{\int }}\mathrm{d}\tau \int \mathrm{d}\mathbf{x}\sum _{\sigma \in \left\{\uparrow ,\downarrow \right\}}\left[{\overline{\psi }}_{\sigma }\left(\mathbf{x}\right)\left(\mathrm{\hslash }{\partial }_{\tau }-{\frac{{\mathrm{\hslash }}^2}{2m}}_{\mathbf{x}}{\nabla }_{\mathbf{x}}^2-{\mu }_{\sigma }\right)\right]{\psi }_{\sigma }\left(\mathbf{x}\right)\hfill \\ & +\underset{0}{\overset{\mathrm{\hslash }\beta }{\int }}\mathrm{d}\tau \int \mathrm{d}\mathbf{x}\int \mathrm{d}\mathbf{y}{\overline{\psi }}_{\uparrow }\left(\mathbf{x}\right){\overline{\psi }}_{\downarrow }\left(\mathbf{y}\right)\mathit{g\delta }\left(\mathbf{x}-\mathbf{y}\right){\psi }_{\downarrow }\left(\mathbf{y}\right){\psi }_{\uparrow }\left(\mathbf{x}\right),\hfill \end{array}$$E7

where $\sigma \in \left\{\uparrow ,\downarrow \right\}$ denotes the spin components of the fermionic spinor fields, the chemical potentials ${\mu }_{\sigma }$ fix the amount of particles of each spin population, and $g$ is the renormalized interaction strength [9, 10], linking the interaction potential to the s-wave scattering length $a_s$:

$$\frac{1}{g}=\frac{m}{4\pi {\mathrm{\hslash }}^2{a}_s}-\int \frac{\mathrm{d}\mathbf{k}}{{\left(2\pi \right)}^3}\frac{m}{{\mathrm{\hslash }}^2{k}^2}.$$E8

For the remainder of the chapter, the units

$$\mathrm{\hslash }=k_B=k_F=2m=1$$E9

will be used, meaning that we work in the natural units of $k_F$, $E_F$, ${\omega }_F=E_F/\mathrm{\hslash }$, and $T_F=E_F/k_B$. Consequentially, the partition function of the ultracold Fermi gas can be written down as

$$\begin{array}{ll}\mathcal{Z}=& \int \mathcal{D}{\overline{\psi }}_{\sigma }\mathcal{D}{\psi }_{\sigma }\text{exp}\left[-\underset{0}{\overset{\beta }{\int }}\mathrm{d}\tau \int \mathrm{d}\mathbf{x}\left(\sum _{\sigma \in \left\{\uparrow ,\downarrow \right\}}{\overline{\psi }}_{\sigma }\left(\mathbf{x},\tau \right)\left(\frac{\partial }{\partial \tau }-{\nabla }_{\mathrm{x}}^2-{\mu }_{\sigma }\right){\psi }_{\sigma }(\mathbf{x},\tau )\right.\right.\hfill \\ & \left.\left.+g{\overline{\psi }}_{\uparrow }\left(\mathbf{x},\tau \right){\overline{\psi }}_{\downarrow }(\mathbf{x},\tau \left){\psi }_{\downarrow }\right(\mathbf{x},\tau \left){\psi }_{\uparrow }\right(\mathbf{x},\tau )\right)\right],\end{array}$$E10

where the label $\sigma$ was explicitly added to the integration measure to show that the path integration is performed also over both spin components of the spinor $\psi$. As noted above, only quadratic path integrals can be solved analytically, meaning that an additional trick is needed² to calculate the above partition sum (10). In the present treatment, this trick will be the Hubbard-Stratonovich transformation.

2.2. Bosonification: the Hubbard-Stratonovich transformation

Using the Hubbard-Stratonovich identity [11, 12, 13, 14],

$$\text{exp}\left(-g\int {\mathrm{d}}^3\mathit{xg}{\overline{\psi }}_{\uparrow }{\overline{\psi }}_{\downarrow }{\psi }_{\downarrow }{\psi }_{\uparrow }\right)=\int \mathcal{D}\overline{\mathrm{\Psi }}\mathcal{D}\mathrm{\Psi }\text{exp}\left(\int {\mathrm{d}}^3\mathbf{x}\left[\frac{{\left|\mathrm{\Psi }\right|}^2}{g}+{\overline{\psi }}_{\uparrow }{\overline{\psi }}_{\downarrow }\mathrm{\Psi }+{\psi }_{\downarrow }{\psi }_{\uparrow }\overline{\mathrm{\Psi }}\right]\right),$$E11

it is possible to rewrite the action in a form that is quadratic in the fermionic fields $\psi$ and $\overline{\psi }$, allowing for the fermionic degrees of freedom to be integrated out. The price of this transformation is the introduction of a new (auxiliary) bosonic field $\mathrm{\Psi }\left(\mathbf{r},\tau \right)$, which can be interpreted as the field of the Cooper pairs that will form the superfluid state. Diagramatically, the Hubbard-Stratonovich identity removes the four-point vertex (quartic interaction term) and replaces it with two three-point vertices (quadratic terms), as illustrated in Figure 1. It is important to note that, although the Hubbard-Stratonovich transformation is an exact identity, further calculations will require approximations for which the choice of collective field (or “channel”) becomes important. Whereas the bosonic pair field is suitable for the superfluid state, it will fail when one tries to use it to take into account interactions in the normal state. It should therefore be pointed out that alternatives exist, notably Kleinert’s variational perturbation theory, in which a classical collective field rather than a quantum collective field is used. This allows for the simultaneous treatment of multiple collective fields [15], for example, the pair field and the density field. For our present purposes, however, it is sufficient to restrict ourselves to the superfluid state and describe it with a single collective field.

Figure 1.

A diagrammatic representation of the different terms in the Hubbard-Stratonovich identity (11).

After applying the Hubbard-Stratonovich identity (11) to expression (10), the partition function becomes

$$\begin{array}{ll}\mathcal{Z}=& \int \mathcal{D}{\overline{\psi }}_{\sigma }\mathcal{D}{\psi }_{\sigma }\int \mathcal{D}\overline{\mathrm{\Psi }}\mathcal{D}\mathrm{\Psi }\text{exp}\left[-\underset{0}{\overset{\beta }{\int }}\mathrm{d}\tau \int \mathrm{d}\mathbf{x}\left(\sum _{\sigma \in \left\{\uparrow ,\downarrow \right\}}{\overline{\psi }}_{\sigma }\left(\mathbf{x},\tau \right)\left(\frac{\partial }{\partial \tau }-{\nabla }_{\mathbf{x}}^2-{\mu }_{\sigma }\right){\psi }_{\sigma }(\mathbf{x},\tau )\right.\right.\hfill \\ & \left.\left.-\frac{{\left|\mathrm{\Psi }\left(\mathbf{x},\tau \right)\right|}^2}{g}-{\overline{\psi }}_{\uparrow }\left(\mathbf{x},\tau \right){\overline{\psi }}_{\downarrow }(\mathbf{x},\tau \left)\mathrm{\Psi }\right(\mathbf{x},\tau )-{\psi }_{\downarrow }(\mathbf{x},\tau \left){\psi }_{\uparrow }\right(\mathbf{x},\tau \left)\overline{\mathrm{\Psi }}\right(\mathbf{x},\tau )\right)\right]\hfill \end{array}$$E12

2.3. The resulting bosonic path integral

Since the path integral over the fermionic fields $\psi$ and $\overline{\psi }$ is now quadratic, it can be performed analytically using formula (4), resulting in the effective bosonic path integral [3]

$$\mathcal{Z}=\int \mathcal{D}\overline{\mathrm{\Psi }}\mathcal{D}\mathrm{\Psi }\text{exp}\left[-\left(-\underset{0}{\overset{\beta }{\int }}\int \mathrm{d}\mathbf{x}\frac{{\left|\mathrm{\Psi }\left(\mathbf{x},\tau \right)\right|}^2}{g}-\mathrm{Tr}\left[\text{ln}\left(-{\mathbb{G}}^{-1}\right)\right]\right)\right],$$E13

where the components of the inverse Green’s function matrix $-{\mathbb{G}}^{-1}$ are given by

$$-{\mathbb{G}}^{-1}\left(\mathbf{x},\tau \right)=\left(\begin{array}{ll}\frac{\partial }{\partial \tau }-{\mathbf{\nabla }}_{\mathbf{x}}^2-{\mu }_{\uparrow }& -\mathrm{\Psi }\left(\mathbf{x},\tau \right)\\ -\overline{\mathrm{\Psi }}\left(\mathbf{x},\tau \right)& \frac{\partial }{\partial \tau }+{\mathbf{\nabla }}_{\mathbf{x}}^2+{\mu }_{\downarrow }\end{array}\right)$$E14

Since $-{\mathbb{G}}^{-1}$ depends on the bosonic field $\mathrm{\Psi }\left(\mathbf{x},\tau \right)$, the action in the exponent is not quadratic, and hence, the remaining bosonic path integral can still not be solved analytically. In order to obtain a workable solution, two different approximations will be considered. First, a mean field approximation (using a constant value for $\mathrm{\Psi }$) will be discussed in Section 3. Subsequently, this mean field theory will form the basis for a finite temperature effective field theory, which also takes into account slow fluctuations of the pair field $\mathrm{\Psi }\left(\mathbf{x},\tau \right)$. This theory will be presented in Section 4.

3. The mean field theory

At first sight, the introduction of the auxiliary bosonic fields $\mathrm{\Psi }\left(\mathbf{x},\tau \right)$ and $\overline{\mathrm{\Psi }}\left(\mathbf{x},\tau \right)$ through the Hubbard-Stratonovich transformation seems to have been of little use; while the transformation enables us to perform the path integrals over the fermionic fields, we end up with path integrals for $\mathrm{\Psi }\left(\mathbf{x},\tau \right)$ and $\overline{\mathrm{\Psi }}\left(\mathbf{x},\tau \right)$ that still cannot be calculated exactly. The advantage of switching to the bosonic pair fields, however, lies in the fact that they allow us to make a physically plausible approximation based on our knowledge of the system. If we want to investigate the superfluid state, we can assume that the most important contribution to the path integral will come from the configuration in which all the bosonic pairs are condensed into the lowest energy state of the system and form a homogeneous superfluid. This assumption is most easily expressed in momentum-frequency representation $\left\{\mathbf{q},m\right\}$:

$$\mathrm{\Psi }\left(\mathbf{q},m\right)\to \sqrt{\mathit{\beta V}}\delta \left(\mathbf{q}\right){\delta }_{m,0}\times \mathrm{\Delta },$$E15

$$\overline{\mathrm{\Psi }}\left(\mathbf{q},m\right)\to \sqrt{\mathit{\beta V}}\delta \left(\mathbf{q}\right){\delta }_{m,0}\times {\mathrm{\Delta }}^{\ast },$$E16

where $m$ characterizes the bosonic Matsubara frequencies ${\tilde{\omega }}_m=2\mathit{m\pi }/\beta$, and $V$ represents the volume of the system. This approximation, which is called the saddle-point approximation for the bosonic path integral, comes down to assuming that the pair field $\mathrm{\Psi }\left(\mathbf{x},\tau \right)$ takes on a constant value $\mathrm{\Delta }$. By applying this approximation to the bosonic path integral in expression (12) (i.e., after performing the Hubbard-Stratonovich transformation but before performing the Grassmann integration over the fermionic fields), the resulting fermionic path integral can be solved analytically using formula (4) to find the saddle-point expression for the partition function:

$${\mathcal{Z}}_{\mathrm{sp}}=\text{exp}\left\{\frac{{\left|\mathrm{\Delta }\right|}^2}{g}-\sum _{\mathbf{k},n}\text{ln}\left[\left(i{\omega }_n-E_{\mathbf{k}}+\zeta \right)\left(-i{\omega }_n-E_{\mathbf{k}}-\zeta \right)\right]\right\}.$$E17

where ${\omega }_n$ are the fermionic Matsubara frequencies ${\omega }_n=\left(2n+1\right)\pi /\beta$. We have also introduced the single-particle excitation energy $E_{\mathbf{k}}=\sqrt{{\xi }_{\mathbf{k}}^2+{\mathrm{\Delta }}^2}$ with ${\xi }_{\mathbf{k}}=k^2-\mu$, and we have defined the average chemical potential $\mu$ and the imbalance chemical potential $\zeta$ as

$$\mu =\frac{{\mu }_{\uparrow }+{\mu }_{\downarrow }}{2}\text{and}\zeta =\frac{{\mu }_{\uparrow }-{\mu }_{\downarrow }}{2}.$$E18

The parameter $\zeta$ determines the population imbalance between the two spin populations. For $\zeta =0$, the numbers of particles of each spin type are equal, while for non-zero values of $\zeta$, there will be more spin-up than spin-down particles or vice versa. The saddle-point partition function can now be rewritten in terms of the saddle-point thermodynamic potential per unit volume ${\mathrm{\Omega }}_{\mathrm{sp}}$ as

$${\mathcal{Z}}_{\mathrm{sp}}=\text{exp}\left\{-\mathit{\beta V}{\mathrm{\Omega }}_{\mathrm{sp}}\right\}.$$E19

After performing the Matsubara summation over $n$ [3] and replacing the sum over $\mathbf{k}$ by a continuous integral in expression (17), we finally find for ${\mathrm{\Omega }}_{\mathrm{sp}}$:

$${\mathrm{\Omega }}_{\mathrm{sp}}=-\frac{{\left|\mathrm{\Delta }\right|}^2}{8\pi k_F{a}_s}-\int \frac{\mathrm{d}\mathbf{k}}{{\left(2\pi \right)}^3}\left\{\frac{1}{\beta }\left[2\text{cosh}\left(\beta E_{\mathbf{k}}\right)+2\text{cosh}\left(\mathit{\beta \zeta }\right)\right]-{\xi }_{\mathbf{k}}-\frac{{\left|\mathrm{\Delta }\right|}^2}{2k^2}\right\}$$E20

The saddle-point value ${\mathrm{\Delta }}_{\mathrm{sp}}$ for the pair field is found through the requirement that ${\mathrm{\Delta }}_{\mathrm{sp}}$ minimizes ${\mathrm{\Omega }}_{\mathrm{sp}}$, which yields the gap equation:

$${\left.\frac{\partial {\mathrm{\Omega }}_{\mathrm{sp}}}{\mathrm{\partial \Delta }}\right|}_{T,\mu ,\zeta }=0$$E21

This is illustrated in Figure 2, which shows the thermodynamic potential ${\mathrm{\Omega }}_{\mathrm{sp}}$ as a function of $\mathrm{\Delta }$ for several values of the imbalance chemical potential $\zeta$. The superfluid state exists when ${\mathrm{\Omega }}_{\mathrm{sp}}$ reaches its minimum at a nonzero value of $\mathrm{\Delta }$. As $\zeta$ is increased, the normal state at $\mathrm{\Delta }=0$ develops and becomes the global minimum above a critical imbalance level. This transition from the superfluid to the normal state under influence of increasing population imbalance is known as the Clogston phase transition [16].

Figure 2.

The thermodynamic potential ${\mathrm{\Omega }}_{\mathrm{sp}}$ in function of $\mathrm{\Delta }$ for several values of the imbalance chemical potential $\zeta$, at temperature $T/T_F=0.01$ and chemical potential $\mu =1.3E_F$. The evolution of the normal state at $\mathrm{\Delta }=0$ as $\zeta$ increases illustrates the Clogston phase transition.

When working with a fixed number of particles, the chemical potential $\mu$ and the imbalance chemical potential $\zeta$ have to be related to the fermion density $n_{\mathrm{sp}}$ and density difference $\delta n_{\mathrm{sp}}$ (between the two spin populations) through the number equations

$$n_{\mathrm{sp}}=-{\left.\frac{\partial {\mathrm{\Omega }}_{\mathrm{sp}}}{\partial \mu }\right|}_{T,\zeta ,\mathrm{\Delta }}$$E22

$$\delta n_{\mathrm{sp}}=-{\left.\frac{\partial {\mathrm{\Omega }}_{\mathrm{sp}}}{\partial \zeta }\right|}_{T,\mu ,\mathrm{\Delta }}$$E23

Since in our units $k_F=1$, the particle density $n_{\mathrm{sp}}$ is fixed by $n_{\mathrm{sp}}=1/\left(3{\pi }^2\right)$. Given the input parameters $\beta$, $\zeta$, and $a_s$, the values $\mathrm{\Delta }$ and $\mu$ can then be found from the coupled set of Eqs. (21) and (22), while (23) fixes $\delta n_{\mathrm{sp}}$ as a function of $\zeta$. Solutions for ${\mathrm{\Delta }}_{\mathrm{sp}}$ and $\mu$ across the BEC-BCS crossover are shown in Figure 3a and b.

Figure 3.

Solutions for the pair field $\mathrm{\Delta }$ and the average chemical potential $\mu$ in function of the interaction strength ${\left(k_F{a}_s\right)}^{-1}$ at temperature $T/T_F=0.01$. The solution for $\mathrm{\Delta }$ is shown for several values of the imbalance chemical potential $\zeta$, illustrating the transition from the superfluid to the normal state under influence of population imbalance.

4. The effective field theory

While the saddle-point approximation is a suitable model for the qualitative description of homogeneous Fermi superfluids, it does not account for the effects of fluctuations of the order parameter, nor does it include any excitations other than the single-particle Bogoliubov excitations. To study the properties and dynamics of non-homogeneous systems, one needs to go beyond the limitations of a mean field theory. In this section, we formulate an effective field theory (EFT) for the pair field $\mathrm{\Psi }\left(\mathbf{r},t\right)$ that can describe nonhomogeneous Fermi superfluids in the BEC-BCS crossover at finite temperatures. To this end, we return to the path integral expression (13) for the partition function, which was obtained after performing the Hubbard-Stratonovich transformation and integrating out the fermionic degrees of freedom. Since the exponent of this partition function only depends on the fields $\mathrm{\Psi }\left(\mathbf{r},t\right)$ and $\overline{\mathrm{\Psi }}\left(\mathbf{r},t\right)$, we can define an effective bosonic action for the pair field given by

$$S_{\mathrm{eff}}=S_B-\mathrm{Tr}\left[\text{ln}\left(-{\mathbb{G}}^{-1}\right)\right],$$E24

where $S_B=-\underset{0}{\overset{\beta }{\int }}\int \mathrm{d}\mathbf{x}\frac{{\left|\mathrm{\Psi }\left(\mathbf{x},\tau \right)\right|}^2}{8}$ is the action for free bosonic fields. The inverse Green’s function matrix $-{\mathbb{G}}^{-1}$ for interacting fermions, which was defined in expression (14), can be separated into its diagonal and off-diagonal components

$$\begin{array}{c}-{\mathbb{G}}^{-1}(x,\tau )=-{\mathbb{G}}_0^{-1}(x,\tau )+\mathbb{F}(x,\tau )\\ =\left(\begin{array}{cc}\frac{\partial }{\partial \tau }-{\nabla }_x^2-{\mu }_{\uparrow }& 0\\ 0& \frac{\partial }{\partial \tau }+{\nabla }_x^2+{\mu }_{\downarrow }\end{array}\right)+\left(\begin{array}{cc}0& -\Psi (x,\tau )\\ -\overline{\Psi }(x,\tau )& 0\end{array}\right),\end{array}$$E25

where $-{\mathbb{G}}_0^{-1}$ describes free fermionic fields, while $\mathbb{F}$ describes the pairing of the fermions. Using this decomposition, we can write the effective bosonic action functional (24) as

$$\begin{array}{c}{S}_{eff}=S_B-Tr\text{[}\ln \left(-{\mathbb{G}}_0^{-1}+\mathbb{F}\right)\text{]}\\ =S_B-Tr\text{[}\ln \left(-{\mathbb{G}}_0^{-1}\right)\text{]}-Tr\text{[}\ln \left(1-{\mathbb{G}}_0\mathbb{F}\right)\text{]}\\ =S_B+S_0+{\displaystyle \sum }_{p=1}^{\infty }\frac{1}{p}Tr\text{[}{\left({\mathbb{G}}_0\mathbb{F}\right)}^p\text{]}.\end{array}$$E26

While, in general, this infinite sum over all powers of the pair field cannot be calculated analytically, there exist many possible approximations that lead to various theoretical treatments of the ultracold Fermi gas. For example, the mean field saddle-point approximation from the previous section can be retrieved by simply setting

$$\mathbb{F}\left(\mathbf{x},\tau \right)\approx {\mathbb{F}}_{\mathrm{sp}}=\left(\begin{array}{ll}0& -\mathrm{\Delta }\\ -\overline{\mathrm{\Delta }}& 0\end{array}\right)$$E27

in (26) and calculating the whole sum over $p$. In the Ginzburg-Landau (GL) treatment for ultracold Fermi gases, the action is approximated by assuming small fluctuations of the pair field $\mathrm{\Psi }\left(\mathbf{x},\tau \right)$ around the normal state $\mathrm{\Psi }=0$. This assumption comes down to keeping only terms up to $p=2$ in the sum in (26) and approximating $\mathbb{F}\left(\mathbf{x},\tau \right)$ by the following gradient expansion

$$\begin{array}{c}{\left.\mathbb{F}\left(\mathbf{x},\tau \right)\approx {\mathbb{F}}_0+\left(\mathbf{x}-{\mathbf{x}}_0\right)\cdot {\mathbf{\nabla }}_{\mathbf{x}}\mathbb{F}\right|}_{{\mathbf{x}}_{\mathbf{0}}}+\frac{1}{2}\sum _{i,j=x,y,z}\left(x_i-x_{0,i}\right)\left(x_j-x_{0,j}\right){\left.\frac{{\partial }^2\mathbb{F}}{\partial x_i\partial x_j}\right|}_{{\mathbf{x}}_0}\\ +\left(\tau -{\tau }_0\right){\left.\frac{\partial \mathbb{F}}{\partial \tau }\right|}_{{\tau }_0}+\frac{1}{2}{\left(\tau -{\tau }_0\right)}^2{\left.\frac{{\partial }^2\mathbb{F}}{\partial {\tau }^2}\right|}_{{\tau }_0},\end{array}$$E28

with ${\mathbb{F}}_0\to 0$. The result is an effective field treatment which is valid close to the critical temperature $T_c$ of the superfluid phase transition. Inspired by the GL formalism, we will now present a beyond saddle-point EFT that is capable of describing Fermi superfluids in the BEC-BCS crossover at finite temperatures. This theory is based on the assumption that the pair field $\mathrm{\Psi }\left(\mathbf{x},\tau \right)$ exhibits slow variations in space and time around a constant bulk value. Since this is a weaker condition than the GL assumption of small variations, it is ultimately expected to lead to a larger applicability domain. The assumption of slow fluctuations is implemented through a gradient expansion of the pair field around its saddle-point value, similar to (28) but with ${\mathbb{F}}_0\to {\mathbb{F}}_{\mathrm{sp}}$. Subsequently, we consider the full infinite sum in (26):

$$\sum _{p=1}^{\infty }\frac{1}{p}\mathrm{Tr}\left[{\left({\mathbb{G}}_0\mathbb{F}\right)}^p\right]=\sum _{p=1}^{\infty }\frac{1}{p}\mathrm{Tr}\left[\underset{p\text{factors}}{\underbrace{{\mathbb{G}}_0\mathbb{F}{\mathbb{G}}_0\mathbb{F}\dots {\mathbb{G}}_0\mathbb{F}}}\right].$$E29

In every term of this sum, we replace (at most) two occurrences of $\mathbb{F}\left(\mathbf{x},\tau \right)$ by its gradient expansion and substitute all remaining factors $\mathbb{F}\left(\mathbf{x},\tau \right)$ by ${\mathbb{F}}_{\mathrm{sp}}$. Afterward, the entire sum over $p$ can be calculated analytically. The result of this calculation, the details of which can be found in [17], is an explicit expression for the Euclidian action functional that governs the dynamics of the pair field $\mathrm{\Psi }\left(\mathbf{x},\tau \right)$ of a three-dimensional (3D) superfluid Fermi gas:

$$\begin{array}{ll}{S}_{\mathrm{EFT}}=& {\int }_0^{\beta }\mathrm{d}\tau \int \mathrm{d}\mathbf{x}\left[\frac{D}{2}\left(\overline{\mathrm{\Psi }}\frac{\mathrm{\partial \Psi }}{\partial \tau }-\frac{\partial \overline{\mathrm{\Psi }}}{\partial \tau }\mathrm{\Psi }\right)+{\mathrm{\Omega }}_s+\frac{C}{2m}\left({\mathbf{\nabla }}_{\mathbf{x}}\overline{\mathrm{\Psi }}\cdot {\mathbf{\nabla }}_{\mathbf{x}}\mathrm{\Psi }\right)-\frac{E}{2m}{\left({\mathbf{\nabla }}_{\mathbf{x}}{\left|\mathrm{\Psi }\right|}^2\right)}^2\right.\hfill \\ & +\left.Q\frac{\partial \overline{\mathrm{\Psi }}}{\partial \tau }\frac{\mathrm{\partial \Psi }}{\partial \tau }-R{\left(\frac{\partial {\left|\mathrm{\Psi }\right|}^2}{\partial \tau }\right)}^2\right].\hfill \end{array}$$E30

The EFT coefficients ${\mathrm{\Omega }}_s$, $C$, $D$, $E$, $Q$ and $R$ are given by

$${\mathrm{\Omega }}_s=-\frac{1}{8\pi k_F{a}_s}{\left|\mathrm{\Delta }\right|}^2-\int \frac{d\mathbf{k}}{{\left(2\pi \right)}^3}\left\{\frac{1}{\beta }\text{ln}\left[2\text{cosh}\left(\beta E_{\mathbf{k}}\right)+2\text{cosh}\left(\mathit{\beta \zeta }\right)\right]-{\xi }_{\mathbf{k}}-\frac{{\left|\mathrm{\Delta }\right|}^2}{2k^2}\right\}$$E31

$$C=\int \frac{d\mathbf{k}}{{\left(2\pi \right)}^3}\frac{k^2}{3m}{f}_2\left(\beta ,E_{\mathbf{k}},\zeta \right)$$E32

$$D=\int \frac{d\mathbf{k}}{{\left(2\pi \right)}^3}\frac{{\xi }_{\mathbf{k}}}{{\left|\mathrm{\Psi }\right|}^2}\left[f_1\left(\beta ,{\xi }_{\mathbf{k}},\zeta \right)-f_1\left(\beta ,E_{\mathbf{k}},\zeta \right)\right]$$E33

$$E=2\int \frac{d\mathbf{k}}{{\left(2\pi \right)}^3}\frac{k^2}{3m}{\xi }_{\mathbf{k}}^2{f}_4\left(\beta ,E_{\mathbf{k}},\zeta \right)$$E34

$$Q=\frac{1}{2{\left|\mathrm{\Psi }\right|}^2}\int \frac{d\mathbf{k}}{{\left(2\pi \right)}^3}\left[f_1\left(\beta ,E_{\mathbf{k}},\zeta \right)-\left(E_{\mathbf{k}}^2+{\xi }_{\mathbf{k}}^2\right)f_2\left(\beta ,E_{\mathbf{k}},\zeta \right)\right]$$E35

$$R=\frac{1}{2{\left|\mathrm{\Psi }\right|}^2}\int \frac{d\mathbf{k}}{{\left(2\pi \right)}^3}\left[\frac{f_1\left(\beta ,E_{\mathbf{k}},\zeta \right)+\left(E_{\mathbf{k}}^2-3{\xi }_{\mathbf{k}}^2\right)f_2\left(\beta ,E_{\mathbf{k}},\zeta \right)}{3{\left|\mathrm{\Psi }\right|}^2}\right.$$

$$+\left.\frac{4\left({\xi }_{\mathbf{k}}^2-2E_{\mathbf{k}}^2\right)}{3}{f}_3\left(\beta ,E_{\mathbf{k}},\zeta \right)+2E_{\mathbf{k}}^2{\left|\mathrm{\Psi }\right|}^2{f}_4(\beta ,E_{\mathbf{k}},\zeta )\right],$$E36

where the functions $f_p\left(\beta ,\epsilon ,\zeta \right)$ are recursively defined as

$$f_1\left(\beta ,\mathrm{ϵ},\zeta \right)=\frac{1}{2\mathrm{ϵ}}\frac{\text{sinh}\left(\beta \mathrm{ϵ}\right)}{\text{cosh}\left(\beta \mathrm{ϵ}\right)+\text{cosh}\left(\mathit{\beta \zeta }\right)}$$E37

$$f_{p+1}\left(\beta ,\mathrm{ϵ},\zeta \right)=-\frac{1}{2p\mathrm{ϵ}}\frac{\partial f_p\left(\beta ,\mathrm{ϵ},\zeta \right)}{\mathrm{\partial ϵ}}$$E38

In general, each of these EFT coefficients depends on the modulus squared of the order parameter ${\left|\mathrm{\Psi }\left(\mathbf{x},\tau \right)\right|}^2$. In practice, however, we will assume that the coefficients associated with the second order derivatives of the pair field can be kept constant and equal to their bulk value, since retaining their full space-time dependence would strictly speaking lead us beyond the second-order approximation of the gradient expansion. This means that in expressions (32), (34), (35), and (36) for the coefficients $C$, $E$, $Q$, and $R$, we set ${\left|\mathrm{\Psi }\left(\mathbf{x},\tau \right)\right|}^2={\left|{\mathrm{\Psi }}_{\infty }\right|}^2$ and $E_{\mathbf{k}}=\sqrt{{\xi }_{\mathbf{k}}^2+{\left|{\mathrm{\Psi }}_{\infty }\right|}^2}$, where ${\left|{\mathrm{\Psi }}_{\infty }\right|}^2$ is the saddle-point value of the pair field for a uniform system. For the coefficients ${\mathrm{\Omega }}_s$ and $D$ on the other hand, the full space-time dependence of ${\left|\mathrm{\Psi }\left(\mathbf{x},\tau \right)\right|}^2$ is preserved.

The effective action functional (30) forms the basis of our EFT description of superfluid Fermi gases. The validity and limitations of the formalism are largely determined by the main assumption that the order parameter varies slowly in space and time, which corresponds to the condition that the pair field should vary over a spatial region much larger than the Cooper pair correlation length. A detailed study of the limitations imposed by this condition was carried out in [18]. In the following chapters, we will demonstrate some of the ways in which the EFT can be employed by applying it to the description of two important topological excitations of the superfluid: dark solitons and vortices.

5. Application 1: Soliton dynamics

In this section, we will use the EFT that was developed in Section 4 to study the properties of dark solitons in Fermi superfluids.

5.1. What is a dark soliton?

Solitons are nonlinear solitary waves that maintain their shape while propagating through a medium at a constant velocity. They are found as the solution of nonlinear wave equations and emerge in a wide variety of physical systems, including optical fibers, classical fluids, and plasmas. More recently, they have also become a subject of interest in superfluid quantum gases [19, 20, 21, 22, 23]. In these systems, solitons appear most often in the form of dark solitons, which are characterized by a localized density dip in the uniform background and a jump in the phase profile of the order parameter. The magnitude of this density dip and phase jump are intrinsically connected to the velocity $v_s$ with which the soliton propagates through the superfluid, as illustrated in Figure 4. The higher the soliton velocity, the smaller the phase jump and soliton depth become. Above a certain critical velocity $v_c$, the phase jump and density dip will disappear completely and a dark soliton solution no longer exists.

Figure 4.

Example of the density profile (upper row) and phase profile (lower row) of a dark soliton for different soliton velocities $v_s$ relative to the critical velocity $v_c$.

5.2. Solution for a one-dimensional dark soliton

For the case of a dark soliton in a one-dimensional (1D) Fermi superfluid with a uniform background, the EFT provides an exact analytical solution for the pair field [24]. To describe the dynamics of the system, it is necessary to move from the imaginary-time action functional (30) to the real-time one, using the formal replacements.

$$\tau \to \mathit{it}$$E39

$$S_{\mathrm{EFT}}\left(\beta \right)\to -i{S}_{\mathrm{EFT}}\left(t_b,t_a\right).$$E40

From the relation between the real-time action functional and the Lagrangian density $\mathcal{L}$,

$$S_{\mathrm{EFT}}\left(t_b,t_a\right)={\int }_{t_a}^{t_b}\mathrm{d}t\int \mathrm{d}\mathbf{x}\mathcal{L},$$E41

we subsequently find the following expression for $\mathcal{L}$:

$$\begin{array}{c}\mathcal{L}=i\frac{D}{2}\left(\overline{\mathrm{\Psi }}\frac{\partial \mathrm{\Psi }}{\partial t}-\frac{\partial \overline{\mathrm{\Psi }}}{\partial t}\mathrm{\Psi }\right)-{\mathrm{\Omega }}_s-\frac{C}{2m}\left({\mathbf{\nabla }}_{\mathbf{x}}\overline{\mathrm{\Psi }}\cdot {\mathbf{\nabla }}_{\mathbf{x}}\mathrm{\Psi }\right)+\frac{E}{2m}{\left({\mathbf{\nabla }}_{\mathbf{x}}{\left|\mathrm{\Psi }\right|}^2\right)}^2\\ +Q\frac{\partial \overline{\mathrm{\Psi }}}{\partial t}\frac{\partial \mathrm{\Psi }}{\partial t}-R{\left(\frac{\partial {\left|\mathrm{\Psi }\right|}^2}{\partial t}\right)}^2,\end{array}$$E42

where the Hamiltonian density $\mathcal{H}$ is defined as

$$\mathcal{H}={\mathrm{\Omega }}_s+\frac{C}{2m}\left({\mathbf{\nabla }}_{\mathbf{x}}\overline{\mathrm{\Psi }}\cdot {\mathbf{\nabla }}_{\mathbf{x}}\mathrm{\Psi }\right)-\frac{E}{2m}{\left({\mathbf{\nabla }}_{\mathbf{x}}{\left|\mathrm{\Psi }\right|}^2\right)}^2+Q\frac{\partial \overline{\mathrm{\Psi }}}{\partial t}\frac{\partial \mathrm{\Psi }}{\partial t}-R{\left(\frac{\partial {\left|\mathrm{\Psi }\right|}^2}{\partial t}\right)}^2.$$E43

As mentioned above, a dark soliton in a superfluid is mainly characterized by a jump in the phase profile and a dip in the amplitude profile of the order parameter. Therefore, it is convenient to write the pair field $\mathrm{\Psi }\left(\mathbf{x},t\right)$ as

$$\mathrm{\Psi }\left(\mathbf{x},t\right)=\mid \mathrm{\Psi }\left(\mathbf{x},t\right)\mid e^{\mathit{i\theta }\left(\mathbf{x},t\right)}.$$E44

Moreover, since a soliton is a localized perturbation, we write the modulus as a product of the constant background value $\mid {\mathrm{\Psi }}_{\infty }\mid$ and a relative amplitude $a\left(\mathbf{x},t\right)$ that modifies the background value at the position of the soliton:

$$\mid \mathrm{\Psi }\left(\mathbf{x},t\right)\mid =\mid {\mathrm{\Psi }}_{\infty }\mid a\left(\mathbf{x},t\right).$$E45

Substituting this form for the pair field in the field Lagrangian (42), we find

$$\begin{array}{ll}\mathcal{L}=& -\kappa \left(a\right)a^2\frac{\partial \theta }{\partial t}-\left[{\mathrm{\Omega }}_s\left(a\right)-{\mathrm{\Omega }}_s\left(a_{\infty }\right)\right]-\frac{1}{2}{\rho }_{\mathrm{qp}}\left(a\right){\left({\mathbf{\nabla }}_{\mathbf{x}}a\right)}^2-\frac{1}{2}{\rho }_{\mathrm{sf}}\left(a\right){\left({\mathbf{\nabla }}_{\mathbf{x}}\theta \right)}^2\hfill \\ & +\left(Q-4R{\left|{\mathrm{\Psi }}_{\infty }\right|}^2{a}^2\right){\left|{\mathrm{\Psi }}_{\infty }\right|}^2{\left(\frac{\partial a}{\partial t}\right)}^2+Q{\left|{\mathrm{\Psi }}_{\infty }\right|}^2{a}^2{\left(\frac{\partial \theta }{\partial t}\right)}^2,\hfill \end{array}$$E46

with

$$\kappa \left(a\right)=D\left(a\right){\left|{\mathrm{\Psi }}_{\infty }\right|}^2,$$E47

$${\rho }_{\mathrm{qp}}\left(a\right)=\frac{C-4E{\left|{\mathrm{\Psi }}_{\infty }\right|}^2{a}^2}{m}{\left|{\mathrm{\Psi }}_{\infty }\right|}^2,$$E48

$${\rho }_{\mathrm{sf}}\left(a\right)=\frac{C}{m}{\left|{\mathrm{\Psi }}_{\infty }\right|}^2{a}^2.$$E49

Here, we added ${\mathrm{\Omega }}_s\left(a_{\infty }\right)$ to the original Lagrangian to obtain a regularized Lagrangian density in which energy values are considered with respect to the energy of the uniform system. The superfluid density ${\rho }_{\mathrm{sf}}$ determines how much the pair condensate resists gradients in its phase field, while the quantum pressure ${\rho }_{\mathrm{qp}}$ is a consequence of the fact that the condensate also resists gradients in the pair density. We will further limit ourselves to a 1D problem in which the soliton propagates with constant speed $v_s$ in the $x$-direction on a uniform background. This assumption can be implemented through the condition that the space-time dependence of the pair field satisfies the relation $f\left(x,t\right)=f\left(x-v_{s}t\right)$. We then perform a change of variables $x^{\prime }=x-v_{s}t$ and $t^{\prime }=t$, corresponding to a transformation to the frame of reference that moves along with the soliton and has its origin at the soliton center. It follows that

$$f\left(x-v_{s}t\right)=f\left(x^{\prime }\right),\frac{\partial }{\partial x}=\frac{\partial }{\partial x^{\prime }},\frac{\partial }{\partial t}=\frac{\partial }{\partial t^{\prime }}-v_s\frac{\partial }{\partial x^{\prime }}.$$E50

If we further drop the primes, the Lagrangian density (46) in the soliton frame of reference can be written as

$$\mathcal{L}=\kappa \left(a\right)a^2{v}_s\frac{\partial \theta }{\partial x}-\left[{\mathrm{\Omega }}_s\left(a\right)-{\mathrm{\Omega }}_s\left(a_{\infty }\right)\right]-\frac{1}{2}{\tilde{\rho }}_{\mathrm{qp}}\left(a\right){\left(\frac{\partial a}{\partial x}\right)}^2-\frac{1}{2}{\tilde{\rho }}_{\mathrm{sf}}\left(a\right){\left(\frac{\partial \theta }{\partial x}\right)}^2.$$E51

with the modified superfluid density and quantum pressure

$${\tilde{\rho }}_{\mathrm{qp}}\left(a\right)=\frac{C-4E{\left|{\mathrm{\Psi }}_{\infty }\right|}^2{a}^2}{m}{\left|{\mathrm{\Psi }}_{\infty }\right|}^2-2\left(Q-4R{\left|{\mathrm{\Psi }}_{\infty }\right|}^2{a}^2\right){\left|{\mathrm{\Psi }}_{\infty }\right|}^2{v}_s^2,$$E52

$${\tilde{\rho }}_{\mathrm{sf}}\left(a\right)=\frac{C}{m}{\left|{\mathrm{\Psi }}_{\infty }\right|}^2{a}^2-2Q{\left|{\mathrm{\Psi }}_{\infty }\right|}^2{a}^2{v}_s^2.$$E53

From the above expression for $\mathcal{L}\left(a,\theta \right)$, we can now find the equations of motion for the relative amplitude field $a\left(x\right)$ and the phase field $\theta \left(x\right)$:

$$\frac{\partial }{\partial t}\left(\frac{\partial \mathcal{L}}{\partial \left({\partial }_{t}a\right)}\right)+\frac{\partial }{\partial x}\left(\frac{\partial \mathcal{L}}{\partial \left({\partial }_{x}a\right)}\right)=\frac{\partial \mathcal{L}}{\partial a},$$E54

$$\frac{\partial }{\partial t}\left(\frac{\partial \mathcal{L}}{\partial \left({\partial }_t\theta \right)}\right)+\frac{\partial }{\partial x}\left(\frac{\partial \mathcal{L}}{\partial \left({\partial }_x\theta \right)}\right)=\frac{\partial \mathcal{L}}{\partial \theta }.$$E55

Starting with the equation for the phase field, we easily find:

$$\frac{\partial }{\partial x}\left(\kappa \left(a\right)a^2{v}_s-{\tilde{\rho }}_{\mathrm{sf}}\left(a\right)\frac{\partial \theta }{\partial x}\right)=0$$E56

$$\iff \frac{\partial \theta }{\partial x}=\frac{\kappa \left(a\right)a^2{v}_s+\alpha }{{\tilde{\rho }}_{\mathrm{sf}}\left(a\right)}.$$E57

The integration constant $\alpha$ can be determined through the boundary condition for a dark soliton:

$$\frac{\partial \theta }{\partial x}\to 0\text{for}x\to \pm \infty .$$E58

which yields $\alpha =-v_s{\kappa }_{\infty }$ with ${\kappa }_{\infty }=\kappa \left(a_{\infty }\right)$ and thus

$$\frac{\partial \theta }{\partial x}=\frac{v_s}{{\tilde{\rho }}_{\mathrm{sf}}\left(a\right)}\left(\kappa \left(a\right)a^2-{\kappa }_{\infty }\right).$$E59

If we set $\theta \left(-\infty \right)=0$, the phase profile of the superfluid is given by

$$\theta \left(x\right)=v_s{\int }_{-\infty }^x\frac{\kappa \left(a\left(x^{\prime }\right)\right)a^2\left(x^{\prime }\right)-{\kappa }_{\infty }}{{\tilde{\rho }}_{\mathrm{sf}}\left(a\left(x^{\prime }\right)\right)}d{x}^{\prime }.$$E60

Next, we derive the equation of motion for $a\left(x\right)$:

$$\frac{\partial }{\partial x}\left(-{\tilde{\rho }}_{\mathrm{qp}}\left(a\right)\frac{\partial a}{\partial x}\right)=\frac{\partial }{\partial a}\left(\kappa \left(a\right)a^2\right)v_s\frac{\partial \theta }{\partial x}-\frac{\partial {\mathrm{\Omega }}_s}{\partial a}-\frac{1}{2}\frac{\partial {\tilde{\rho }}_{\mathrm{qp}}}{\partial a}{\left(\frac{\partial a}{\partial x}\right)}^2-\frac{1}{2}\frac{\partial {\tilde{\rho }}_{\mathrm{sf}}}{\partial a}{\left(\frac{\partial \theta }{\partial x}\right)}^2.$$E61

Inserting the solution for the derivative of the phase field (59) and defining

$$X\left(a\right)={\mathrm{\Omega }}_s\left(a\right)-{\mathrm{\Omega }}_s\left(a_{\infty }\right),$$E62

$$Y\left(a\right)=\frac{{\left(\kappa \left(a\right)a^2-{\kappa }_{\infty }\right)}^2}{2{\tilde{\rho }}_{\mathrm{sf}}\left(a\right)},$$E63

we find

$$\frac{1}{2}\frac{\partial {\tilde{\rho }}_{\mathrm{qp}}}{\partial a}{\left(\frac{\partial a}{\partial x}\right)}^2+{\tilde{\rho }}_{\mathrm{qp}}\left(a\right)\frac{{\partial }^{2}a}{\partial x^2}=\frac{\partial }{\partial a}\left(X\left(a\right)-v_s^{2}Y\left(a\right)\right).$$E64

While the above equation does not allow for a straightforward solution for $a$ as a function of the position $x$, it can be solved for $x$ as a function of $a$ instead. Using the boundary conditions for a dark soliton

$${\left.\frac{\partial a}{\partial x}\right|}_{x\to \pm \infty }=0\text{and}{\left.a\left(x\right)\right|}_{x\to \pm \infty }=1,$$E65

we find that (64) can be integrated, yielding:

$$\frac{1}{2}{\tilde{\rho }}_{\mathrm{qp}}\left(a\right){\left(\frac{\partial a}{\partial x}\right)}^2=X\left(a\right)-v_s^{2}Y\left(a\right),$$E66

$$\iff {\left(\frac{\partial x}{\partial a}\right)}^2=\frac{1}{2}\frac{{\tilde{\rho }}_{\mathrm{qp}}\left(a\right)}{X\left(a\right)-v_s^{2}Y\left(a\right)},$$E67

$$\iff x=\pm \frac{1}{\sqrt{2}}{\int }_{a_0}^a\frac{\sqrt{{\tilde{\rho }}_{\mathrm{qp}}\left(a^{\prime }\right)}}{\sqrt{X\left(a^{\prime }\right)-v_s^{2}Y\left(a^{\prime }\right)}}d{a}^{\prime }.$$E68

Here, $a_0=a\left(x=0\right)$ is the relative amplitude at the center of the soliton, which is found as the solution of

$$X\left(a_0\right)-v_s^{2}Y\left(a_0\right)=0.$$E69

For given values of the interaction parameter ${\left(k_F{a}_s\right)}^{-1}$, the temperature $T/T_F$, the imbalance chemical potential $\zeta$, and the soliton velocity $v_s$, formulae (60) and (68) allow us to calculate the complete pair field profile of the dark soliton. For example, the soliton density and phase profiles in Figure 4 were calculated using the above expressions.

5.3. Dark solitons in imbalanced Fermi gases

The dark soliton solution derived in the previous section has been employed in the description of various soliton phenomena in superfluid Fermi gases. For instance, adding a small two-dimensional perturbation to the exact 1D solution allows for a description of the snake instability mechanism [25], which makes the soliton decay into vortices if the radial width of the system is too large [23, 26]. We have also studied collisions between dark solitons by numerically evolving two counter-propagating 1D solitons in time [27]. As an example of an application, we will give a short description of the influence of spin-imbalance on dark solitons, a topic that was studied in detail in [18].

In ultracold Fermi gases, the amount of atoms in each spin population can be tuned experimentally, allowing for the possibility of having unequal amounts of spin-up and spin-down particles [28, 29]. In that case, when particles of different spin type pair up and form a superfluid state, an excess of unpaired particles will remain in the normal state, which in turn can have interesting effects on other phenomena in the system, including dark solitons. In the context of the EFT, we control the population imbalance by setting the value of the imbalance chemical potential $\zeta$, defined in (18). Figure 5a and b shows respectively the fermion particle density $n\left(x\right)$ and spin-population density difference $\mathit{\delta n}\left(x\right)$ (both with respect to the bulk density $n_{\infty }$) along a stationary dark soliton for ${\left(k_F{a}_s\right)}^{-1}=0$ (unitarity), $T=0.1T_F$, and for different values of $\zeta$. The density and density difference profiles are calculated using formulas (22) and (23) in a mean-field local density approximation. From the left figure, we observe that as we raise the imbalance chemical potential, the fermion density at the soliton center increases and the soliton broadens. However, we also know that, for a stationary dark soliton, the pair density at the center is always zero (as shown in the upper left panel of Figure 4), which means that the particles filling up the soliton are unpaired particles. This is confirmed by the right figure, which shows that the density difference between spin-up and spin-down particles in the soliton center increases with $\zeta$. The same effects are observed across the whole BEC-BCS crossover.

Figure 5.

Fermion density (left figure) and density difference (right figure) profiles of a dark soliton for ${\left(k_F{a}_S\right)}^{-1}=0$ at temperature $T/T_F=0.1$, for different values of the imbalance chemical potential $\zeta$. The densities are given with respect to the bulk density $n_{\infty }$.

As the imbalance between the spin components in the Fermi gas increases, so does the amount of unpaired particles that cannot participate in the superfluid state of pairs. While some of these normal state particles can coexist with the pair condensate as a thermal gas, it is energetically favorable for the remaining excess to be spatially separated from the superfluid. In this context, the soliton dip is a very suitable location to accommodate the excess particles and consequently fills up with an increasing amount of unpaired particles as the imbalance gets higher. Also, the broadening of the soliton with increasing imbalance might be a way of providing the system with more space to store the excess component. The fact that a dark soliton in an imbalanced superfluid Fermi gas has to drag along additional particles changes its effective mass, which in turn influences its general dynamical properties [18]. Moreover, since a soliton plane provides more space to accommodate the excess component than a vortex core, the presence of spin imbalance has been found to stabilize dark solitons with respect to the snake instability [25].

6. Application 2: the vortex structure

As a second application, the time-independent version of the theory is considered in order to derive the stable vortex structure. For the description of the vortex, the quantum velocity field $\mathbf{v}$ will be used, defined as:

$$\mathbf{v}=\frac{\mathrm{\hslash }}{m}{\mathbf{\nabla }}_{\mathbf{x}}\theta ,$$E70

where $\theta$ is the phase field from the hydrodynamical description (44). In the time-independent case, the action (30) reduces to the free energy (times the inverse temperature), which is given by:

$$F=\int \mathrm{d}\mathbf{r}\mathcal{F}\left(a,{\mathbf{\nabla }}_{\mathbf{x}}a,\mathbf{x}\right)\text{with}\mathcal{F}\left(a,{\mathbf{\nabla }}_{\mathbf{x}}a,\mathbf{x}\right)=X\left(a\right)+\frac{1}{8}{\rho }_{\mathrm{sf}}\left(a\right)v^2\left(\mathbf{x}\right)+\frac{1}{2}{\rho }_{\mathrm{qp}}\left(a\right){\left|{\mathbf{\nabla }}_{\mathbf{x}}a\right|}^2.$$E71

The free energy was written in a more compact³ form using the hydrodynamical description (48), (49), (62) and (70). As an application of the effective field theory, the general structure of a superfluid vortex will be numerically determined and compared with the commonly used variational hyperbolic tangent. A more detailed description on vortices in superfluids and their behavior can be found in [30].

6.1. What is a vortex?

Both in the classical and the quantum sense, a vortex is defined as a line in the fluid around which there is a circulating flow. In order to quantify this rotation around an axis, the circulation $\kappa$ is defined as:

$$\kappa ={\oint }_{\gamma }\mathbf{v}\left(\mathbf{r}\right)\cdot \mathrm{d}\mathbf{s},$$E72

where $\gamma$ is a closed contour and $\mathbf{v}$ the superfluid velocity field (70). A distinct feature of superfluids⁴ is that the circulation $\kappa$ is only allowed to take on values which are integer multiples of the circulation quantum $h/m$. In superfluids, circulation is always carried by quantized vortices.

This quantization of the circulation can be derived using the definition of the velocity field (70). Upon substitution, the circulation (72) can be written as:

$$\kappa =\frac{\mathrm{\hslash }}{m}{\oint }_{\gamma }{\mathbf{\nabla }}_{\mathbf{x}}\theta \cdot \mathrm{d}\mathbf{s}=n\frac{h}{m}\text{with}n\in \mathrm{\mathbb{Z}},$$E73

where the gradient theorem was used together with the fact that the phase field $\theta$ is a periodic function (period $2\pi$).

As the bulk superfluid itself is irrotational, any loop with nonzero circulation must encircle a node in the superfluid order parameter. As a consequence, the superfluid pair density must go to zero along the entire vortex line, resulting in a vortex “core” region with a radius comparable to the healing length. Important to note is that vortices of a single circulation quantum are energetically more favorable than multiply quantized vortices in a homogeneous condensate (which is the type of condensate that will be considered in this chapter) [9]. For the remainder of this application, only singly quantized vortices will thus be studied.

6.2. About the structure of a quantum vortex

The most natural coordinate system to describe vortices are the polar coordinates $\mathbf{x}=\left(r,\varphi \right)$. The origin of the polar coordinates will be chosen in the center of the vortex (at the point where the superfluid density reaches zero). In order to derive the vortex structure, a set of boundary conditions is required. In the radial direction, the boundary conditions are then given by⁵:

$$a\left(r\to 0\right)=0\text{and}a\left(r\to \infty \right)=1,$$E74

meaning that the superfluid density relaxes to the bulk value away from the vortex.

We factorize the amplitude function in a radial and an angular part⁶:

$$a\left(r,\varphi \right)=f\left(r\right)\mathrm{\Phi }\left(\varphi \right).$$E75

Since the structure is periodic, the general solution for $\mathrm{\Phi }\left(\varphi \right)$ is thus given by:

$$\mathrm{\Phi }\left(\varphi \right)=\sum _{n=-\infty }^{\infty }{a}_n{e}^{\mathrm{i}\mathit{n\varphi }},$$E76

leading to a basis of angular modes for the vortex structure. In order to find the lowest energy state, one usually restricts the problem to one of the many possible modes:

$$\mathrm{\Phi }\left(\varphi \right)=e^{\mathrm{i}\mathit{n\varphi }}\text{with}n\in \mathrm{\mathbb{Z}},$$E77

which results in the velocity field and circulation (using (70) and (72)) for a single mode given by:

$$\mathbf{v}\left(r\right)=n\frac{\mathrm{\hslash }}{\mathit{mr}}{\mathbf{e}}_{\varphi }\Rightarrow \kappa =n\frac{h}{m},$$E78

where the velocity field diverges in the point where the superfluid vanishes. It was noted before that for our case, the most energetic vortex states are those with the least circulation quanta. Since the object of interest is the vortex structure with a minimal free energy, the value of $n$ will be restricted to $n=\pm 1$. The state with $n=1$ is known as the “vortex,” where the state with $n=-1$ is known as the “anti-vortex.” This means that the vortex velocity field is given by⁷:

$$\mathbf{v}\left(r\right)=\pm \frac{\mathrm{\hslash }}{\mathit{mr}}{\mathbf{e}}_{\varphi },$$E79

where the “+” sign is for vortices and the “-” sign for anti-vortices.

Currently, there is no analytical solution available for the full vortex structure $f\left(r\right)$. Calculations including vortices are therefore either done numerically (for the exact structure) or variationally. One way to numerically find the minimal structure is by writing down the equations of motion (the Euler-Lagrange equations for the free energy (71)), which is analogous to what was done for the soliton in the previous section. Directly solving the equations of motion, however, is a numerical challenge due to the divergence of the velocity field in the center of the vortex. A second numerical method is briefly discussed further on. The disadvantage of the full numerical approach is that it takes time. As an alternative, it is possible to work with a variational model. By working with a variational model, it is possible to retain a fair amount of accuracy while gaining several orders of magnitude in computational speed. The usage of variational models is discussed in the next subsection. A disadvantage of using variational models is however that a certain structure is proposed, meaning the variational guess can be wrong in certain situations. When using variational models, one should consequently always check the validity of the model and the range of application.

6.3. A variational model for the vortex core

In order to speed up the vortex calculations, a variational model can be used to describe the vortex structure. First of all, the variational model should meet the required boundary conditions (74). Second, the variational model should contain the necessary information to describe the vortex physics. For example, in liquid helium, the vortex core sizes are of the order of nanometers [33], meaning that the vortex core structure will not play a prominent role in the vortex physics; in this case, a simple hollow cylinder is already a good variational model for the vortex core. For vortices in ultracold gases on the other hand, the vortex core size is of the order of micrometers [34], meaning that its structure becomes important⁸; a simple cylindric hole will no longer capture the entire vortex physics. In order to provide a more detailed description, different variational models are available [9, 30].

The variational model that will be discussed here is the hyperbolic tangent model:

$$f\left(r\right)=\text{tanh}\left(\frac{r}{\sqrt{2}\xi }\right),$$E80

where the quantity $\xi$ is defined as the healing length. The hyperbolic tangent (80) is the exact solution of the Gross-Pitaevskii equation in 1D for a condensate with a hard wall boundary [9, 35]. Since the variational model describes the healing from a hole in the condensate, it is expected that this model will also sufficiently describe the vortex physics. The merit of the presented effective field theory in Section 4 is that an analytical solution can be derived for the vortex healing length $\xi$; this will be done in the remainder of this subsection. Using the definitions (71), the free energy of the variational vortex structure is given by:

$$F=\underset{0}{\overset{\infty }{\int }}r\mathrm{d}r\left[X\left(a\right)+\frac{A}{2r^2}{\text{tanh}}^2\left(\frac{r}{\sqrt{2}\xi }\right)+\frac{{\rho }_{\mathrm{qp}}\left(\text{tanh}\left(\frac{r}{\sqrt{2}\xi }\right)\right)}{4{\xi }^2{\text{cosh}}^4\left(\frac{r}{\sqrt{2}\xi }\right)}\right],$$E81

where the value of the constant $A$ is defined as:

$$A=2C{\left|{\mathrm{\Delta }}_{\infty }\right|}^2.$$E82

The second term in the integrand of (81) causes a divergence, since

$$\underset{R\to \infty }{\text{lim}}\frac{A}{2}\underset{0}{\overset{R}{\int }}\frac{1}{r}{\text{tanh}}^2\left(\frac{r}{\sqrt{2}\xi }\right)\mathrm{d}r\propto \text{log}\left(R\right)$$E83

diverges logarithmically with increasing radius of the integration domain. The physical reason is clear: the velocity profile of a vortex decays as $1/r$, so that the kinetic energy of the superflow will grow as the logarithm of the container size. However, the derivative with respect to $\xi$ of this kinetic energy of the superflow does not diverge. This can be seen by first switching to a dimensionless variable $x=r/\xi$:

$$F={\xi }^2\underset{0}{\overset{\infty }{\int }}\mathit{xX}\left(a\right)\mathrm{d}x+\frac{A}{2}\underset{R\to \infty }{\text{lim}}\left(\underset{0}{\overset{R/\xi }{\int }}{\text{tanh}}^2\left(\frac{x}{\sqrt{2}}\right)\frac{\mathrm{d}x}{x}\right)+\underset{0}{\overset{\infty }{\int }}\frac{{\rho }_{\mathit{qp}}\left(\text{tanh}\left(x/\sqrt{2}\right)\right)}{4{\text{cosh}}^4\left(x/\sqrt{2}\right)}\mathrm{d}x.$$E84

The last term no longer contains a dependency on $\xi$, so its derivative with respect to $\xi$ vanishes. We obtain

$$\frac{\mathrm{d}F}{\mathrm{d}\xi }=2\xi \underset{0}{\overset{\infty }{\int }}\mathit{xX}\left(a\right)\mathrm{d}x+\frac{A}{2}\underset{R\to \infty }{\text{lim}}\frac{\mathrm{d}}{\mathrm{d}\xi }\left[\underset{0}{\overset{R/\xi }{\int }}{\text{tanh}}^2\left(\frac{x}{\sqrt{2}}\right)\frac{\mathrm{d}x}{x}\right]$$E85

The remaining derivative now acts on the boundary of the integration domain. Applying

$$\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}\xi }\underset{0}{\overset{R/\xi }{\int }}g\left(x\right)\mathrm{d}x=\underset{\mathrm{\Delta }\xi \to 0}{\text{lim}}\frac{1}{\mathrm{\Delta }\xi }\left(\underset{0}{\overset{R/\left(\xi +\mathrm{\Delta }\xi \right)}{\int }}g\left(x\right)\mathrm{d}x-\underset{0}{\overset{R/\xi }{\int }}g\left(x\right)\mathrm{d}x\right)\\ =\underset{\mathrm{\Delta }\xi \to 0}{\text{lim}}\frac{1}{\mathrm{\Delta }\xi }\left(\frac{R}{\xi +\mathrm{\Delta }\xi }-\frac{R}{\xi }\right)g\left(\frac{R}{\xi }\right)=-\frac{R}{{\xi }^2}g\left(\frac{R}{\xi }\right)\end{array}$$E86

to (85) we get

$$\frac{\mathrm{d}F}{\mathrm{d}\xi }=2\xi \underset{0}{\overset{\infty }{\int }}\mathit{xX}\left(a\right)\mathrm{d}x-\frac{A}{2\xi }\underset{R\to \infty }{\text{lim}}\left[{\text{tanh}}^2\left(\frac{R}{\sqrt{2}\xi }\right)\right],$$E87

so that

$$\frac{\mathrm{d}F}{\mathrm{d}\xi }=0\iff 2\xi \underset{0}{\overset{\infty }{\int }}\mathit{xX}\left(a\right)\mathrm{d}x=\frac{A}{2\xi }.$$E88

With $A=2C{\left|{\mathrm{\Delta }}_{\infty }\right|}^2$ as above, and

$$B=\underset{0}{\overset{\infty }{\int }}\mathit{xX}\left(a\right)\mathrm{d}x$$E89

we find a closed form result

$$\xi =\frac{1}{2}\sqrt{\frac{A}{B}}.$$E90

The formula for the healing length (90) can also be plotted, this is done in Figure 6. In both limits, the healing length shows to be in a good agreement with the exact limits.

Figure 6.

The vortex variational healing length (90) throughout the BEC-BCS crossover for the case $\beta =100$ and $\zeta =0$. The dotted lines yield the exact solutions in the deep BEC [36] and BCS [37] limits. This plot made using the same data as in [38].

6.4. Comparison to the exact (numerical) solution

In order to check the validity of the variational model (80) (and thus the results it produces) is, the variational structure should be compared with the exact vortex structure. This exact vortex structure is easily obtained by a direct minimization of the free energy functional (71). As mentioned before, the direct minimization of the free energy is more suitable for the calculation of the vortex structures; the reason why this method is preferential lies in the fact that the velocity field diverges for $r\to 0$. The divergence of the velocity field in the origin will be strongly pronounced when solving the equations of motion. While in the case of a direct minimization, the same divergence will have less impact on the solution.

The numerical method that was used in order to determine the exact vortex structure for a given set of parameters $\left(\beta ,a_s,\zeta \right)$ is discussed in full detail in [38]. In a nutshell, this method comes down to making a discretized version of the vortex structure: $\left\{f_1,f_2,\dots ,f_N\right\}$, where $f_1=0$ and $f_N=1$ due to the boundary conditions. During the minimization procedure, a program runs through the list of points $\left\{f_n|n\in \left\{2,3,\dots ,N-1\right\}\right\}$, where it suggests a (random) new value; if the new value results in a lower energy, it is accepted as the new value of the vortex structure. The minimization program continues to run until a certain tolerance is reached and the structure is not changing any more.

Once the exact structure is obtained, it can be analyzed and compared to the variational vortex structure. As an example, we can look at the relative difference in the free energy throughout the BEC-BCS crossover for different temperatures and polarizations. From the plots shown in Figure 7, it can be seen that the difference in free energy is around the order of 1%; this seems to suggest that the variational guess is a good guess.⁹ Moreover, the results in Figure 7 allow to provide an error bar on energy calculations using the variational structure. This error bar on the energy is useful for example when making phase diagrams including vortex structures. In order to be sure whether the variational model is indeed a good description of the vortex hole, other parameters were also tested and discussed in [38]. The conclusion from the numerical analysis was that the variational model is indeed a good fit for describing the vortex structure.

Figure 7.

The relative energy difference between the exact and variational solutions throughout the BEC-BCS crossover for different values of the temperature $\beta$ and polarization $\zeta$. These results were also shown and discussed in [38].

7. Concluding section

In this chapter, an effective field theory for the description of dilute fermionic superfluids was derived. The main advantages of an effective field theory are the gain in computational speed and the fact that it allows analytic solutions for dark solitons and the variational healing length of the vortex structure. Both the gain in computational speed and the availability of an analytic starting point contribute to the possibility to study several soliton/vortex phenomena throughout the entire BEC-BCS crossover at finite temperatures $\beta$ for a given polarization $\zeta$ within a reasonable computational time span.

On the subject of soliton dynamics, we specifically looked at 1D dark solitons, for which an exact analytical solution was derived. Using this solution, the effect of spin-imbalance on the soliton properties was studied, revealing that the unpaired particles of the excess component mainly reside inside the soliton core. Additionally, the EFT has also been employed in the study of the snake instability of dark solitons [25] and the dynamics of dark soliton collisions [27] in imbalanced superfluid Fermi gases.

For vortices, the structure of a vortex was studied, for which unfortunately no analytical solution is available at the moment. Using a variational model, an analytical solution for the vortex healing length was derived. The variational model was compared with the exact solution. From this analysis, the variational model was found to be a good fit for the exact vortex structure. Other EFT research on vortices includes the behavior of vortices in multiband systems [39] and the study of the “vortex charge” [40].

Acknowledgments

The authors gratefully acknowledge the financial support provided by the Fund for Scientific Research Flanders (FWO), through the FWO project: G042915 N (Superfluidity and superconductivity in multi-component quantum condensates). One of us (N.V.) acknowledges a Ph.D. fellowship of the University of Antwerp (2014BAPDOCPROEX167). One of us (W.V.A.) acknowledges a Ph.D. fellowship from the FWO (1123317 N). We also acknowledge financial support from the Research Fund (BOF-GOA) of the University of Antwerp.