# Instability of superfluid Fermi gases induced by a roton-like density mode in optical lattices

###### Abstract

We study the stability of superfluid Fermi gases in deep optical lattices in the BCS–Bose-Einstein condensation (BEC) crossover at zero temperature. Within the tight-binding attractive Hubbard model, we calculate the spectrum of the low-energy Anderson-Bogoliubov (AB) mode as well as the single-particle excitations in the presence of superfluid flow in order to determine the critical velocities. To obtain the spectrum of the AB mode, we calculate the density response function in the generalized random-phase approximation applying the Green’s function formalism developed by Côté and Griffin to the Hubbard model. We find that the spectrum of the AB mode is separated from the particle-hole continuum having the characteristic rotonlike minimum at short wavelength due to the strong charge-density-wave fluctuations. The energy of the rotonlike minimum decreases with increasing the lattice velocity and it reaches zero at the critical velocity which is smaller than the pair breaking velocity. This indicates that the superfluid state is energetically unstable due to the spontaneous emission of the short-wavelength rotonlike excitations of the AB mode instead due to pair-breaking. We determine the critical velocities as functions of the interaction strength across the BCS-BEC crossover regime.

###### pacs:

## I introduction

The recent realization of superfluidity in Fermi gases Regal (); Bartenstein (); Zwierlein (); Kinast (); Bourdel (); CChin () has opened a new research frontier in ultracold atoms Giorgini (). A great experimental advantage of this system is the ability in controlling atomic interactions using a Feshbach resonance Ketterle (). This allows us to access the crossover between the Bardeen-Cooper-Shrieffer (BCS)-type superfluidity and Bose-Einstein condensation (BEC) of bound molecules, which is referred to as the BCS-BEC crossover Eagles (); Leggett (); Nozieres (); SadeMelo (); Holland (); Timmermans (); Ohashi (); Tamaki (). The study of superfluid Fermi gases in the BCS-BEC crossover is expected to offer new insights into the phenomena of superfluidity and superconductivity, which can be applied in various fields such as condensed-matter physics, nuclear physics, and particle physics.

One of the most dramatic features of a superfluid system is the dissipationless superfluid flow Pitaevskii1 (). In particular, critical velocities of superfluid flow have attracted much interest in various systems such as superfluid He Tilley (), superfluid He Vollhardt (), and atomic Bose-Einstein condensates Raman (); Onofrio (). It is well known that the underlying mechanisms for the instability of dissipationless flow are different in the BCS and BEC regions in a uniform system. Namely, the instability of BCS-type superfluids is considered to be dominated by Cooper pair breaking Tinkham (), whereas the instability of Bose superfluids is induced by spontaneous emission of phonon excitations Landau (). It is of interest to study how the mechanism of the instability in superfluid Fermi gases changes in the BCS-BEC crossover.

Recently, Miller et al. investigated experimentally the stability of superfluid flow in Fermi gases in shallow one-dimensional (1D) optical lattices across the BCS-BEC crossover Miller (). They measured superfluid critical velocities, at which the number of condensed atoms starts to decrease, by moving the optical lattice potential through the atomic cloud for different values of interatomic interaction and lattice depth Miller (). The measured critical velocities showed a crossover behavior between the BCS and BEC regimes taking a maximum value at the crossover regime Miller (). Critical velocities in superfluid Fermi gases in the BCS-BEC crossover have been also addressed theoretically in several papers Combescot (); Spuntarelli (); Pitaevskii2 (). The observed crossover behavior of the critical velocities has been predicted in Refs. Combescot (); Spuntarelli (). However, most of the theoretical papers are limited within a uniform system Combescot () or a system in the presence of a single potential barrier Spuntarelli (), which cannot be directly compared to the experiment using optical lattices in Ref. Miller (). In Ref. Pitaevskii2 (), sound propagation in superfluid Fermi gases in optical lattices has been studied using the hydrodynamic approximation. However, microscopic calculation of the critical velocities of superfluid Fermi gases in optical lattices has not been worked out yet.

In this paper, we study the stability and critical velocities of superfluid Fermi gases in deep one-dimensional, two-dimensional (2D), and three-dimensional (3D) optical lattices in the BCS-BEC crossover at zero temperature. We apply the generalized random-phase approximation (GRPA) developed by Côté and Griffin Cote () to the attractive tight-binding Hubbard model in order to calculate the excitation spectra in the presence of a moving optical lattice. For the stability of Fermi gases in the BCS-BEC crossover, two kinds of excitations play crucial roles. One is the single-particle excitation which arises when Cooper pairs are broken. The other is the collective density-fluctuation mode, the so-called Anderson-Bogoliubov (AB) mode Anderson1 (); Bogoliubov (). In a uniform system, the single-particle excitation induces the instability of superfluid flow in the BCS regime, while phonon excitation of the Bogoliubov mode, which corresponds to the AB mode in the BCS regime, induces the instability in the BEC regime Combescot (); Spuntarelli (); Miller (). We find that in deep 1D, 2D, and 3D optical lattices, the excitation spectrum of the AB mode has a characteristic rotonlike structure and lies below the particle-hole continuum due to the strong charge-density-wave (CDW) fluctuation. The energy of the rotonlike minimum decreases with increasing the superfluid velocity and it reaches zero before the particle-hole continuum does, i.e., before pair breaking occurs. As a result, in contrast to the uniform case, the instability of superfluid flow in 1D, 2D, and 3D optical lattices is induced by the rotonlike excitations of the AB mode rather than by pair-breakings. We calculate the critical velocities at which spontaneous emission of the rotonlike excitations occurs as functions of the interaction strength in the entire BCS-BEC crossover regime.

This paper is organized as follows. In Sec. II, we present the model and formalism. We introduce the tight-binding Hubbard model and the Green’s function formalism for the GRPA. In Sec. III, we present the results for the stability and the critical velocities of superfluid Fermi gases in 1D, 2D, and 3D optical lattices. We calculate the excitation spectra and determine critical velocities as functions of the attractive interaction. We summarize our results in Sec. IV.

## Ii Model and formalism

In this section, we summarize the Green’s function formalism applied to an attractive Hubbard model. This is necessary for the calculation of response functions in GRPA. The excitation spectra of collective modes can be obtained as the poles of the response functions. Since our major interest is in the AB mode, we calculate the density response function assuming an external field coupled with density. To discuss the stability of superfluid states, we extend the previous work for the ground state Belkhir (); Koponen () to the current-carrying states.

### ii.1 Green’s function formalism

We consider two-component atomic superfluid Fermi gases with equal populations loaded into optical lattices. We suppose that the optical lattice potential is moving with a constant velocity in the laboratory frame. If the velocity of the lattice potential does not exceed the critical velocity, the Fermi gas remains stable in the laboratory frame due to its superfluidity. This situation can be described equivalently in the frame fixed with respect to the lattice potential as a superfluid Fermi gas flowing with a constant quasimomentum , where is the mass of a fermion. In the following, we describe the system in the frame fixed with respect to the lattice potential. Namely, we assume a time-independent lattice potential and a supercurrent with the quasimomentum .

We assume that the optical lattice potential is sufficiently deep so that the tight-binding approximation is valid. Thus, the system can be described by a single-band Hubbard model as (we set )

(1) |

where is the annihilation operator of a fermion on the th site with pseudospin . Here, is the nearest-neighbor hopping energy, is the on-site interaction energy, and is the chemical potential. We assume an attractive interaction between atoms ().

In order to calculate the density response function, we introduce a fictitious time-dependent external field which is coupled with the density. The Hamiltonian with the external field is given by

(2) | |||||

(3) |

where is the number operator. The density response function is obtained by taking a functional derivative of the single-particle Green’s function by the external field. This will be carried out in Sec. II.3.

We use the imaginary time Green’s function technique KadanoffBaym (). The Heisenberg representations of the annihilation and creation operators in the imaginary time are defined as,

(4) | |||||

(5) |

We introduce the normal and anomalous single-particle Green’s functions, respectively, as Gorkov ()

(6) | |||||

(7) |

where represents the time-ordering operator with respect to . Using the Nambu representation Nambu () with

(8) |

the single-particle Green’s function can be written in the matrix form as

(9) | |||||

We note that in the absence of the external field, the Green’s function at equal sites and imaginary times is given by

(10) | |||||

(11) |

where and is the pair annihilation operator which is related to the wave function of Cooper pairs, as we discuss below.

From the equations of motion for and , we obtain the equation for the matrix Green’s function in Eq. (9), as

(12) |

where is the Pauli matrix

(13) |

The non-interacting Green’s function satisfies

(14) |

From Eqs. (12) and (14), the Green’s function satisfies the Dyson equation

(15) | |||||

where and is the temperature. In Eq. (15), the self-energy is given by

(16) |

Here, we introduced the inverse matrix Green’s function , which satisfies

(17) |

Using Eq. (17), Eq. (15) can be simplified as

(18) |

To calculate the self-energy in Eq. (16), we use the Hartree-Fock-Gor’kov (HFG) approximation Cote (); Gorkov ()

(19) |

Thus, the self-energy in the HFG approximation is given by

(20) | |||||

In the presence of supercurrent with velocity , Cooper pairs are Bose-condensed into the state with the center-of-mass quasimomentum . Since the anomalous Green’s function can be regarded as the wave function of Cooper pairs Gorkov (), it can be written as

(21) |

where is the superfluid gap and is the location of the th site. The exponential factor on the right-hand side of Eq. (21) describes the supercurrent with quasimomentum .

In the presence of the supercurrent, the normal and anomalous Green’s functions can be written as

(22) | |||||

(23) |

respectively, where is the relative coordinate and is the center-of-mass coordinate of the Cooper pair. Here, and are functions of .

To eliminate the phase factors associated with the supercurrent, it is convenient to introduce an operator and a matrix Green’s function as,

(24) | |||||

(25) | |||||

where the matrix for the unitary transformation between and is given by

(26) |

Using Eq. (11) and (25), in the absence of the external field, reduces to

(27) |

In the following sections, we derive equations for .

### ii.2 Equilibrium Green’s function

In this section, we calculate the equilibrium Green’s function in the absence of the external field within the HFG approximation introduced in Eq. (19). Since in Eq. (25) is a function of , we define the Fourier transform of the Green’s function as

(28) |

where is the number of lattice sites and is the Fermi Matsubara frequency. We note that cannot be expanded as Eq. (28) because the phase factor in Eq. (23) describing the supercurrent depends on the center-of-mass coordinate .

From Eq. (15), we obtain the Dyson equation in Fourier space in the absence of the external field as

(29) |

Here, is the Fourier transform of . In Eq. (29), the unperturbed Green’s function is given by

(30) |

where is the kinetic energy, is the index for spatial dimension, and is the lattice constant. From Eq. (19), we obtain the self-energy in Eq. (29) as

(31) |

In deriving Eq. (31), we shifted the chemical potential by the Hartree-Fock energy , where is the average number of atoms per site.

By solving Eq. (29) with the self-energy in Eq. (31), we obtain the single-particle Green’s function as

(32) |

where

(33) | |||||

(34) | |||||

(35) | |||||

(36) | |||||

(37) |

Here, the single-particle excitation energy is given by

(38) |

where , , and . Equation (38) explicitly shows that depend on the superfluid velocity . The single-particle excitation spectrum in Eq. (38) is shown in Fig. 1 for different superfluid velocities. In Fig. 1, the energy gap becomes smaller as increases. When the energy gap reaches zero, pair breaking occurs, i.e., , where is the pair-breaking velocity Rodriguez ().

We determine the superfluid gap and the chemical potential by solving self-consistently the number equation

(39) |

and the gap equation

(40) |

Equations (39) and (40) are obtained from the diagonal and off-diagonal elements of the Green’s function in Eq. (32). Here, is the Fermi distribution function. This scheme of solving Eqs. (39) and (40) self-consistently interpolates the weak-coupling BCS limit and strong-coupling BEC limit at low temperature when the fluctuation effect due to pairs with finite center of mass momenta can be neglected Leggett (). Throughout the work including the calculation of and , we assume .

In Eq. (40), depends on superfluid velocity via , , and . To explicitly show this, we plot in 3D at as a function of in Fig. 2. Note that when at , the superfluid gap does not vanish (), but the superfluid state is destabilized due to pair breaking.

### ii.3 Response function

In this section, we calculate the density response function. The density response function can be derived by taking a functional derivative of the density by the external field as

(41) |

Here, we introduce the three-point correlation function as

(42) |

where . Using Eq. (42), the density response function can be written as

(43) | |||||

where is the density fluctuation operator and

(44) |

In deriving Eq. (43), we used the functional differentiation of the single-particle Green’s function by the external field as

(45) |

We note that Eq. (45) is valid within the linear-response regime.

We derive the equation of motion for the three-point correlation function . Differentiating Eq. (17) with respect to the external field, we obtain

(46) | |||||

Using the HFG approximation in Eq. (19), the three-point correlation function satisfies

(47) | |||||

where the lowest-order correlation function is given by . Thus, the density response function can be obtained by solving Eq. (47) which is referred to as the GRPA equation Cote ().

To see the diagrammatic structure of Eq. (47) more clearly, it is useful to rewrite Eq. (47) in terms of the irreducible correlation function

(48) |

Using Eq. (48), Eq. (47) reduces to

(49) |

It is clear from Eq. (48) that includes the ladder diagrams. On the other hand, Eq. (49) includes the bubble diagrams which lead to the random-phase approximation (RPA) Cote (). In a homogeneous system, the contribution from the bubble diagrams for an attractive interaction can be neglected Cote (). However, in our lattice system, since the bubble diagrams induce the instability due to the CDW fluctuation, it is crucial for the analysis of the stability of the system to keep the bubble diagrams in Eq. (49). We compare the excitation spectra with and without the contribution from the bubble diagrams and give a detailed discussion of the effects of the CDW fluctuations on the stability of the system in Sec. III.

We solve Eqs. (48) and (49) to calculate the density response function in Eq. (41). We define the Fourier transform of as

(50) |

where is the Bose Matsubara frequency. From Eq. (43), the Fourier component of the density response function is given by .

The Fourier transforms of Eqs. (48) and (49) are represented as

(51) |

and

(52) |

respectively, where

(53) |

In order to rewrite Eqs. (51) and (52) in simpler forms, we define a column vector as

(54) |

Here, we have used the notation . The same notation is adapted for and . In addition, we define a matrix as

(55) |

where

(56) |

From Eq. (56), describes a single bubble diagram of a particle-hole excitation. Thus, Eqs. (51) and (52) can be written in the matrix forms as

(57) | |||||

(58) |

respectively. Solving Eq. (57), we obtain

(59) |

Here, is the unit matrix. After the analytic continuation (we take the limit after the calculation), we obtain the density response function as

(60) |

The excitation spectrum of the AB mode is obtained from the pole of Eq. (60).