Leggett collective excitations in a two-band Fermi superfluid at finite temperatures

The Leggett collective excitations for a two-band Fermi gas with s-wave pairing and Josephson interband coupling in the BCS-BEC crossover at finite temperatures are investigated within the Gaussian pair fluctuation approach. Eigenfrequencies and damping factors for Leggett modes are determined in a nonperturbative way, using the analytic continuation of the fluctuation propagator through a branch cut in the complex frequency plane, as in Phys. Rev. Lett. 122, 093403 (2019). The treatment is performed beyond the low-energy expansion, which is necessary when the collective excitation energy reaches the pair-breaking continuum edge. The results are applied in particular to cold atomic gases at the orbital Feshbach resonance and in a regime far from BEC, which can be relevant for future experiments.


I. INTRODUCTION
The collective excitations predicted by Leggett [1] are specific for multiband neutral or charged superfluids. Physically, they describe oscillations of the relative phase of different superfluid band-components. An existence of this branch of collective excitations necessarily follows from the simple reasoning. In a one-band BCS superfluid or superconductor, there can exist two branches of collective excitations: the phononic (Goldstone) branch of phase oscillations and the pair-breaking (Higgs) branch of amplitude oscillations. In a two-band superfluid with a Josephson interband coupling, there still exists a single gapless Goldstone branch, which corresponds to oscillations of the common phase of the two superfluid band components. For completeness, a branch of phase collective excitations must exist which correspond to the motion of the relative phase of the two components and has a finite gap due to the Josephson coupling. These oscillations, by definition, are absent in a one-band-component superfluid/superconducting system. For cold Fermi gases which are not in the BCS regime, the picture is more complicated, because eigenmodes of collective excitations are not purely phase or amplitude, as obtained in the present work.
Leggett modes were experimentally detected in multiband superconductors [2][3][4]. The two-band superfluidity and Leggett collective modes have been considered theoretically for superconducting MgB 2 [5] and for condensed Fermi gases in the BCS-BEC crossover [6,7]. It was not clear for a long time whether the two-band superfluidity and Leggett collective excitations can be experimentally achieved in cold atomic gases. The idea of the two-band superfluid scenario based on the so-called orbital Feshbach resonance (OFR) was theoretically proposed by R. Zhang et al. [8]. OFR has been later successfully realized in a condensed 173 Yb gas [9,10]. This stimulated theoretical efforts to describe the two-band superfluid states [11,12] and Leggett modes [13,14]. BCS-BEC crossover in a two-band atomic superfluid has been studied at the mean-field level and with fluctuations pointing to unique features of the two-component condensate in recent works [15,16]. The Leggett mode frequencies in cold Fermi gases with OFR were first calculated in Ref. [13] at zero temperature, solving the Gaussian pair fluctuation propagator for undamped modes. The work [14] reprersents an alternative method exploiting the densitydensity response function. In the recent work [17], the massive Leggett mode and gapless phonon mode of the two-band Fermi superfluid are studied for the Sarma phase.
For nonzero temperatures, two-band superfluidity and Leggett modes in cold Fermi gases were investigated earlier in Ref. [18] using effective field theory (EFT), which assumes pair fields slowly varying in time and space. The EFT is applicable only for low-energy excitations, which is not always the case for Leggett modes. Moreover, it is inapplicable for energies in the vicinity of the pair-breaking continuum edge [19]. Therefore the study of Leggett collective excitations at nonzero temperature beyond the low-energy approximation is timely and relevant.
In the present work, we consider Leggett collective modes in a two-band Fermi superfluid using the Gaussian pair fluctuation (GPF) effective action for a twoband system derived in Ref. [18] within the path integral formalism. The frequencies and the damping factors for Leggett modes are determined through complex poles of the fluctuation propagator, similarly to Refs. [19,20], where this method has been applied to pair-breaking and phononic collective excitations in ultracold Fermi gases. In the literature, the damping factor of collective excitations is frequently determined through a second-order perturbation expression with a real eigenfrequency. This is well substantiated only when damping is relatively small with respect to the frequency. Here, the complex poles are found as proposed by P. Nozières [21], so that the real and imaginary part of the complex excitation fre-quency are determined mutually consistently and beyond the aforesaid perturbation approach. In this context, the present calculation can be classified as nonperturbative. The GPF approach is also not restricted to the weak interband coupling regime and is valid in the whole temperature range below the transition temperature T c .
The spectra of Leggett modes are analyzed as a function of the coupling parameters, temperature, and the detuning factor δ, which characterizes the band offset between the two bands. Here, numeric results are presented for the long-wavelength limit, focusing first on a superfluid Fermi gas in the BEC regime, which is relevant for the recent experimental achievements [9,10]. It was shown that the whole BCS-BEC crossover can be realized for superfluid Fermi gases with OFR [13]. Consequently, we also consider Leggett collective modes in the BCS-BEC crossover regime.

II. GPF METHOD FOR TWO BANDS
We use the model fermionic action with s-wave intraband and interband pairing channels proposed in Ref. [18], expressed in Grassmann variables ψ σ,j , ψ σ,j , where σ and j are, respectively, the spin-component and band-component indices. Assuming = 1, this action, with the inverse temperature β, reads: The chemical potentials µ σ,j can be, in general, different for the two band and two spin components. The interaction term U (r, τ ) describes the contact interactions between fermions for both intraband and interband contact interactions: The contributions for the two intraband Cooper pairing channels are with g 1 , g 2 , while those with g 3 , g 4 are related to the interband fermion-fermion scattering. After the Hubbard-Stratonovich transformation with auxiliary bosonic fields as described in [18] we arrive at the effective bosonic action, with the Josephson interband coupling strength γ = (8π/m) (1/g 3 − 1/g 4 ). In should be noted that the Josephson coupling refers to the exchange of Cooper pairs between the two bands. Hence it is a coupling term between two partial condensates of the system, as distinct from the cross-band pairing, i. e., formation of interband Cooper pairs, which is not considered here.
The effective action S (j) ef f for each band depends on the bosonic pair fields Ψ j , Ψ j in the same way as in the one-band theory, e. g., [22]: with the inverse Nambu tensor: Collective excitations are small Gaussian pair fluctuations {φ j , ϕ j } about the saddle-point gap solution ∆ j , ∆ j , The saddle-point gap is determined for a two-band system by the coupled set of gap equations. The gap equations are determined using the standard renormalization of the coupling constants for the contact interaction [22] which removes ultraviolet divergences and results in swave scattering lengths a j : The coupled gap equations are then given by: with the function depending on the Bogoliubov excitation energy E We have denoted here µ j = (µ ↑,j + µ ↓,j ) /2 and ζ j = (µ ↑,j − µ ↓,j ) /2. When setting equal scattering lengths for the two band components, the effective action derived in [18] is reduced to the model Hamiltonian for OFR of Ref. [14].
The GPF action is a quadratic form of fluctuation coordinates in the two-band superfluid, Here, ϕ q,n are 4-dimensional Nambu vectors: where the superscripts indicate the band components. The inverse GPF propagator for the two-band system M 2b (q, iΩ n ) depending on the pair momentum q and the Matsubara frequency Ω n = 2πn/β can be written in the form where κ ≡ mγ/ (4π), I is the unit 2 × 2 matrix, and M (j) (q, iΩ n ) is the inverse GPF propagator for the j-th band component, with the matrix elements [23][24][25] (for simplicity, we drop here the band-component index j): and Our model includes also the particular case of a Fermi superfluid near OFR, choosing equal scattering lengths for the open and closed channels a 1 = a 2 ≡ a and accounting for the band offset through the chemical potentials by µ 2 = µ 1 − δ/2 with the detuning parameter δ. The scattering length a and the interband coupling strength γ are then related to the scattering lengths for singlet and triplet channels (a − , a + ) by: The set of gap equations (8) has two solutions: the inphase solution with ∆ 2 /∆ 1 > 0 and the out-of phase solution with ∆ 2 /∆ 1 < 0 [11,14]. For γ > 0, the in-phase solution is stable, and the out-of-phase one is metastable. It is easy to see from (8) that when the sign of γ is changed (γ → −γ) at the same 1/a j , the in-phase and out-of-phase solutions are swapped keeping the same magnitudes |∆ j | as for γ > 0. The derivation above briefly describes the analytic procedure developed in Ref. [18] which can be found there in detail and which is literally applied in the present work. To summarize, the GPF effective action is the secondorder term of the series for the effective pair field action (3) in powers of fluctuation coordinates. The structure of the GPF action derived in [18] transparently reproduces known one-band GPF actions for each band component coupled through the Josephson coupling. For the description of collective excitations, the GPF approach is in fact equivalent to the Random Phase Approximation (RPA) [26]. It is reliable as long as fluctuation coordinates are sufficiently small, when higher-order terms of the fluctuation expansion can be neglected.

III. DETERMINATION OF EIGENMODES
The spectrum of collective excitations can be determined using the same method as for a one-band Fermi superfluid [19,20]. We solve the equation for the determinant of the inverse fluctuation propagator analytically continued from the set of Matsubara frequencies to the complex frequency plane (iΩ n → z), Formally determined, this equation has no complex roots. They can be found however following the procedure developed by P. Nozières [21]. A function f (z) having a branch cut at the real axis can be analytically continued to the lower semi-panel, Im (z) < 0, using the spectral function determined at the real axis z = ω, which can be analytically continued to the complex z plane, ρ f (ω) → ρ f (z). The analytic continuation for f (z) is then: The complex eigenfrequencies of collective excitations z q ≡ ω q − iΓ q /2, accounting for both the frequency and the damping factor, are obtained as roots of (17) with the analytically continued determinant detM 2b (q, z) related toM 2b (q, z) in accordance with (19).
In general, the spectral function ρ f (ω) contains several non-analytic points which determine several windows for the analytic continuation. Within the interpretation of Refs. [19,20], eigenfrequencies which lie in different windows describe peaks of the pair response function in the corresponding frequency intervals.
The analytic continuation method is complementary to the approach of Ref. [14] where frequencies of collective modes are determined numerically from peak positions of the density-density response function. The complex poles of the GPF propagator give us a semianalytic solution with both eigenfrequencies and damping factors of collective excitations, describing in a consistent way both long-lived and damped modes. Also, the results presented here extend the study of Leggett excitations to nonzero temperatures.
In a two-band superfluid, there exist hybrid phononic branches with sound velocities described by EFT [18], and two pair-breaking branches [19] with frequencies strongly pinned to the two pair-breaking continuum edges 2∆ j . In the present treatment, we focus on Leggett collective modes.
Further on, chemical potentials are assumed to be equal for different spin-components in both bands: µ ↑,j = µ ↓,j = µ j . First, we consider the Fermi superfluid in the BEC regime, which is typical for Fermi gases realized using OFR. The Leggett collective modes for a Fermi gas with OFR are studied at T = 0 in Refs. [13,14]. Here, the treatment is performed for nonzero temperatures. Remarkably, although the starting fermionic model action of our work [18] is different from that of Ref. [14], they lead to the same effective bosonic actions.
We choose units such that the fermion mass m = 1/2 and the total particle density of two band components n = 2/ 3π 2 following Ref. [18]. In these units, the Fermi wave vector of free fermions in two bands for zero detuning is k F ≡ 3π 2 n/2 1/3 = 1, and the corresponding Fermi energy E F = 1. The equation of state is calculated in the present work within the mean-field approximation. It is adequate in the BCS regime, as long as a contribution of fluctuations to the equation of state can be neglected. For stronger couplings, fluctuations beyond the mean-field solution for the pair field lead to a significant reduction of the superfluid transition temperature, and, consequently, to disagreement of calculated quantities with experiment. In order to improve the meanfield results to be more accurate predictions, we use the scaling of the temperature dependence of Leggett mode frequencies expressing them as function of the relative temperature T /T c . Also the equation of state can be written down in terms of T /µ j and µ j /∆ j in addition to (T /E F , k F a j ), as in Refs. [19,20]. Because fluctuations lead to scaling of parameters of state, the dimensionless ratios (T /T c , T /µ j , µ j /∆ j ) appear to be independent on an equation of state, what makes the mean-field approximation for background parameters more predictive, and results in an adequate qualitative description of collective excitations. The scaled picture of the eigenfrequencies can be distorted (like a conformal map) with respect to that obtained a more precise equation of state, but all crossings and anticrossings of different excitation branches must survive. There are of course limitations even for this qualitative description. First, the mean-field approximation is not justified near the superfluid transition temperature, where there is a critical regime driven by strong fluctuations. Thus, for quantitative agreement with experiments in the whole temperature range below T c , taking fluctuations into account is still necessary. Second, anharmonic three-phonon and four-phonon processes, which are beyond GPF, can bring a relatively important contribution to damping at low temperatures. This analysis, as well as the detailed study of other (phononic and pair-breaking) collective excitations in a two-band Fermi superfluid, is beyond the scope of the present work. In Fig. 1, long-wavelength Leggett mode frequencies ω L with the momentum q → 0 are plotted as a function of the detuning factor δ for several values of T /T c using the same values of scattering lengths as in Ref. [14]: 1/k F a + = 1 and a − /a + = 0.8. This leads to the inverse scattering length 1/k F a = 1. 125 and the interband coupling γ = 0.125, which correspond to µ (T c ) /T c | δ=0 ≈ −0.794 and µ/∆| T =0,δ=0 ≈ −0.949. Because γ > 0, the stable background solution for the gaps is the in-phase solution: ∆ 2 /∆ 1 > 0. For T = 0, the obtained Leggett mode frequencies precisely match the results of Ref. [14].
The pair-breaking continuum edge for q = 0 is determined as a minimal value ω b = min (ω b,1 , ω b,2 ) from the pair-breaking edges for the two bands (j = 1, 2): In Fig. 1, these two pair-breaking edges are shown by thin dashed and dotted curves. The true pair-breaking edge ω b corresponds to the lowest of these curves. Physically, ω b indicates the minimal threshold energy above which collective excitations become damped due to decay into fermion pairs. Within the analytic continuation method described above (see also for details Refs. [19,20]), the two pair-breaking edges determine, in general, three windows for the analytic continuation: the The window C does not provide roots for Leggett modes, because it lies completely in the pair-breaking continuum for both band components. The window A generates the undamped solution with the Leggett mode frequency ω The chosen values of the scattering lengths in Fig. 1 are related to the BEC regime, where the chemical potentials for both bands are negative, even at zero detuning. Therefore there are no pair-breaking collective excitations in this regime [19]. As can be seen from the figure, the frequencies of Leggett modes smoothly cross the pair-breaking continuum edge ω b . They acquire only a relatively small damping factor when ω L > ω b but do not vanish. Therefore the visible diminution of the Leggett mode peaks for the response function in the continuum obtained in Ref. [14] can be attributed to a decrease of the spectral weight of Leggett modes rather than to damping.
In Fig. 2, we show the temperature dependence of long-wavelength Leggett mode frequencies and pairbreaking continuum edges for a Fermi superfluid using the same parameters as in Fig. 1, with two values of detuning δ = 0 and δ = 4E F . As distinct from the BCS/unitarity cases, where the Leggett mode frequencies tend to zero at T → T c , in the BEC regime they remain finite when approaching the transition temperature. This difference occurs due to the following reasons. First, previous theoretical studies of Leggett collective excitations were based on a perturbative weak-coupling approach, which is focused on the BCS regime but can fail at strong coupling. In the weak-coupling approach, ω 2 L ∝ ∆ 1 ∆ 2 and therefore the Leggett mode frequency turns to zero at T c . The present treatment does not assume the weak-coupling approximation, and therefore the aforesaid trend is not necessarily fulfilled. Moreover, in the BEC regime, as soon as µ < 0, the pair-breaking continuum edge exceeds the value 2∆ 2 according to (20), and hence ω b,j = 0 at the transition temperature, which favors the survival of Leggett modes. As can be seen from Figs. 1 (b) and 2, the Leggett mode frequencies cross the pair-breaking continuum edge only at a sufficiently large detuning δ c . The value δ c slightly decreases when increasing temperature [see Fig. 1 (b)] but it does not turn to zero. This behavior is specific for the sufficiently far BEC regime, but this not always the case for weaker couplings, where the Leggett mode frequency may reach ω b , as shown below. We can conclude that the BEC regime is promising for detection of Leggett collective excitations. We consider also Leggett modes far from the BEC regime, and with different scattering lengths for open and closed channels, which is more general than the case of OFR where a 1 = a 2 . This setup is not yet reached experimentally, but can represent an interest for future experiments. Fig. 3 shows Leggett mode frequencies at q = 0 as a function of temperature for the interband coupling γ = 0.02 (a) and γ = 0.1 (b). The inverse scattering lengths are taken here 1/a 1 = 0 and 1/a 2 = −0.5. These parameters give us the ratios µ (T c ) /T c | δ=0 ≈ 1.50, µ/∆ 1 | T =0,δ=0 ≈ 0.847, µ/∆ 2 | T =0,δ=0 ≈ 2.11 for the interband coupling γ = 0.02 (panel a) and µ (T c ) /T c | δ=0 ≈ 1.38, µ/∆ 1 | T =0,δ=0 ≈ 0.777, µ/∆ 2 | T =0,δ=0 ≈ 1.59 for γ = 0.1 (panel b). The obtained solutions are compared with the Leggett mode frequency ω (EF T ) L obtained using the low-frequency expansion of the effective action [18] and with the energy of the pair-breaking continuum edge 2∆ 2 .
At low temperatures T ≪ T c , only one root of the dis- persion equation exists, with the frequency denoted here as ω (A) L < 2∆ 2 and shown as a black curve in Fig. 3. For q = 0, the Leggett mode corresponding to this root is undamped, because the quasiparticle-quasiparticle and quasihole-quasihole terms [with iΩ n ± (E k + E k+q )] in the GPF matrix elements (13) and (14) cannot contribute to the damping of collective excitations with energy below 2∆ 2 , while quasiparticle-quasihole terms [with iΩ n ± (E k − E k+q )] do not contribute at all to the matrix elements when q = 0. This picture is quite similar to that of pair-breaking collective modes [19].
As the temperature is increased, the Leggett mode frequency ω (A) L at q = 0 tends to 2∆ 2 , remaining undamped. It does not cross the pair-breaking continuum edge. The underlying physics of this behavior consists in the avoided crossing with the pair-breaking collective excitations, because different eigenmodes have different energies. This anticrossing exists when the chemical potential µ 2 > 0, particularly in the BCS/unitarity regime. It does not exist in the BEC case considered above because of the absence of pair-breaking collective excitations.
The obtained anticrossing of the Leggett mode frequency is drastically different from the solution ω (EF T ) L predicted by the low-energy expansion shown as a green curve in Fig. 3. The latter one does not capture the interplay of the Leggett collective mode with the pairbreaking continuum edge and hence crosses the value ω = 2∆ 2 without any feature.
At sufficiently high temperatures, a second eigenfrequency ω (B) L appears, which is strongly damped in the considered regime far from BEC. The frequency and the damping factor for this root are shown in Fig. 3, respectively, by solid and dashed red curves. In order to interpret it properly, we recall the analytic continuation of the matrix elements through the branch cut. For q = 0, there are 4 non-analytic points on the real axis: L /2 appears when the analytic continuation is performed through the window B, ω 1 < ω < min (ω 2 , ω 3 ). Strictly speaking, this solution is physically relevant (e. g., for the spectral response) only in the window B, where its frequency and damping are shown by thick red curves in Fig. 3. We can formally extend the analytic continuation through the window ω 1 < ω < min (ω 2 , ω 3 ) to the whole lower half-plane, as in Ref. [19]: this solution is shown by thin curves. We can see that the second root appears at lower temperatures starting from a finite ω (B) L with zero damping. As temperature rises, the damping factor of this mode rapidly increases.
We can see that the full solution of the reduced dispersion equation corresponding to Leggett modes can be close to the result of the low-frequency expansion ω (EF T ) L only when the temperature is sufficiently low. The lowfrequency expansion thus becomes inapplicable when the Leggett mode frequency approaches the range close to the pair-breaking continuum edge.
It is a complicated question whether the damped eigenfrequency ω L /2 can be attributed to the pairbreaking collective excitation, because the latter one in a one-band system reveals at small q strong pinning to the pair-breaking continuum edge. Also it hardly can be a purely phase Leggett collective mode, because amplitude fluctuations must bring a substantial contribution to collective excitations in the continuum. In general, this solution is a hybrid of phase and amplitude fluctuations. In order to clarify the physical interpretation of the obtained solutions, it is useful to consider also the spectral weight functions for the fluctuation propagator at a nonzero momentum, since the pair-breaking modes disappear at q = 0.

IV. COLLECTIVE EXCITATIONS FOR A NONZERO MOMENTUM IN SPECTRAL WEIGHT FUNCTIONS
In the literature, phononic and Leggett modes are associated with phase collective excitations, while pairbreaking modes are attributed to amplitude excitations. This is correct only in the far BCS limit, where the phase and amplitude modes can be accurately identified through the poles of the fluctuation propagator. In the general case, amplitude and phase oscillations are coupled, and one cannot extract them explicitly. Thus both amplitude and phase excitations bring contributions to all modes, and we can see fingerprints of all modes in all response functions. The dominating terms remain however the same as in the BCS case.
Consequently, also an alternative method can be used to determine spectra of collective excitations in a superfluid Fermi gas, namely through the spectral weight functions for the GPF propagator [27], which reveals the presence of poles via peaks. The expressions for the spectral weight functions are determined in Appendix A. Here, we discuss relative phase-phase and amplitude-amplitude spectral weight functions, because they are the most representative to extract necessary physical information on collective excitations. The spectral functions (multiplied by q 2 ) are shown as contour plots in Figs. 4 and 5.  behavior of collective excitations in the unitarity/BCS regime, with the inverse scattering lengths 1/k F a 1 = 0, 1/k F a 2 = −0.5 and the interband coupling γ = 0.1. Since the Leggett modes below the pair-breaking continuum edge (i. e., with ω < ω b,2 ) have a very low damping even at nonzero momentum, the spectral weight functions are plotted with the complex argument z = ω + iγ where a relatively small damping parameter γ = 0.03E F is added in order to visualize Leggett modes. The spectra of collective excitations are shown as functions of momentum and temperature. The phase-phase spectral weight functions in Figs. 4 (a, c) exhibit two branches of collective excitations: Leggett and phononic modes. These branches are also visible in Figs. 4 (b, d ), which show amplitude-amplitude spectral weight functions, but with smaller magnitudes. This confirms that these branches are dominated by phase fluctuations. This is especially expressed for phononic excitations, whereas Leggett modes in the unitarity/BCS regime contain a non-negligible part of amplitude fluctuations, as distinct from the far BCS limit, where they are purely phase excitations.
The pair-breaking excitations are not manifested in the phase spectral weights, since they are essentially governed by amplitude fluctuations. They are revealed as finite-width peaks denoted as PB 1 and PB 2 in Figs.  4 (b, d ). The pair-breaking excitation branch for the "weak" (with 1/k F a 2 = −0.5) band PB 2 is clearly visible only for sufficiently low temperatures, here at T 0.5T c . When temperature rises, this pair-breaking excitation peak shifts to higher frequencies away from the pairbreaking continuum edge and gradually broadens. As can be seen from Fig. 4 (d ), this shift shows an avoided crossing with Leggett modes, whose frequencies tend to 2∆ 2 from below. This behavior of pair-breaking and Leggett modes makes clear the interpretation of the finite-width peak above 2∆ 2 in Fig. 3: we can attribute it to damped pair-breaking modes, unpinned from 2∆ 2 due to anticrossing with the Leggett branch. Moreover, due to this anticrossing, pair-breaking modes in a two-band Fermi gas become visible even in the zero-momentum limit, as distinct from the one-band system.
At large momenta, the Leggett branch of collective excitations in Fig. 4 (a) is continued through the pairbreaking continuum edge without anticrossing with the pair-breaking branch, since the latter one dissolves at sufficiently large q. Instead of this, the Leggett excitation branch exhibits anticrossing with the phononic branch when passing the edge. As one can see from the comparison of Figs. 4 (a) and (c), it has more expressed phase character above 2∆ 2 , because the amplitude spectral weight for this branch diminishes. We can conclude therefore that the Leggett branch above the pairbreaking continuum edge acquires a significant contribution of phase fluctuations coming from the phononic branch, and becomes a mix of phononic and Leggett modes.
In Figs 4 (c, d ), we can observe also the upper pair-breaking frequency ω b,1 = 2∆ 1 . In agreement with the above results for the zero-momentum limit, there exist a pair-breaking collective excitation branch close to ω b,1 , which is absent in the phase-phase spectral weight function, Fig. 4 (c), and is clearly seen in the amplitudeamplitude spectrum, Fig. 4 (d ). The spectral weight functions for a two-band system in the BEC regime with 1/k F a 1 = 1, 1/k F a 2 = 0.5 and γ = 0.1 plotted in Fig. 5 show a different behavior of the collective excitations. In the BEC regime the pairbreaking collective excitations do not exist. [19], hence there are no peaks corresponding to these branches in the spectral weight functions. As can be seen from Figs. 5 (c, d ), the Leggett branch in the BEC regime can approach the pair-breaking continuum edge 2∆ 2 only at a temperature rather close to T c . Since pair-breaking excitations are absent in the BEC regime, the Leggett branch can cross the pair-breaking continuum edge, becoming damped,which can also be observed in Figs. 5 (c, d ). For lower temperatures, e. g., T = 0.5T c as shown in Figs.  5 (a, b), both phononic and Leggett modes vary similarly to each other and to 2∆ 2 without crossing. This behavior is simpler than in the unitarity/BCS regimes, but represents an interest, because the BEC regime has been realized in experiments [9,10]. As can be seen from

V. CONCLUSIONS
In our preceding works [19,20], we have determined for a one-band system the frequency and the damping factor for phononic and pair-breaking collective excitations in a nonperturbative way using the analytic continuation for the Gaussian fluctuation propagator. Here, this method has been straightforwardly extended for two-band Fermi superfluids, which reveal Leggett collective excitations absent in one-band systems.
The behavior of the phase collective excitations for interacting two-band Fermi gases was already studied in the limiting cases T = 0 and T → T c [6,7,14]. They mix with each other, resulting in one phononic and one Leggett branch, common for the whole two-band system. The question of existence of the two pair-breaking amplitude branches in a two-band Fermi gas is more subtle. In a one-band Fermi gas with s-wave pairing, the GPF approximation reveals up to two branches of collective excitations: the phononic modes [20], and except the BEC regime, the pair-breaking modes [19]. Intuitively, it could be expected that the number of branches would be conserved when the two band components interact. It was not however a priori clear whether the second pair-breaking mode for a weak band survives. We have found that the GPF approximation for the twoband Fermi gas predicts up to four branches of collective excitations. In the unitarity/BCS regimes, there exist one phononic branch, one Leggett branch and two pairbreaking branches in a two-band Fermi gas. Phononic and Leggett modes consist mainly of phase fluctuations, and pair-breaking modes are related to amplitude fluctuations. The first two branches represent a strong mix of the fluctuations for both band components, while the pair-breaking modes of the two band components are essentially pinned to the two pair-breaking edges. They therefore couple only weakly to each other. However, the pair-breaking collective excitations can strongly interact with Leggett modes, revealing an avoided crossing near the pair-breaking continuum edge. Moreover, the interaction with Leggett excitations reveals the pair-breaking modes even in the long-wavelength limit, where they vanish in a one-band system.
We have found that in the BEC regime for a Fermi superfluid prepared using OFR, the Leggett mode as a function of detuning does not vanish when merging with the pair-breaking continuum, but becomes damped. There is a drastic difference between the behavior of Leggett collective excitations near the pair-breaking continuum edge in the BEC regime and in the crossover regime other than BEC. In the BEC regime, the Leggett mode passes the edge without any feature. Moreover, within the present nonperturbative treatment, Leggett mode frequencies in the BEC regime do not turn to zero at the transition temperature, as distinct from weak-coupling results. On the contrary, far from BEC, the Leggett mode frequency avoids crossing with the edge. This different behavior can be attributed to the interplay of Leggett and pair-breaking collective excitations, which exist in the BCS-BEC crossover sufficiently close to BCS, but disappear in the BEC regime. We can conclude that strong-coupling two-band Fermi superfluids, particularly near the orbital Feshbach resonance, are favorable for the experimental observation of Leggett collective excitations. At present, only the BEC case has been experimentally achieved for two-band atomic gases, but coupling parameters of a Fermi gas close to OFR can be experimentally tuned through the whole BCS-BEC crossover region and even into the BCS regime [28]. This makes the present study relevant for subsequent experiments.
Collective modes of the pair field can be visible through the density response due to the coupling between pair and density fields [14,20,29]. An experimental observation of Leggett collective excitations is a challenging problem, because they, in general, have relatively small spectral weights or high decay rates. Leggett modes have been detected in a two-band superconductor MgB 2 using Raman [2] and tunneling [3] spectroscopy, and, recently, using intense terahertz lights pulses [4]. In the experimentally realized two-band 173 Yb cold atomic gases with OFR, Leggett collective excitations can be hardly resolved, because they are strongly damped in the experimentally achieved stable out-of-phase state [13]. For two-band Fermi gases with both large singlet and triplet scattering lengths, the whole BCS-BEC crossover, which is favorable for long-sought Leggett modes, is promising for future experiments.
The GPF formalism exploited in the present work is based on the model action functional, which is suitable for charge-neutral Fermi gases in the BCS-BEC crossover rather than for superconductors. It is possible to reformulate the method using a BCS-like model, so that it will allow us to calculate in the same way both frequencies and damping rates of collective excitations in multiband superconductors, as MgB 2 and iron-based multiband materials. two-band system The spectral weight functions for the GPF effective action S (2b) GP F are determined using the generating functional P [η, η] which depends on the auxiliary field variables: where, following the notations of Eq.(11), η q,n ,η q,n are 4-dimensional vectors, similarly to ϕ q,n ,φ q,n . This generating functional is calculated analytically using a linear shift of the pair field variables, which results in the expression: The spectral weight functions in the Matsubara representation are then determined by the matrix elements of the GPF propagator, which is the inverse of M 2b : The matrix elements of the GPF propagator are explicitly given by: 12 +M            This matrix gives us, in particular, the intraband and interband spectral weight functions: (1) (2) −1 Analytically continuing the Matsubara frequencies to the complex z plane, we determine the spectral weight functions for the total and relative responses, χ (±) (q, z) = 1 π Im φ (1) q,n ±φ (2) q,n ϕ (1) q,n ± ϕ (2) q,n iΩn→z , (A9) which are expressed through the matrix elements of the GPF propagator using (A5) to (A8): These functions are useful to distinguish between the contributions of the total and relative responses to the analytic solutions of the dispersion equation det M 2b (q, z) = 0. An even more clear identification of different modes is possible using spectral weight functions for total and relative amplitude and phase responses. These spectral weight functions are determined using the inverse GPF propagator in the basis of the amplitude and phase field coordinates. The inverse GPF propagator for a one-band system is determined as in Ref. [23]: Q (q, z) = Q 1,1 (q, z) Q 1,2 (q, z) Q 2,1 (q, z) Q 2,2 (q, z) (A11) with the matrix elements: In this basis, the diagonal matrix elements Q 1,1 and Q 2,2 correspond to amplitude and phase collective excitations, respectively. The non-diagonal matrix elements describe their mixing. Obviously, det Q = det M.