Diversified vortex phase diagram for a rotating trapped two-band Fermi gas in the BCS-BEC crossover

Quantum gases constitute a remarkable testing ground for the theory of the superfluid state and its various macroscopic excitations, such as vortices and solitons. Vortices are stabilized in an atomic gas by rotating the gas, just as vortices in superconductors can be stabilized by a magnetic field: the Coriolis force acts on a particle in a rotating frame of reference in the same way as the Lorentz force acts on a charged particle in a magnetic field [1]. Consequently, vortices and many-vortex states in rotating trapped superfluid Fermi gases have become a subject of an intense experimental [2–5] and theoretical [6–15] research during last two decades.

Since recently, this line of research has become of particular interest in multicomponent quantum systems. For example, multiband superconductors, such as MgB₂, were intensely studied both experimentally and theoretically over the last decade. However, interaction parameters can hardly be tuned in solid-state systems. Contrary to superconductors, interactions in ultracold Fermi gases can be controlled, and a broad range of regimes from BCS to BEC can be realized. This stimulated a theoretical interest to multiband quantum gases in anticipation of future experiments [16–18]. Recently, a two-band Fermi superfluid has been successfully created in an ultracold gas of ¹⁷³Yb atoms [19, 20] using the orbital Feshbach resonance predicted in [21] (see also [22]). This makes a theoretical investigation of different phenomena in multiband Fermi gases (vortices, solitons, etc) timely and important. Vortices in multiband superconducting and superfluid systems are particularly interesting due to a rich variety of observable phenomena, such as fractional vortex states that occur when winding numbers in different band-components of the condensate are not equal. Fractional vortices in multiband superconductors have been investigated fairly thoroughly [23–25]. Concerning the BCS–BEC crossover, there was a widely spread opinion during a long time that superconductors cannot be realized away from the BCS regime. However recently the BCS–BEC crossover has been successfully reached in superconducting FeSe [26]. Moreover, vortex matter in multiband superconductors in the BCS-BEC crossover regime has also attracted much attention [27]. Vortex states in multiband quantum atomic gases have been studied to a far lesser extent.

Therefore, in this work, the subject of our interest are fractional vortices in two-band Fermi gases of ultracold atoms in the BCS–BEC crossover. Although there is some similarity between superconductors and condensed atomic Fermi gases, the analogy is not complete. For one difference, cold gases certainly require an independent treatment by specific methods suitable in the entire BCS–BEC crossover range.

Recently, the stability of different vortex states in a rotating trapped one-band Fermi gas has been theoretically studied in [14, 15] using, respectively, the coarse-grained Bogoliubov–de Gennes (BdG) theory [28] (first applied to atomic Fermi gases in [29]) and the recently developed finite-temperature effective field theory (EFT) [30, 31]. The finite-temperature EFT results agree with the results of the BdG theory and experiment for different manifestations: collective excitations and vortices [15, 31], and solitons [32, 33]. The finite-temperature EFT is aimed to find analytic results whenever possible. For example, for dark solitons in condensed Fermi gases, the finite-temperature EFT provides exact analytic solutions of the soliton equation of motion [32], while the BdG equations for the same problem have been solved only numerically.

Besides our studies, there are several modifications of the EFT of condensed Fermi gases described in different publications and related to different ranges of external parameters (e.g., temperature and scattering length). They are developed either for the close vicinity to the critical temperature [34] or for the case T = 0 (e.g., [35–37]). As analyzed in [31], the present finite-temperature EFT agrees with preceding works in all these limiting cases.

At present, experimental data on vortices in two-band superfluid atomic Fermi gases are still lacking, in spite of the expected new physics stemming from the interband interactions in such a system. We report the first such theoretical study, to pave the way for future experiments with two- and multiband atomic Fermi superfluids. This work builds on the research performed in [15], with an extension to two-band fermionic systems. The main goal of the present investigation is to reveal novel vortex phenomena which can appear in a two-band Fermi gas, that are arguably easy to verify experimentally. More specifically, we study the evolution of equilibrium vortex states when varying the temperature and the interband coupling strength, as well as the frequency of rotation, to identify regions of stability for fractional vortices, clusters of non-composite vortices, and multivortex states. Two variants are considered: (1) the 'canonical ensemble' case when numbers of particles in each band are fixed separately, and (2) the 'grand canonical ensemble' case, when the numbers of particles are determined from the common chemical potential.

The paper is organized as follows. In section 2, the used method and approximations are described. In section 3, the parameters of state and the vortex phase diagrams are analyzed in detail. Our results are summarized in section 4.

We apply the finite-temperature EFT developed in [31] to vortices and many-vortex states in a two-band Fermi gas with s-wave pairing. The fermion system is confined to a cylindrically symmetric parabolic trapping potential with confinement frequency ω₀. This setup can be relevant for ultracold fermionic atoms in elongated traps, as, e.g., in [2–5]. The stabilization of vortices is achieved by rotating the fermionic system with an angular frequency ω. The rotation is incorporated in the EFT in the same manner as in [15].

The details of the finite-temperature EFT with a discussion of the approximations used and the range of applicability of the method are readily available in [31]. The incorporation of rotation in the EFT and the extension of the formalism to a two-band rotating Fermi gas have been already described in [15]. Consequently, in this section we reproduce the formalism only briefly, because detailed derivations, discussions and proofs can be found in [15, 31].

The treatment of rotating trapped Fermi gases is performed within the path-integral formalism in the space of anticommuting fermion fields, see, e.g., [34, 35]. Throughout the treatment, we use units such that ℏ = 1, the fermion mass for a 'strong' band m₁ = 1/2, and the Fermi energy for a non-interacting fermion system E_F = 1. Note that there is a difference between the Fermi energy for a Fermi gas in bulk and in a trapping potential. In the present system of units, E_F is the global rather than a local parameter, i.e., E_F is not coordinate-dependent. For example, the Fermi energy for a gas trapped to a 3D parabolic confinement potential with confinement frequencies ω_x, ω_y, ω_z is ${E}_{F}\approx {\hslash }{(3{\omega }_{x}{\omega }_{y}{\omega }_{z}N)}^{1/3}$ when N is sufficiently large. The parameter having the dimensionality of the Fermi wave vector is formally determined as ${k}_{F}\equiv \sqrt{2{{mE}}_{F}}/{\hslash }$.

The thermodynamic quantities are calculated on the basis of the partition function,

$$\begin{eqnarray}&&{ \mathcal Z }\propto \int { \mathcal D }[\bar{\psi },\psi ]{{\rm{e}}}^{-S}.\end{eqnarray} \tag{ 1 }$$

Here, S is the action functional for a two-band fermionic system [31]:

$$\begin{eqnarray}&&S=\displaystyle \sum _{j=1,2}{S}_{0,j}+{\int }_{0}^{\beta }{\rm{d}}\tau \int {\rm{d}}{\bf{r}}\,U({\bf{r}},\tau ),\end{eqnarray} \tag{ 2 }$$

where $\beta =1/({k}_{B}T),T$ is the temperature, k_B is the Boltzmann constant, and S_{0, j} is the free-fermion action for the jth band,

$$\begin{eqnarray}&&{S}_{0,j}={\int }_{0}^{\beta }{\rm{d}}\tau \int {\rm{d}}{\bf{r}}\displaystyle \sum _{\sigma =\uparrow ,\downarrow }{\bar{\psi }}_{\sigma ,j}\left(\displaystyle \frac{\partial }{\partial \tau }+{H}_{\sigma ,j}\right){\psi }_{\sigma ,j},\end{eqnarray} \tag{ 3 }$$

with the one-particle Hamiltonian in the rotating frame of reference, which allows for independent populations in different bands and different spins σ, given by

$$\begin{eqnarray}&&{H}_{\sigma ,j}=-\displaystyle \frac{{({\rm{\nabla }}-{\rm{i}}{{\bf{A}}}_{j}({\bf{r}}))}^{2}}{2{m}_{j}}+\displaystyle \frac{{m}_{j}({\omega }_{0}^{2}-{\omega }^{2})}{2}({x}^{2}+{y}^{2})-{\mu }_{\sigma ,j},\end{eqnarray} \tag{ 4 }$$

where ${{\bf{A}}}_{j}({\bf{r}})$ is the rotational vector potential ${{\bf{A}}}_{j}({\bf{r}})={m}_{j}[{\boldsymbol{\omega }}\times {\bf{r}}]$, and ${\boldsymbol{\omega }}\equiv \omega {{\bf{e}}}_{z}$ is the rotation vector. The fermion–fermion interaction $U({\bf{r}},\tau )$ for a two-band system assumes both intraband and interband contact interactions,

$$\begin{eqnarray}U & = & \displaystyle \sum _{j=1,2}{g}_{{jj}}{\bar{\psi }}_{\uparrow ,j}{\bar{\psi }}_{\downarrow ,j}{\psi }_{\downarrow ,j}{\psi }_{\uparrow ,j}+{g}_{12}^{(a)}({\bar{\psi }}_{\uparrow ,1}{\psi }_{\uparrow ,1}{\bar{\psi }}_{\downarrow ,2}{\psi }_{\downarrow ,2}+{\bar{\psi }}_{\downarrow ,1}{\psi }_{\downarrow ,1}{\bar{\psi }}_{\uparrow ,2}{\psi }_{\uparrow ,2})\\ & & +{g}_{12}^{(p)}({\bar{\psi }}_{\uparrow ,1}{\psi }_{\uparrow ,1}{\bar{\psi }}_{\uparrow ,2}{\psi }_{\uparrow ,2}+{\bar{\psi }}_{\downarrow ,1}{\psi }_{\downarrow ,1}{\bar{\psi }}_{\downarrow ,2}{\psi }_{\downarrow ,2}),\end{eqnarray} \tag{ 5 }$$

with intraband coupling constants g_jj(j = 1, 2) and interband coupling constants ${g}_{12}^{(a)}$ and ${g}_{12}^{(p)}$ which describe interband contact interactions of density fluctuations with antiparallel and parallel spins, respectively.

The finite-temperature EFT action for a one- and two-band system of interacting fermions with s-wave pairing has been formulated in [31]. Here, details of the derivation of the effective field action are given in the appendix. After the Hubbard–Stratonovich transformation using bosonic pair fields Ψ₁, Ψ₂, the integration over fermion fields and the gradient expansion (assuming slow variation of the pair field in time and space), the resulting action functional for a two-band system ${S}_{{\rm{EFT}}}^{(2{b})}$ takes the form:

$$\begin{eqnarray}&&{S}_{{\rm{EFT}}}^{(2{b})}=\displaystyle \sum _{j=1,2}{S}_{{\rm{EFT}}}^{(j)}-{\int }_{0}^{\beta }{\rm{d}}\tau \int {\rm{d}}{\bf{r}}\,\displaystyle \frac{\sqrt{{m}_{1}{m}_{2}}}{4\pi }\gamma ({\bar{{\rm{\Psi }}}}_{1}{{\rm{\Psi }}}_{2}+{\bar{{\rm{\Psi }}}}_{2}{{\rm{\Psi }}}_{1}),\end{eqnarray} \tag{ 6 }$$

where γ is the strength of the interband coupling expressed through the interband scattering lengths,

$$\begin{eqnarray}&&\gamma =2\left(\displaystyle \frac{1}{{a}_{s,12}^{(a)}}-\displaystyle \frac{1}{{a}_{s,12}^{(p)}}\right).\end{eqnarray} \tag{ 7 }$$

The one-band effective field action ${S}_{{\rm{EFT}}}^{(j)}$ for the jth band is determined by:

$$\begin{eqnarray}{S}_{{\rm{EFT}}}^{(j)} & = & {\displaystyle \int }_{0}^{\beta }{\rm{d}}\tau \displaystyle \int {\rm{d}}{\bf{r}}\left[{{\rm{\Omega }}}_{s,j}({| {{\rm{\Psi }}}_{j}| }^{2})+\displaystyle \frac{{D}_{j}}{2}\left({\bar{{\rm{\Psi }}}}_{j}\displaystyle \frac{\partial {{\rm{\Psi }}}_{j}}{\partial \tau }-\displaystyle \frac{\partial {\bar{{\rm{\Psi }}}}_{j}}{\partial \tau }{{\rm{\Psi }}}_{j}\right)\right.\\ & & +{Q}_{j}\displaystyle \frac{\partial {\bar{{\rm{\Psi }}}}_{j}}{\partial \tau }\displaystyle \frac{\partial {{\rm{\Psi }}}_{j}}{\partial \tau }-\displaystyle \frac{{R}_{j}}{2{| {{\rm{\Psi }}}_{j}| }^{2}}{\left(\displaystyle \frac{\partial ({| {{\rm{\Psi }}}_{j}| }^{2})}{\partial \tau }\right)}^{2}+{C}_{j}({{\rm{\nabla }}}_{{\bf{r}}}{\bar{{\rm{\Psi }}}}_{j}\cdot {{\rm{\nabla }}}_{{\bf{r}}}{{\rm{\Psi }}}_{j})-{E}_{j}{({{\rm{\nabla }}}_{{\bf{r}}}({| {{\rm{\Psi }}}_{j}| }^{2}))}^{2}\\ & & +{{\rm{i}}{D}}_{j}{{\bf{A}}}_{j}\cdot ({\bar{{\rm{\Psi }}}}_{j}{{\rm{\nabla }}}_{{\bf{r}}}{{\rm{\Psi }}}_{j}-{{\rm{\Psi }}}_{j}{{\rm{\nabla }}}_{{\bf{r}}}{\bar{{\rm{\Psi }}}}_{j})].\end{eqnarray} \tag{ 8 }$$

The coordinate-dependent thermodynamic potential ${{\rm{\Omega }}}_{s,j}({| {{\rm{\Psi }}}_{j}| }^{2})$ and the coefficients of the gradient expansion given in [15, 31] and explained in detail in the appendix. The rotation is incorporated in the effective action as described in [15]. It leads to the appearance of the linear gradient term ${{\rm{i}}{D}}_{j}{{\bf{A}}}_{j}\cdot ({\bar{{\rm{\Psi }}}}_{j}{{\rm{\nabla }}}_{{\bf{r}}}{{\rm{\Psi }}}_{j}-{{\rm{\Psi }}}_{j}{{\rm{\nabla }}}_{{\bf{r}}}{\bar{{\rm{\Psi }}}}_{j})$ in the effective bosonic action, which vanishes in the absence of rotation. As shown in [15], the obtained extension of the effective bosonic action to a rotating fermion system is in agreement with results of the functional renormalization group theory [38, 39].

3.1. Conditions for vortex stability

One of the most interesting consequences of rotation in a two-band Fermi system is the formation of fractional vortices. The fractional vortices can be energetically favorable only at rather small interband coupling strength, because the Josephson coupling energy (the last term in (6)) penalizes phase differences between the condensates in different bands. However, two-band Fermi gases with vanishingly small (and even zero) interband coupling can be prepared in practice—for example, when the two condensates are spatially separated. Therefore in the present work we study vortices in two-band Fermi gases with a particular attention to the range of small Josephson coupling strengths.

The stability of different vortex states is analyzed using the free energy ${F}^{(2b)}$, which is obtained from the effective action (8) assuming pair fields to be stationary (time-independent),

$$\begin{eqnarray}&&{F}^{(2b)}=\displaystyle \sum _{j=1,2}{F}^{(j)}-\int {\rm{d}}{\bf{r}}\,\displaystyle \frac{\sqrt{{m}_{1}{m}_{2}}}{4\pi }\gamma ({\bar{{\rm{\Psi }}}}_{1}{{\rm{\Psi }}}_{2}+{\bar{{\rm{\Psi }}}}_{2}{{\rm{\Psi }}}_{1}),\end{eqnarray} \tag{ 9 }$$

where ${F}^{(j)}$ is the contribution to the free energy from the jth band-component of the condensate:

$$\begin{eqnarray}{F}^{(j)} & = & \displaystyle \int {\rm{d}}{\bf{r}}[{{\rm{\Omega }}}_{s,j}({| {{\rm{\Psi }}}_{j}| }^{2})+{C}_{j}({{\rm{\nabla }}}_{{\bf{r}}}{\bar{{\rm{\Psi }}}}_{j}\cdot {{\rm{\nabla }}}_{{\bf{r}}}{{\rm{\Psi }}}_{j})\\ & & -{E}_{j}{({{\rm{\nabla }}}_{{\bf{r}}}({| {{\rm{\Psi }}}_{j}| }^{2}))}^{2}+{{\rm{i}}{D}}_{j}{{\bf{A}}}_{j}\cdot ({\bar{{\rm{\Psi }}}}_{j}{{\rm{\nabla }}}_{{\bf{r}}}{{\rm{\Psi }}}_{j}-{{\rm{\Psi }}}_{j}{{\rm{\nabla }}}_{{\bf{r}}}{\bar{{\rm{\Psi }}}}_{j})].\end{eqnarray} \tag{ 10 }$$

The pair field ${{\rm{\Psi }}}_{j}({\bf{r}})$ for a particular vortex state with n vortices in the jth band-component of the condensate is described here by the variational function, which is a product of the background pair field ${{\rm{\Delta }}}_{j}(r)$ and the vortex factor:

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{j}({\bf{r}})\equiv {{\rm{\Psi }}}_{j}({\bf{r}},\{{a}_{\nu ,j},{\theta }_{\nu ,j}\})={{\rm{\Delta }}}_{j}(r)\displaystyle \sum _{\nu =1}^{n}{a}_{\nu ,j}({\bf{r}}){{\rm{e}}}^{{\rm{i}}{\theta }_{\nu ,j}({\bf{r}})},\end{eqnarray} \tag{ 11 }$$

where ${a}_{\nu ,j}({\bf{r}})$ and ${\theta }_{\nu ,j}({\bf{r}})$ are, respectively, the relative amplitude and the phase for the νth vortex in the jth band component of the condensate. The phase for a vortex ${\theta }_{\nu ,j}({\bf{r}})$ coincides with the polar angle calculated in the frame of reference related to the vortex center. The vortex amplitudes in the present work, as well as in [15], are modeled by variational functions,

$$\begin{eqnarray}&&{a}_{\nu ,j}({\bf{r}})=\tanh \left(\displaystyle \frac{| {\bf{r}}-{{\bf{r}}}_{\nu ,j}| }{\sqrt{2}{\xi }_{j}}\right).\end{eqnarray} \tag{ 12 }$$

The parameter ξ_j has the sense of the vortex healing length. The healing lengths and the positions of vortex centers r_ν,j are variational parameters for vortex states. Their optimal values are found by the minimization of the free energy (9) after substituting the variational function (11) there. The conditions of stability for vortex states are determined here, as in [15], by considering the difference between the free energies with and without vortices:

$$\begin{eqnarray}&&\delta F\equiv F[\{{{\rm{\Psi }}}_{j}({\bf{r}})\}]-F[\{{{\rm{\Delta }}}_{j}(r)\}].\end{eqnarray} \tag{ 13 }$$

The parametric ranges of existence for vortex configurations are determined by comparison of free energies corresponding to different vortex states, including also the state without vortices, where ${{\rm{\Psi }}}_{j}({\bf{r}})={{\rm{\Delta }}}_{j}(r)$, and choosing a vortex state with the lowest free energy.

The background distribution of the pair field ${{\rm{\Delta }}}_{j}(r)$ without vortices is determined using normalization conditions for the fermion density. We consider here two possible constraints of the numbers of fermions in two bands. Under the first condition, total numbers of particles in each band are fixed. The alternative condition is the 'grand canonical' setting, when the populations in the two bands have been relaxed to values determined through equal chemical potentials in the two band-components of the condensate. At present, the ensemble with fixed numbers of particles seems to be more easily achieved in experiments with atomic Fermi gases than the grand canonical ensemble with equal chemical potentials. However, it is interesting to consider them both, in order to understand the influence of these constraints on vortex states.

In order to make clear the difference between the two aforesaid ensembles, we describe the normalization conditions for trapped two-band Fermi condensates in detail. Since we consider systems where the scale of the confinement potential is sufficiently large with respect to the vortex size, the normalization conditions are applied within the local density approximation through coordinate-dependent thermodynamic parameters. For a one-band Fermi gas trapped in a cylindrically symmetric confinement potential, the normalization condition in the local density approximation reads,

$$\begin{eqnarray}&&2\pi {\int }_{0}^{\infty }n(\beta ,\mu (r),{\rm{\Delta }}(r)){r}{\rm{d}}{r}=N,\end{eqnarray} \tag{ 14 }$$

where $n(\beta ,\mu (r),{\rm{\Delta }}(r))$ is the fermion density depending on the radial variable $r=\sqrt{{x}^{2}+{y}^{2}}$ through the background pair field ${\rm{\Delta }}(r)$ (we neglect here the feedback of vortices to the density normalization) and the chemical potential, according to (4),

$$\begin{eqnarray*}&&\mu (r)=\mu (0)-\displaystyle \frac{m}{2}({\omega }_{0}^{2}-{\omega }^{2}){r}^{2}.\end{eqnarray*}$$

The coordinate-dependent gap parameter is determined through the local gap equation [31],

$$\begin{eqnarray}&&\int \displaystyle \frac{{\rm{d}}{\bf{k}}}{{(2\pi )}^{3}}\left(\displaystyle \frac{\sinh \beta {E}_{{\bf{k}}}}{2{E}_{{\bf{k}}}(\cosh \beta {E}_{{\bf{k}}}+\cosh \beta \zeta )}-\displaystyle \frac{m}{{k}^{2}}\right)+\displaystyle \frac{m}{4\pi {a}_{s}}=0,\end{eqnarray} \tag{ 15 }$$

consistently with the normalization condition (14). Here, ${E}_{{\bf{k}}}=\sqrt{{\left(\tfrac{{k}^{2}}{2m}-\mu (r)\right)}^{2}+{\rm{\Delta }}{(r)}^{2}}$ is the Bogoliubov excitation energy, $\mu (r)=({\mu }_{\uparrow }+{\mu }_{\downarrow })/2$ and $\zeta =({\mu }_{\uparrow }-{\mu }_{\downarrow })/2$ are, respectively, the averaged chemical potential and the difference of chemical potentials for the 'spin up' and 'spin down' species.

The fermion density n entering normalization condition (14) can be calculated in different approximations, e.g., the mean-field or Gaussian pair fluctuation approximation, as in [40]. In the present work, we restrict the normalization of the density by the mean-field approximation, because account of fluctuations in the background parameters can lead only to a rescaling of phase diagrams, without changing their shape and sequences of areas of stability for vortex states. It should be noted that the gap equation has the same form (15) with and without account of fluctuations [35, 40].

For a two-band system, coordinate-dependent parameters of state are determined from two normalization conditions,

$$\begin{eqnarray}&&2\pi {\int }_{0}^{\infty }{n}_{j}(\beta ,{\mu }_{j}(r),{{\rm{\Delta }}}_{j}(r)){r}{\rm{d}}{r}={N}_{j}\quad (j=1,2)\end{eqnarray} \tag{ 16 }$$

together with the coupled set of two gap equations for the two-band Fermi gas given by equation (28) of [31],

$$\begin{eqnarray}&&\left[\displaystyle \int \displaystyle \frac{{\rm{d}}{\bf{k}}}{{(2\pi )}^{3}}\left(\displaystyle \frac{1}{2{E}_{{\bf{k}},j}}\displaystyle \frac{\sinh (\beta {E}_{{\bf{k}},j})}{\cosh (\beta {E}_{{\bf{k}},j})+\cosh (\beta {\zeta }_{j})}-\displaystyle \frac{{m}_{j}}{{k}^{2}}\right)+\displaystyle \frac{{m}_{j}}{4\pi {a}_{s,j}}\right]{{\rm{\Delta }}}_{j}+\displaystyle \frac{\sqrt{{m}_{1}{m}_{2}}}{4\pi }\gamma {{\rm{\Delta }}}_{3-j}=0.\end{eqnarray} \tag{ 17 }$$

The difference between two aforesaid ensembles consists in the following. Under the first condition, numbers of particles N_j in (16) are fixed. The chemical potentials ${\mu }_{j}(r)$ are determined from these two equations, taking into account (17) which in fact depends on both μ₁ and μ₂. In the 'grand canonical' setting, N₁ and N₂ are not fixed separately. Instead, we apply (17) with the normalization condition for the total number of particles N ≡ N₁ + N₂,

$$\begin{eqnarray}&&2\pi \displaystyle \sum _{j=1,2}{\int }_{0}^{\infty }{n}_{j}(\beta ,\mu (r),{{\rm{\Delta }}}_{j}(r)){r}{\rm{d}}{r}=N,\end{eqnarray} \tag{ 18 }$$

where $\mu (r)$ is the same for both band-components of the condensate. Consequences of this difference, in particular the depletion of the 'weak' band at low temperatures for the 'grand canonical' setting, were considered in [31].

The free energy and the normalization integrals are calculated including both the superfluid core and the normal phase which are spatially separated. It is assumed here that they both rotate with the frequency imposed by a stirring field. A treatment beyond this approximation may lead only to a marginal correction to the obtained results without changing the picture, because we consider vortices near the trap center, i.e., far from the boundary between the superfluid and normal phases. The same assumption has been substantiated and applied in [14, 15].

The term 'fractional vortex' means that winding numbers for vortex states in two bands are not the same. In order to classify different vortex configurations, we label the 'strong' band (with a larger inverse scattering length 1/a_s) by index 1, and the 'weak' band by index 2. Following [23] we use the notation $({n}_{1},{n}_{2},{n}_{c})$ for the classification of vortex states, where n₁ and n₂ are, respectively, the winding numbers in the first and second band-components of the condensate, and n_c is the number of composite vortices, i.e., integer vortices pinned to each other in two bands. Correspondingly, the clouds in two band-components of the condensate classified by the set of indices $({n}_{1},{n}_{2},{n}_{c})$ look as or multivortex configurations corresponding to a cluster of n₁ vortices in the first band and n₂ vortices in the second band, out of which n_c vortices coincide.

In the present treatment, we do not consider superfluids corresponding to two bands rotating relatively to each other, because the rotation frequency is determined by a common stirring field. Therefore the vortex clusters in different bands are in rest with respect to each other, forming a stable configuration. The relative positions between vortices in two bands (both radial and angular) are determined minimizing the free energy with respect to positions of all vortices.

Several resulting vortex clouds are exemplified in figure 1. Analogous clouds for a two-band superconductor can be found, e.g., in [23]. The figure shows the spatial distribution of the amplitude of the pair field $| {{\rm{\Psi }}}_{j}({\bf{r}})|$ in the 'strong' and 'weak' band-components of the condensate for several stable vortex states, both integer and fractional ones. The third column of figure 1 shows the spatial behavior of the total particle density n = n₁ + n₂. The particle density is necessarily modulated in vortices due to the variation of the gap parameters. The modulation of the particle density is experimentally observable. As can be seen from figure 1, fractional vortex configurations can be experimentally resolved from integer vortex states due to a difference between vortex patterns in the two bands. In the present work, we assume that the size of the trap is much larger than the healing lengths of vortices in both band-components. This favors non-composite vortices, except for the case when they are positioned in the center of the trap. Therefore, for the vortex states shown in figure 1, only the integer $(1,1,1)$ state is composite, and the other states are non-composite. It is worth noting however that the distinction between composite and non-composite vortices depends on convention. The interband coupling results in an attraction between vortices in different band-components of the condensate. This attraction tends to minimize a non-compensated phase shift between the pair fields of the two band-components. It leads to the binding of vortices in two bands. When the distance between vortex cores in the two band-components of the condensate is of the same order as the healing lengths (or less), it may become difficult to resolve experimentally. In experiments, such bound vortices can be observed as elongated. In this case, they are non-composite formally, but close to composite vortices from the experimental point of view.

Zoom In Zoom Out Reset image size

3.2. Vortex phase diagrams

3.2.1. Phase diagrams in variables $(\gamma ,\omega )$

Analyzing the areas of stability for vortex configurations, we consider first the evolution of the vortex states as a function of the interband coupling strength γ. When the interband coupling is sufficiently strong, only integer vortex states can survive, because the Josephson coupling contribution to the free energy is positive and proportional to the relative phase difference for the two band-components of the condensate. Therefore we are particularly focused on the range of relatively small γ, where stable fractional vortices can exist.

In order to study this evolution of stable vortex configuration when varying γ, we plot in figures 2 and 3 the vortex phase diagrams in the variable space $({\mathrm{log}}_{10}\gamma ,\omega )$. The logarithmic scale for γ is chosen because of our particular attention to the range of weak interband couplings to track the transitions between fractional and integer vortices. Figure 2 shows the areas of stability for vortex states assuming numbers of particles per unit length in each band fixed, N₁ = N₂ = 500. Figure 3 represents the analogous picture for the grand canonical ensemble of fermions with the total number of particles per unit length N = 10³.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figures 2 and 3, when varying γ at fixed temperatures, a rich variety of stable vortex configurations, both integer and fractional, can be seen in the vortex phase diagrams. The sequence of stable vortex configurations with an increasing interband coupling is physically quite transparent: fractional vortices turn to integer ones at a critical coupling strength which only slightly depends on the rotation frequency. On the contrary, when the 'weak' band is in the BCS regime, the dependence of the critical interband coupling strength on the rotation frequency is rather complicated. As one can see from figures 2 and 3 (especially clearly in panels (e, f), but also in other phase diagrams), each vortex configuration except $(3,3,0)$ may appear twice, at low and high-rotation frequencies. This trend in the sequences of the vortex stability areas for $(\gamma ,\omega )$ phase diagrams is related to the fact that each 'phase boundary' in the $(T,\omega )$ phase diagrams exhibits a bend-over dependence on ω. This bend-over dependence is manifested in figures 2 and 3 through an appearance of two distinct areas for given vortex configurations (except when the upper-frequency areas lie in close vicinity to the maximal rotation frequency ${\omega }_{\max }(T)$ and therefore are not seen). Note that the range of γ for which fractional vortices can exist is rather small, and that the vortex phase diagrams in figures 2 and 3 are plotted in logarithmic scale for γ. This explains the result that the phase boundaries weakly depend on γ in this range.

In the case of sufficiently weak interband coupling, the fractional vortices are energetically more favorable with respect to integer vortices, except maybe the low-temperature and high-rotation-frequency areas, where the configuration $(3,3,0)$ and higher (in our notations, this denotes three integer non-composite vortices) is energetically more favorable with respect to fractional vortex states. This result does not ensure however that the three-vortex integer configuration cannot be suppressed by a fractional configuration of a higher order if we will account for more vortices. When the interband coupling is a bit stronger, both fractional and integer vortex configurations can survive. For a still stronger interband coupling, we can observe only integer vortices and vortex systems. As a rule, stable fractional vortex configurations appear in such a way that vortices survive in the 'strong' band, when the amplitude of the order parameter in the 'weak' band becomes sufficiently small.

3.2.2. Phase diagrams in variables $(T,\omega )$

Next, we analyze the vortex phase diagrams in the variables $(T,\omega )$ which show areas of stability for different vortex states in a two-band rotating Fermi gas with several inverse scattering lengths and different interband coupling strengths. Figure 4 shows the vortex phase diagrams for the ensemble with the fixed numbers of particles per unit length (along the z direction) N₁ = N₂ = 500. In figure 5, the same fermionic system with N ≡ N₁ + N₂ = 1000 is considered for the thermalized state with a common chemical potential. In these figures, panels (a)–(c) show vortex states at a rather small interband coupling strength γ = 10⁻⁴. The panels (d)–(f) represent the vortex phase diagrams for a larger (but still rather small) interband coupling γ = 10⁻². The left-hand panels (a) and (d) of figures 4 and 5 correspond to the case when the 'strong' band is at unitarity, and 'weak' band is in the BCS regime—with, respectively, 1/a_s,1 = 0 and 1/a_s,2 = −0.5. The central panels (b) and (e) show the case of stronger couplings, with, 1/a_s,1 = 0.5 for the 'strong' band and 1/a_s,2 = 0 for the 'weak' band. Finally, the right-hand panels (c) and (f) show vortex phase diagrams for 1/a_s,1 = 1 for the 'strong' band and 1/a_s,2 = 0 for the 'weak' band, i.e., when the 'strong' band is in the BEC regime.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As can be seen from figures 2 and 3, the value γ = 10⁻⁴ represents the specific interesting case, because it lies in a rather narrow range of γ where fractional vortices turn to integer vortices. The other value, γ = 10⁻², is chosen to show phase diagrams where only integer vortex configurations survive. Higher interband coupling strengths do not lead to substantial qualitative changes of phase diagrams, because fractional vortices can be energetically favorable only at weak interband coupling. This explains the choice of parameters made for figures 4 and 5.

In the same way as for a one-band Fermi gas, the vortex phase diagrams are restricted to the area of the superfluid state, i.e., to temperatures T below the critical temperature ${T}_{c}(\omega )$. The critical temperature is frequency-dependent and tends to zero when ω tends to ω₀. The value γ = 10⁻⁴ for a weak interband coupling is chosen here because, on the one hand, it is sufficiently small so that the vortex phase diagrams can contain fractional vortex states, and, on the other hand, this value, although being small, is sufficient to yield rich stability areas for both fractional and integer vortex states. The other set of phase diagrams has been calculated for a relatively small but larger interband coupling strength γ = 10⁻², at which only integer vortex states survive. Similarly to vortex phase diagrams for a one-band Fermi gas [12, 14, 15], the boundaries of the areas of stability for different vortex configurations exhibit a non-monotonic dependence on the rotation frequency and a reentrant behavior as a function of the temperature. The explanation of this bend-over behavior of critical rotation frequencies is the same for two-band and one-band Fermi gases. As discussed in [12, 15], it is related to a decrease of the radius of the superfluid core for a trapped Fermi gas when the rotation frequency becomes close to the confinement frequency. The lower critical rotation frequency corresponds to the ordinary threshold for vortex formation, similarly to that in a superconductor in an external magnetic field, when increasing the field strength. The other (upper) critical rotation frequency appears when the radius of the superfluid core is comparable with the healing length.

When comparing the vortex phase diagrams for ensembles with fixed numbers of particles and with the common chemical potential to each other, one can see that the 'grand canonical' ensemble is more favorable for fractional vortices than that with fixed numbers of particles. This difference is explained by the effect of a partial depletion of a 'weak' band. As found in [31], this depletion is a feedback of the gap parameter to the relative band populations. As a result, the areas of stability for fractional vortices in the case of a common chemical potential (figure 5) are wider than in the case of fixed equal numbers of particles in each band. Moreover, as can be seen comparing figures 4(c) and 5(c), in the BEC regime this depletion can lead to a complete vanishing of integer vortex phases.

Distinct manifestations of the interband coupling in a two-band Fermi gas are related to the temperature T_c,2 equal to the critical temperature of the superfluid phase transition in the 'weak' band-condensate component in the absence of the interband coupling. First, we find a kink of the 'phase boundaries' between different vortex configurations in figures 4 and 5, positioned at T_c,2. This feature is absent in the one-band system considered in [12, 14, 33]. It is directly related to a non-monotonic, resonant peak-shape temperature dependence of the vortex healing length for the 'weak' band ξ₂, analyzed in [31, 41, 42] and termed 'hidden criticality' in [42]. Although there is no true criticality at T ≈ T_c,2, many parameters of state demonstrate non-monotonic behavior at this temperature. The peak position for ξ₂ lies at the critical temperature for the 'weak' band in the absence of the interband coupling, so that it is a fingerprint of the weak-band criticality, which is lost in a two-band system. In figure 6, the evolution of this 'hidden criticality' point is shown as a function of $(T,\omega )$ for the inverse scattering lengths $1/({k}_{F}{a}_{s,1})=0,\ 1/({k}_{F}{a}_{s,2})=-0.5$, the numbers of fermions per unit length N₁ = N₂ = 500, and the interband coupling strength γ = 10⁻². The peak temperature of the healing length ${\xi }_{2}$ diminishes when the rotation frequency rises, with an increasing peak magnitude and a decreasing width.

Zoom In Zoom Out Reset image size

The seemingly unusual sequence of phases in the phase diagrams for weak interband coupling (upper panels of figures 4 and 5) can be transparently understood using simple physical reasoning. First, let us temporarily 'turn off' the interband interaction and consider each band-component independently as a one-band system. For a one-band trapped Fermi gas, this behavior was independently analyzed in [12, 14, 15], using different methods: BdG theory and EFT. These two theories completely agree in predictions on the behavior of vortex states.

The evolution of the number of vortices in a given band-condensate, as a function of the rotation rate at a fixed temperature is non-monotonic. The healing length, which is a parameter characterizing the vortex size, is an increasing function of the rotation frequency. As established in [12, 14, 15], the number of vortices increases when rotation frequency increases as long as the vortex size remains substantially smaller than the size of the superfluid core. When the rotation frequency is further increased, an upper critical rotation frequency appears such that the vortex state turns back to the superfluid state for higher rotation frequencies, because the radius of the superfluid state reduces when the rotation frequency becomes close to the maximal rotation frequency at which superfluidity vanishes. It was also found in [15] that the boundary between the states with different numbers of vortices behaves similarly to the critical rotation frequency for a single vortex. For example, the lower critical rotation frequency for a vortex pair is higher than the lower critical rotation frequency for a single vortex. On the contrary, the upper critical rotation frequency for a vortex pair is lower than the upper critical rotation frequency for a single vortex.

Regarding the variation of the number of vortices as a function of temperature at a finite rotation rate (sufficient for stable vortices at T = 0), the number of vortices monotonically decreases with increasing temperature, because the free energy of the system with vortices becomes higher when the vortex size increases (as an increasing function of temperature). For a two-band system without interband coupling, similarly to the one-band system, the number of vortices in each band-component of the condensate must be a monotonically decreasing function of temperature. Therefore vortices in the 'weak' band in the absence of the interband coupling cannot exist above the critical temperature T_c,2 < T_c,1, while vortices in the 'strong' band can survive above T_c,2 and vanish at higher temperatures close to T_c,1, as can be seen from figures 4 and 5.

The monotonic dependence of vortex numbers as a function of temperature can be broken by interband coupling, because in the weakly-coupled two-band system the healing length for the 'weak' band-component exhibits a pronounced peak near T_c,2, as shown in [31]. In particular, the healing length for the 'weak' band reaches a local maximum at T ≈ T_c,2 and again reduces when temperature slightly exceeds T_c,2. This favors an increasing number of vortices in the 'weak' band-component of the condensate in a relatively narrow temperature area above T_c,2. Since the effects related to 'hidden criticality' are very transparent and necessarily follow from the interband coupling, which is a highly controllable parameter (see, e.g., [21, 22]), the sequence of phases obtained in our work can be experimentally observable and serve as a clear 'smoking gun' evidence for multiband physics in superfluids.

The resonant dependence of the healing length for the 'weak' band influences the free energy of the two-band Fermi gas, as can be seen from figure 7. This figure shows the free energy difference δF for an integer vortex, denoted as $\delta F(1,1)$ for $1/({k}_{F}{a}_{s,1})=0,\ 1/({k}_{F}{a}_{s,2})=-0.5,$ and γ = 10⁻². Thus both the plot for the healing length in figure 6 and the free energy shown in figure 7 correspond to the same values of parameters as the vortex phase diagram in figure 4(d). The area $(T,\omega )$ for $T\lt {T}_{c}(\omega )$ can be subdivided to two areas by the contour indicating δF = 0 shown in the figure explicitly. These areas correspond to δF < 0 and δF > 0 (respectively, inside and outside the contour for δF = 0). As can be seen from figure 7, the behavior of contour lines for the free energy in these two areas is quite different. Also the free energy exhibits a non-monotonic behavior of isoenergetic contours which contain kinks along the path of the peak for the healing length ξ₂. The same kinks are manifested in the contour lines of the vortex phase diagrams.

Zoom In Zoom Out Reset image size

Second, the peak-shape dependence of the healing length in the 'weak' band-component of the condensate ξ₂ on the temperature leads to the fact that ξ₂ decreases in a certain temperature range above T_c,2 to values of the same order as ξ₁. The small healing length ξ₂ facilitates vortex stabilization in the 'weak' band-component of the condensate. As a result, as can be seen from figures 4(a)–(c) and 5(a) and (b), an area of integer vortex states appears above T_c,2 (see, e.g., figure 6 which corresponds to the same set of parameters as figure 4(a)). Note that this area of integer vortex states above T_c,2 is a consequence of a non-zero interband coupling, because at γ = 0 there is no condensate in the 'weak' band. The appearance of this resonant regime is not present in vortex phase diagrams for a one-band Fermi condensate. For comparison, in a one-band system, the sequence of different vortex configurations is such that the areas of stability for higher winding numbers lie completely inside the areas of stability for lower winding numbers [12, 15] in phase diagrams. For a two-band system, this ordering of vortex states can be violated. Namely, at T < T_c,2 we observe the usual sequence of stability areas, where integer vortex configurations change to fractional states with an increasing temperature. For T > T_c,2, the 'anomalous' sequence becomes possible, where fractional states change to integer states when temperature rises. When temperature further increases towards T_c, the sequence of stability areas becomes usual again. The appearance of the 'anomalous' sequence of stability areas for fractional and integer vortex states is arguably the most striking result of the interband coupling for two-band rotating Fermi gases.

3.3. Relevance to experiments

In this work, we have applied the finite-temperature EFT to investigate integer and fractional vortices and multivortex states in rotating two-band Fermi gases throughout the BCS–BEC crossover. As distinct from the one-band system, a rich spectrum of vortex states is realizable in a two-band Fermi gas. Fractional vortices (the states with different winding numbers in the two band-components of the condensate) are stable in this system for sufficiently weak interband couplings. When the interband coupling strength γ exceeds a critical value, which is dependent on the frequency of rotation ω, only integer vortices and vortex clusters can be found. Note that integer vortices in many-vortex clusters are not necessarily composite: the vortex centers in 'strong' and 'weak' band may reside at different distances from the trap center.

The phase boundaries between different stable vortex configurations in the $(T,\omega )$ vortex phase diagrams depend non-monotonically on the rotation frequency, turning to zero both at lower and upper critical rotation frequencies at T = 0. Correspondingly, they exhibit a bend-over temperature dependence, quite similar to those obtained in preceding works for a one-band Fermi gas.

We have characterized the difference between vortex phase diagrams obtained for two-band Fermi condensates in two different regimes: the regime of fixed numbers of particles for each band-components of the condensate, and the grand canonical ensemble when these numbers of particles are determined through the common chemical potential. The depletion of the 'weak' band-component of the condensate in the grand canonical ensemble can lead to a substantial expansion of stability areas with respect to those at fixed numbers of particles. This difference can be experimentally verified.

A striking difference appears between the evolution of vortex configurations as a function of γ in the BCS and BEC regimes. A rich variety of fractional vortex states exist at weak interband couplings in a wide range of rotational frequencies, which turn to integer vortex states at certain γ. The boundaries between different vortex states strongly depend on γ in the BCS regime. On the contrary, in the BEC regime, these boundaries rather weakly depend on the interband coupling strength.

The obtained manifestations of the non-monotonic behavior of the healing length in the 'weak' band-component of the condensate through kinks of phase boundaries and the 'anomalous' resonant sequence of stability areas in vortex phase diagrams can be experimentally accessible. Usually, the 'hidden criticality' phenomena and other effects of the interband coupling cannot be easily resolved experimentally. Therefore, the aforesaid manifestations can represent an effective and transparent 'smoking gun' for the interband coupling in trapped atomic Fermi superfluids.

We thank C A R Sá de Melo and N Verhelst for valuable discussions. This work has been supported by the Research Foundation-Flanders (FWO-Vl), project nrs. G.0115.12N, G.0119.12N, G.0122.12N, G.0429.15N, G.0666.16N, by the Scientific Research Network of the Flemish Research Foundation, WO.033.09N, and by the Research Fund of the University of Antwerp.

The exact effective bosonic action for interacting fermions is obtained introducing auxiliary Hubbard–Stratonovich fields and integrating the partition function out fermion fields as described in [31]. The resulting effective bosonic action for a two-band system reads,

$$\begin{eqnarray}&&{S}_{{\rm{eff}}}^{(2b)}=\displaystyle \sum _{j=1,2}{S}_{{\rm{eff}}}^{(j)}-{\int }_{0}^{\beta }{\rm{d}}\tau \int {\rm{d}}{\bf{r}}\,\displaystyle \frac{\sqrt{{m}_{1}{m}_{2}}\gamma }{4\pi }({\bar{{\rm{\Psi }}}}_{1}{{\rm{\Psi }}}_{2}+{\bar{{\rm{\Psi }}}}_{2}{{\rm{\Psi }}}_{1}),\end{eqnarray} \tag{ A1 }$$

where ${S}_{{\rm{eff}}}^{(j)}$ is the effective field action for a single band-component of the condensate determined by (dropping the band index j):

$$\begin{eqnarray}&&{S}_{{\rm{eff}}}={S}_{B}-\mathrm{Tr}\,\mathrm{ln}[-{{\mathbb{G}}}^{-1}].\end{eqnarray} \tag{ A2 }$$

We follow here the notations of [34], where ${{\mathbb{G}}}^{-1}({\bf{r}},\tau )={{\mathbb{G}}}_{0}^{-1}({\bf{r}},\tau )-{\mathbb{F}}({\bf{r}},\tau )$ is the inverse Nambu tensor. It is subdivided to the free-fermion inverse Nambu tensor ${{\mathbb{G}}}_{0}^{-1}$ and the matrix ${\mathbb{F}}$ proportional to the pair field Ψ:

$$\begin{eqnarray}{{\mathbb{G}}}_{0}^{-1}({\bf{r}},\tau )=\left(\begin{array}{cc}-\tfrac{\partial }{\partial \tau }-{\hat{H}}_{\uparrow } & 0\\ 0 & -\tfrac{\partial }{\partial \tau }+{\hat{H}}_{\downarrow }\end{array}\right),\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}{\mathbb{F}}({\bf{r}},\tau )=\left(\begin{array}{cc}0 & -{\rm{\Psi }}({\bf{r}},\tau )\\ -\bar{{\rm{\Psi }}}({\bf{r}},\tau ) & 0\end{array}\right).\end{eqnarray} \tag{ A4 }$$

Next, the effective action (A2) is expanded as a Taylor series in powers of the pair field:

$$\begin{eqnarray}&&{S}_{{\rm{eff}}}={S}_{B}-\mathrm{Tr}\,\mathrm{ln}[-{{\mathbb{G}}}_{0}^{-1}]+\displaystyle \sum _{p=1}^{\infty }\displaystyle \frac{1}{p}\mathrm{Tr}[{({{\mathbb{G}}}_{0}{\mathbb{F}})}^{p}].\end{eqnarray} \tag{ A5 }$$

In more detail, the trace $\mathrm{Tr}[{({{\mathbb{G}}}_{0}{\mathbb{F}})}^{p}]$ is written as:

$$\begin{eqnarray}\mathrm{Tr}[{({{\mathbb{G}}}_{0}{\mathbb{F}})}^{p}] & = & \displaystyle \int {\rm{d}}{\tau }_{1}\ldots {\rm{d}}{\tau }_{p}\displaystyle \int {\rm{d}}{{\bf{r}}}_{1}\ldots {\rm{d}}{{\bf{r}}}_{p}\mathrm{Tr}[{\mathbb{F}}({{\bf{r}}}_{1},{\tau }_{1}){{\mathbb{G}}}_{0}({{\bf{r}}}_{1}-{{\bf{r}}}_{2},{\tau }_{1}-{\tau }_{2})\\ & & \times {\mathbb{F}}({{\bf{r}}}_{2},{\tau }_{2}){{\mathbb{G}}}_{0}({{\bf{r}}}_{2}-{{\bf{r}}}_{3},{\tau }_{2}-{\tau }_{3})\ldots {\mathbb{F}}({{\bf{r}}}_{p},{\tau }_{p}){{\mathbb{G}}}_{0}({{\bf{r}}}_{p}-{{\bf{r}}}_{1},{\tau }_{p}-{\tau }_{1})].\end{eqnarray} \tag{ A6 }$$

The free-fermion Green's function matrix ${{\mathbb{G}}}_{0}({\bf{r}},\tau )$ is expressed through the Fourier representation,

$$\begin{eqnarray}&&{{\mathbb{G}}}_{0}({\bf{r}},\tau )=\displaystyle \frac{1}{\beta V}\displaystyle \sum _{{\bf{k}},n}{{\mathbb{G}}}_{0}({\bf{k}},n){{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}-{\rm{i}}{\omega }_{n}\tau },\end{eqnarray} \tag{ A7a }$$

$$\begin{eqnarray}{{\mathbb{G}}}_{0}({\bf{k}},n)=\left(\begin{array}{cc}\tfrac{1}{{\rm{i}}({\omega }_{n}-{\rm{i}}\zeta )-{\xi }_{{\bf{k}}}} & 0\\ 0 & \tfrac{1}{{\rm{i}}({\omega }_{n}-{\rm{i}}\zeta )+{\xi }_{{\bf{k}}}}\end{array}\right),\end{eqnarray} \tag{ A7b }$$

with the free-fermion energies

$$\begin{eqnarray}&&{\xi }_{{\bf{k}}}=\displaystyle \frac{{k}^{2}}{2m}-\mu ,\end{eqnarray} \tag{ A8 }$$

and the fermion Matsubara frequencies

$$\begin{eqnarray}&&{\omega }_{n}=\displaystyle \frac{(2n+1)\pi }{\beta }.\end{eqnarray} \tag{ A9 }$$

For spin-imbalanced fermions, μ in (A8) is the averaged chemical potential $\mu =({\mu }_{\uparrow }+{\mu }_{\downarrow })/2$. The chemical potential imbalance parameter is $\zeta =({\mu }_{\uparrow }-{\mu }_{\downarrow })/2$.

The action (A5) is still exact. Further, various approximations are possible. The crudest one is the saddle-point (mean-field) approximation, where ${\rm{\Psi }}({{\bf{r}}}_{n},{\tau }_{n})$ is replaced to a constant value Δ which realizes the least action principle for S_eff. The frequently used approach beyond the saddle-point approximation accounts for Gaussian pair fluctuations about the saddle-point value [34]. It assumes that fluctuations are small with respect to Δ. We apply the alternative method which does not use an assumption of small fluctuations but is related to conditions when the pair fields slowly vary in time and space. Many approximations (in particular, the Ginzburg–Landau and Gross–Pitaevskii theories) use the same assumption. Within this scheme, the first-step extension of the saddle-point approximation for fermions non-uniformly distributed in space (e.g., for fermions trapped to a confinement potential) is the local field approximation (LFA). Within this approximation, we set the space and time variables for each ${\mathbb{F}}({{\bf{r}}}_{n},{\tau }_{n})$ in (A6) to be the same, e.g., ${\mathbb{F}}({{\bf{r}}}_{n},{\tau }_{n})={\mathbb{F}}({{\bf{r}}}_{1},{\tau }_{1})$. This leads to the exact summation over p, resulting in the LFA effective action,

$$\begin{eqnarray*}&&{S}_{{\rm{LFA}}}={\int }_{0}^{\beta }{\rm{d}}\tau \int {\rm{d}}{\bf{r}}\,{{\rm{\Omega }}}_{s}({| {\rm{\Psi }}({\bf{r}},\tau )| }^{2}),\end{eqnarray*}$$

where ${{\rm{\Omega }}}_{s}({| {\rm{\Psi }}({\bf{r}},\tau )| }^{2})$ has the same form as the saddle-point thermodynamic potential but with the coordinate-dependent squared modulus of the order parameter (and also with a coordinate-dependent chemical potential):

$$\begin{eqnarray}{{\rm{\Omega }}}_{s}({| {\rm{\Psi }}| }^{2}) & = & -\displaystyle \int \displaystyle \frac{{\rm{d}}{\bf{k}}}{{(2\pi )}^{3}}\left[\displaystyle \frac{1}{\beta }\mathrm{ln}(2\cosh \beta \zeta +2\cosh \beta {E}_{{\bf{k}}})-{\xi }_{{\bf{k}}}-\displaystyle \frac{m{| {\rm{\Psi }}| }^{2}}{{k}^{2}}\right]\\ & & -\displaystyle \frac{m{| {\rm{\Psi }}| }^{2}}{4\pi {a}_{s}}.\end{eqnarray} \tag{ A10 }$$

Here ${E}_{{\bf{k}}}=\sqrt{{\xi }_{{\bf{k}}}^{2}+{| {\rm{\Psi }}| }^{2}}$ is the Bogoliubov excitation energy.

The next step within this scheme beyond the saddle-point approximation is the gradient expansion of the pair fields ${\rm{\Psi }}({{\bf{r}}}_{n},{\tau }_{n}),\bar{{\rm{\Psi }}}({{\bf{r}}}_{n},{\tau }_{n})$ about one common point, e.g., about $({{\bf{r}}}_{1},{\tau }_{1})$. There is no difference which of the fields is chosen as the background, because the trace (A6) is invariant with respect to cyclic permutations of pair fields. The gradient expansion is then:

$$\begin{eqnarray}{\rm{\Psi }}({{\bf{r}}}_{n},{\tau }_{n}) & = & {\rm{\Psi }}({{\bf{r}}}_{1},{\tau }_{1})+({\tau }_{n}-{\tau }_{1})\displaystyle \frac{\partial {\rm{\Psi }}({{\bf{r}}}_{1},{\tau }_{1})}{\partial {\tau }_{1}}+\displaystyle \frac{1}{2}{({\tau }_{n}-{\tau }_{1})}^{2}\displaystyle \frac{{\partial }^{2}{\rm{\Psi }}({{\bf{r}}}_{1},{\tau }_{1})}{\partial {\tau }_{1}^{2}}\\ & & +({{\bf{r}}}_{n}-{{\bf{r}}}_{1})\cdot {\rm{\nabla }}{\rm{\Psi }}({{\bf{r}}}_{1},{\tau }_{1})+\displaystyle \frac{1}{2}\displaystyle \sum _{i,j=1}^{3}({x}_{n,i}-{x}_{1,i})({x}_{n,j}-{x}_{1,j})\displaystyle \frac{{\partial }^{2}{\rm{\Psi }}({{\bf{r}}}_{1},{\tau }_{1})}{\partial {x}_{1,i}\partial {x}_{1,j}}+...\end{eqnarray} \tag{ A11 }$$

and the same for the conjugated pair field. Here, we restrict this series up to the second order derivatives in time and space. Obviously, the zeroth-order term of the gradient expansion corresponds to the aforesaid LF approximation. Also for other terms of (A11) substituted in (A5), the whole sum over p is collected analytically (for each term of the gradient expansion separately). As a result, we arrive at the EFT action functional derived in [31] (denoting $w\equiv {| {\rm{\Psi }}| }^{2}$):

$$\begin{eqnarray}{S}_{{\rm{EFT}}} & = & {\displaystyle \int }_{0}^{\beta }{\rm{d}}\tau \displaystyle \int {\rm{d}}{\bf{r}}\left[{{\rm{\Omega }}}_{s}(w)+\displaystyle \frac{D(w)}{2}\left(\bar{{\rm{\Psi }}}\displaystyle \frac{\partial {\rm{\Psi }}}{\partial \tau }-\displaystyle \frac{\partial \bar{{\rm{\Psi }}}}{\partial \tau }{\rm{\Psi }}\right)\right.\\ & & +Q(w)\displaystyle \frac{\partial \bar{{\rm{\Psi }}}}{\partial \tau }\displaystyle \frac{\partial {\rm{\Psi }}}{\partial \tau }-\displaystyle \frac{R(w)}{2w}{\left(\displaystyle \frac{\partial w}{\partial \tau }\right)}^{2}\\ & & \left.+\displaystyle \frac{C(w)}{2m}({{\rm{\nabla }}}_{{\bf{r}}}\bar{{\rm{\Psi }}}\cdot {{\rm{\nabla }}}_{{\bf{r}}}{\rm{\Psi }})-\displaystyle \frac{E(w)}{2m}{({{\rm{\nabla }}}_{{\bf{r}}}w)}^{2}\right],\end{eqnarray} \tag{ A12 }$$

with the coefficients:

$$\begin{eqnarray}&&C=\int \displaystyle \frac{{\rm{d}}{\bf{k}}}{{(2\pi )}^{3}}\displaystyle \frac{{k}^{2}}{3m}{f}_{2}(\beta ,{E}_{{\bf{k}}},\zeta ),\end{eqnarray} \tag{ A13 }$$

$$\begin{eqnarray}&&D=\int \displaystyle \frac{{\rm{d}}{\bf{k}}}{{(2\pi )}^{3}}\displaystyle \frac{{\xi }_{{\bf{k}}}}{w}[{f}_{1}(\beta ,{\xi }_{{\bf{k}}},\zeta )-{f}_{1}(\beta ,{E}_{{\bf{k}}},\zeta )],\end{eqnarray} \tag{ A14 }$$

$$\begin{eqnarray}&&E=2\int \displaystyle \frac{{\rm{d}}{\bf{k}}}{{(2\pi )}^{3}}\displaystyle \frac{{k}^{2}}{3m}{\xi }_{{\bf{k}}}^{2}\,{f}_{4}(\beta ,{E}_{{\bf{k}}},\zeta ),\end{eqnarray} \tag{ A15 }$$

$$\begin{eqnarray}&&Q=\displaystyle \frac{1}{2w}\int \displaystyle \frac{{\rm{d}}{\bf{k}}}{{(2\pi )}^{3}}[{f}_{1}(\beta ,{\xi }_{{\bf{k}}},\zeta )-({E}_{{\bf{k}}}^{2}+{\xi }_{{\bf{k}}}^{2}){f}_{2}(\beta ,{E}_{{\bf{k}}},\zeta )],\end{eqnarray} \tag{ A16 }$$

$$\begin{eqnarray}R & = & \displaystyle \int \displaystyle \frac{{\rm{d}}{\bf{k}}}{{(2\pi )}^{3}}\left[\displaystyle \frac{{f}_{1}(\beta ,{E}_{{\bf{k}}},\zeta )+({E}_{{\bf{k}}}^{2}-3{\xi }_{{\bf{k}}}^{2}){f}_{2}(\beta ,{E}_{{\bf{k}}},\zeta )}{3w}\right.\\ & & \left.+\displaystyle \frac{4({\xi }_{{\bf{k}}}^{2}-2{E}_{{\bf{k}}}^{2})}{3}{f}_{3}(\beta ,{E}_{{\bf{k}}},\zeta )+2{E}_{{\bf{k}}}^{2}{{wf}}_{4}(\beta ,{E}_{{\bf{k}}},\zeta )\right].\end{eqnarray} \tag{ A17 }$$

The functions ${f}_{s}(\beta ,\varepsilon ,\zeta )$ are the sums over fermion Matsubara frequencies:

$$\begin{eqnarray}&&{f}_{s}(\beta ,E,\zeta )=\displaystyle \frac{1}{\beta }\displaystyle \sum _{n=-\infty }^{\infty }\displaystyle \frac{1}{{({\omega }_{n}+{\rm{i}}\zeta )}^{2}+{\varepsilon }^{2}}.\end{eqnarray} \tag{ A18 }$$

In [31], they are analytically determined using the recurrence relations:

$$\begin{eqnarray}&&{f}_{1}(\beta ,\varepsilon ,\zeta )=\displaystyle \frac{1}{2\varepsilon }\displaystyle \frac{\sinh (\beta \varepsilon )}{\cosh (\beta \varepsilon )+\cosh (\beta \zeta )},\end{eqnarray} \tag{ A19 }$$

$$\begin{eqnarray}&&{f}_{s+1}(\beta ,\varepsilon ,\zeta )=-\displaystyle \frac{1}{2s\varepsilon }\displaystyle \frac{\partial {f}_{s}(\beta ,\varepsilon ,\zeta )}{\partial \varepsilon }.\end{eqnarray} \tag{ A20 }$$

The range of validity of EFT has been carefully discussed in our preceding work [33]. The incorporation of rotation in the gradient expansion is described in [15]. The gradient expansion for the two-band fermion system is performed for each band-component exactly in the same way as for the one-band system.