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

We report the equilibrium vortex phase diagram of a rotating two-band Fermi gas confined to a cylindrically symmetric parabolic trapping potential, using the recently developed finite-temperature effective field theory [Phys. Rev. A $\bf{94}$, 023620 (2016)]. A non-monotonic resonant dependence of the free energy as a function of the temperature and the rotation frequency is revealed for a two-band superfluid. We particularly focus on novel features that appear as a result of interband interactions and can be experimentally resolved. The resonant dependence of the free energy is directly manifested in vortex phase diagrams, where areas of stability for both integer and fractional vortex states are found. The study embraces the BCS-BEC crossover regime and the entire temperature range below the critical temperature $T_{c}$. Significantly different behavior of vortex matter as a function of the interband coupling is revealed in the BCS and BEC regimes.


I. INTRODUCTION
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][3][4][5] and theoretical [6][7][8][9][10][11][12][13][14][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 2 , 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][17][18]. Recently, a two-band Fermi superfluid has been successfully created in an ultracold gas of 173 Yb atoms [19,20] using the orbital Feshbach resonance predicted in Ref. [21] (see also Ref. [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 [25][26][27]. 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 [23]. Moreover, vortex matter in multiband superconductors in the BCS-BEC crossover regime has also attracted much attention [24]. 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 Refs. [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 Bogoliubov -de Gennes equations for the same problem have been solved only numerically.
Besides our studies, there are several modifications of the effective field theory 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][36][37]).
As analyzed in Ref. [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 Ref. [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 Sec. II, the used method and approximations are described. In Sec. III, the parameters of state and the vortex phase diagrams are analyzed in detail. Our results are summarized in Sec. IV.

II. METHOD
We apply the finite-temperature EFT developed in Ref. [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 ω 0 . This setup can be relevant for ultracold fermionic atoms in elongated traps, as, e. g., in Refs. [2][3][4][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 Ref. [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 Ref. [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 Ref. [15]. Consequently, in this section we reproduce the formalism only briefly, because detailed derivations, discussions and proofs can be found in Refs. [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., Refs. [34,35]. Throughout the treatment, we use units such that = 1, the fermion mass for a "strong" band m 1 = 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 ≈ (3ω x ω y ω 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 ≡ √ 2mE F / .
The thermodynamic quantities are calculated on the basis of the partition function, where β = 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 j-th band, with the one-particle Hamiltonian in the rotating frame of reference, which allows for independent populations in different bands and different spins σ, given by where A j (r) is the rotational vector potential A j (r) = m j [ω × r], and ω ≡ ωe z is the rotation vector. The fermion-fermion interaction U (r, τ ) for a two-band system assumes both intraband and interband contact interactions, with intraband coupling constants g jj (j = 1, 2) and interband coupling constants g (a) 12 and g (p) 12 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 Ref. [31]. Here, details of the derivation of the effective field action are given in Appendix A. After the Hubbard-Stratonovich transformation using bosonic pair fields Ψ 1 , Ψ 2 , 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 (2b) EF T takes the form: where γ is the strength of the interband coupling expressed through the interband scattering lengths, (j) EF T for the j-th band is determined by: The coordinate-dependent thermodynamic potential Ω s,j |Ψ j | 2 and the coefficients of the gradient expansion given in Refs. [15,31] and explained in detail in Appendix A. The rotation is incorporated in the effective action as described in Ref. [15]. It leads to the appearance of the linear gradient term iD j A j · Ψ j ∇ r Ψ j − Ψ j ∇ rΨj in the effective bosonic action, which vanishes in the absence of rotation. As shown in Ref. [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].

A. 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), where F (j) is the contribution to the free energy from the j-th band-component of the condensate: The pair field Ψ j (r) for a particular vortex state with n vortices in the j-th bandcomponent of the condensate is described here by the variational function, which is a product of the background pair field ∆ j (r) and the vortex factor: where a ν,j (r) and θ ν,j (r) are, respectively, the relative amplitude and the phase for the ν-th vortex in the j-th band component of the condensate. The phase for a vortex θ ν,j (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 Ref. [15], are modeled by variational functions, 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 Ref. [15], by considering the difference between the free energies with and without vortices: 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 Ψ j (r) = ∆ j (r), and choosing a vortex state with the lowest free energy.
The background distribution of the pair field ∆ 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, where n (β, µ (r) , ∆ (r)) is the fermion density depending on the radial variable r = x 2 + y 2 through the background pair field ∆ (r) (we neglect here the feedback of vortices to the density normalization) and the chemical potential, according to (4), The coordinate-dependent gap parameter is determined through the local gap equation [31], consistently with the normalization condition (14).
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 Ref. [42]. 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,42].
For a two-band system, coordinate-dependent parameters of state are determined from two normalization conditions, together with the coupled set of two gap equations for the two-band Fermi gas given by Eq.
(28) of Ref. [31], 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 µ j (r) are determined from these two equations, taking into account (17) which in fact depends on both µ 1 and µ 2 . In the "grand canonical" setting, N 1 and N 2 are not fixed separately.
Instead, we apply (17) with the normalization condition for the total number of particles where µ (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 Ref. [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 Refs. [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 Ref. [25] we use the notation (n 1 , n 2 , n c ) for the classification of vortex states, where n 1 and n 2 are, respectively, the winding numbers in the first and second bandcomponents 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 bandcomponents 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 1 vortices in the first band and n 2 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 Fig. 1. Analogous clouds for a two-band superconductor can be found, e. g., in Ref. [25]. The figure shows the spatial distribution of the amplitude of the pair field |Ψ j (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 Fig. 1 shows the spatial behavior of the total particle density n = n 1 + n 2 .
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 Fig. 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 In order to study this evolution of stable vortex configuration when varying γ, we plot in band. Finally, the right-hand panels (c,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.
As can be seen from Figs 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 (ω). The critical temperature is frequency-dependent and tends to zero when ω tends to ω 0 . The value γ = 10 −4 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 −2 , 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 bendover behavior of critical rotation frequencies is the same for two-band and one-band Fermi gases.
As discussed in Refs. [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 Ref. [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 (Fig. 5) are wider than in the case of fixed equal numbers of particles in each band. Moreover, as can be seen comparing Fig. 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 Figs. 4 and 5, positioned at T c,2 . This feature is absent in the one-band system considered in Refs. [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 ξ 2 , analyzed in Refs. [31,40,41] and termed "hidden criticality" in Ref. [41]. 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 ξ 2 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 Fig. 6, the evolution of this "hidden criticality" point is shown as a function of (T, ω) 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 1 = N 2 = 500, and the interband coupling strength γ = 10 −2 . The peak temperature of the healing length ξ 2 diminishes when the rotation frequency rises, with an increasing peak magnitude and a decreasing width.
The seemingly unusual sequence of phases in the phase diagrams for weak interband coupling (upper panels of Figs. 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 Refs. [12,14,15], using different methods: Bogoliubov -de Gennes 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 Refs. [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 Ref. [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 Figs. 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 Fig. 7. This figure shows the free energy difference δF for an integer vortex, denoted as δF (1, 1) for 1/ (k F a s,1 ) = 0, 1/ (k F a s,2 ) = −0.5, and γ = 10 −2 . Thus both the plot for the healing length in Fig. 6 and the free energy shown in Fig. 7 correspond to the same values of parameters as the vortex phase diagram in Fig. 4 (d ). The area (T, ω) for T < T c (ω) 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 Fig.   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 ξ 2 . The same kinks are manifested in the contour lines of the vortex phase diagrams.
Second, the peak-shape dependence of the healing length in the "weak" band-component of the condensate ξ 2 on the temperature leads to the fact that ξ 2 decreases in a certain temperature range above T c,2 to values of the same order as ξ 1 . The small healing length ξ 2 facilitates vortex stabilization in the "weak" band-component of the condensate. As a result, as can be seen from Figs. 4 (a, b, c) and 5 (a, b), an area of integer vortex states appears above T c,2 [see, e. g., Fig. 6 which corresponds to the same set of parameters as Fig. 4 (a)]. Note that this area of integer vortex states above T c,2 is a consequence of a nonzero 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.
Recently, a new kind of orbital Feshbach resonance has been experimentally achieved in ultracold 173 Yb fermionic atoms [19,20]. As shown in Refs. [43][44][45], the many-body Hamiltonian describing this system with an orbital Feshbach resonance is exactly analogous to that of a two-band s-wave fermionic superfluid with contact interaction. We can map the parameters corresponding to these experiments on the present theory. According to Refs. [19,20,43], the orbital-singlet and orbital-triplet scattering lengths are, respectively, a s+ ≈ 1900a 0 and a s− ≈ 200a 0 , where a 0 is the Bohr radius. In the experimental conditions of Ref. [19], an ultracold gas of N = 6 × 10 4 fermionic 173 Yb atoms is confined to a cigar-like trap with ω x = 2π × 13 Hz, ω y = 2π × 188 Hz, and ω z = 2π × 138 Hz. This gives us the Fermi energy E F = (3ω x ω y ω z N) 1/3 ≈ 2.604 × 10 −30 J, and the parameter k F ≈ 1.1597 × 10 7 m −1 .
Using the set of units described in Sec. II and following to Ref. [43], we find the dimensionless parameters corresponding to the experiment [19] to be 1/a s,1 = 1/a s,2 ≈ 4.50, and γ ≈ 3.64.
This strong interband coupling only allows for integer vortices.
The scattering lengths are highly controllable for ultracold Fermi gases. One can therefore expect that low values of the interband coupling strength, which are needed for fractional vortices in multiband superfluids, are realizable in future experiments.

IV. CONCLUSIONS
In this work, we have applied the finite temperature effective field theory 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, ω) 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 twoband 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. 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 Ref. [31]. The resulting effective bosonic action for a two-band system reads, where S (j) ef f is the effective field action for a single band-component of the condensate determined by (dropping the band index j): We follow here the notations of Ref. [34], where G −1 (r, τ ) = G −1 0 (r, τ ) − F (r, τ ) is the inverse Nambu tensor. It is subdivided to the free-fermion inverse Nambu tensor G −1 0 and the matrix F proportional to the pair field Ψ: Next, the effective action (A2) is expanded as a Taylor series in powers of the pair field: In more detail, the trace Tr [(G 0 F) p ] is written as: The free-fermion Green's function matrix G 0 (r, τ ) is expressed through the Fourier representation, with the free-fermion energies and the fermion Matsubara frequencies For spin-imbalanced fermions, µ in (A8) is the averaged chemical potential µ = (µ ↑ + µ ↓ ) /2.
The action (A5) is still exact. Further, various approximations are possible. The crudest one is the saddle-point (mean-field) approximation, where Ψ (r n , τ n ) is replaced to a constant value ∆ which realizes the least action principle for S ef f . 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 F (r n , τ n ) in (A6) to be the same, e. g., F (r n , τ n ) = F (r 1 , τ 1 ). This leads to the exact summation over p, resulting in the LFA effective action, where Ω s |Ψ (r, τ )| 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): Here E k = ξ 2 k + |Ψ| 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 Ψ (r n , τ n ) ,Ψ (r n , τ n ) about one common point, e. g., about (r 1 , τ 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: Ψ (r n , τ n ) = Ψ (r 1 , τ 1 ) + (τ n − τ 1 ) ∂Ψ (r 1 , τ 1 ) ∂τ 1 + 1 2 (τ n − τ 1 ) 2 ∂ 2 Ψ (r 1 , τ 1 ) ∂τ 2 1 + (r n − r 1 ) · ∇Ψ (r 1 , τ 1 ) + 1 2 3 i,j=1 (x n,i − x 1,i )(x n,j − x 1,j ) ∂ 2 Ψ (r 1 , τ 1 ) ∂x 1,i ∂x 1,j + ... (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 Ref. [31] (denoting w ≡ |Ψ| 2 ):