#### Abstract

A low-energy theory of the Nambu-Goldstone excitation spectrum and the corresponding speed of sound of an interacting Fermi mixture of Lithium-6 and Potassium-40 atoms in a two-dimensional optical lattice at finite temperatures with the Fulde-Ferrell order parameter has been formulated. It is assumed that the two-species interacting Fermi gas is described by the one-band Hubbard Hamiltonian with an attractive on-site interaction. The discussion is restricted to the BCS side of the Feshbach resonance where the Fermi atoms exhibit superfluidity. The quartic on-site interaction is decoupled via a Hubbard-Stratonovich transformation by introducing a four-component boson field which mediates the Hubbard interaction. A functional integral technique and a Legendre transform are used to give a systematic derivation of the Schwinger-Dyson equations for the generalized single-particle Green’s function and the Bethe-Salpeter equation for the two-particle Green’s function and the associated collective modes. The numerical solution of the Bethe-Salpeter equation in the generalized random phase approximation shows that there exist two distinct sound velocities in the long-wavelength limit. In addition to low-energy (Goldstone) mode, the two-species Fermi gas has a superfluid phase revealed by two roton-like minima in the asymmetric collective-mode energy.

#### 1. Introduction

Optical lattices are formed by the interference of counter propagating laser beams. If the laser beams have equal frequencies, the gases of ultracold alkali atoms can be trapped in periodic potentials (microtraps) created by standing waves of laser light. Because of the Stark effect the ground-state alkali atoms couple to the electromagnetic field via an induced electric dipole moment. From theoretical point of view, the simplest approach to the trapped fermions is the tight-binding approximation, which requires sufficiently deep lattice potential. In the tight-binding limit, two alkali atoms of opposite pseudospins on the same site have an interaction energy , while the probability to tunnel to a neighboring site is given by the hopping parameters. The hopping parameters as well as the interaction energy depend on the depth of the lattice potential and can be tuned by varying the intensity of the laser beams. We assume that the interacting fermions are in a sufficiently deep periodic lattice potential described by the Hubbard Hamiltonian. We restrict the discussion to the case of atoms confined to the lowest-energy band (single-band Hubbard model), with two possible states described by pseudospins . We consider different amounts of and atoms in each state (, ) achieved by considering different chemical potentials and . There are atoms distributed along sites, and the corresponding filling factors are smaller than unity. The Hubbard Hamiltonian is defined as follows: where is the single electron hopping integral and is the density operator on site . The Fermi operator () creates (destroys) a fermion on the lattice site with pseudospin projection . The symbol means sum over nearest-neighbor sites of the two-dimensional lattice. The first term in (1) is the usual kinetic energy term in a tight-binding approximation. All numerical calculations will be performed assuming that the hopping (tunneling) ratio . In our notation the strength of the on-site interaction is positive, but the negative sign in front of the interaction corresponds to the Hubbard model with an attractive interaction. In the presence of an (effective) attractive interaction between the fermions, no matter how weak it is, the alkali atoms form bound pairs, also called the Cooper pairs. As a result, the system becomes unstable against the formation of a new many-body superfluid ground state. The superfluid ground state comes from the symmetry breaking and it is characterized by a nonzero order parameter, which in the population-balanced case is assumed to be a constant in space . Physically, it describes superfluid state of Cooper pairs with zero momentum. Superfluid state of Cooper pairs with nonzero momentum occurs in population-imbalanced case between a fermion with momentum and spin and a fermion with momentum and spin . As a result, the pair momentum is . A finite pairing momentum implies a position-dependent phase of the order parameter, which in the Fulde-Ferrell [1] (FF) case varies as a single plane wave , where is a real quantity. The order parameter also can be a combination of two plane waves as in the case of the Larkin-Ovchinnikov [2, 3] (LO) superfluid states. In both cases we are dealing with a spontaneous translational symmetry breaking and with an inhomogeneous superfluid state. When continuous and global symmetries are spontaneously broken the collective modes known as the Nambu-Goldstone (NG) modes appear [4, 5].

Generally speaking, the single-particle excitations manifest themselves as poles of the single-particle Green’s function, , while the two-particle (collective) excitations could be related to the poles of the two-particle Green’s function, . The poles of these Green’s functions are defined by the solutions of the Schwinger-Dyson (SD) equation [6, 7] and the Bethe-Salpeter (BS) equation [8], respectively. Here, is the free single-particle propagator, is the fermion self-energy, is the BS kernel, and the two-particle free propagator is a product of two fully dressed single-particle Green’s functions. Since the fermion self-energy depends on the two-particle Green’s function, the positions of both poles must be obtained by solving the SD and BS equations self-consistently.

Instead of solving the SD and BS equations self-consistently, it is widely accepted that the single-particle dispersion can be obtained in the mean-field approximation or by solving the Bogoliubov-de Gennes (BdG) equations in a self-consistent fashion, while the generalized random phase approximation (GRPA) is the one that can provide the collective excitations in a weak-coupling regime. In the GRPA, the single-particle excitations are replaced with those obtained by diagonalizing the Hartree-Fock (HF) Hamiltonian, while the collective modes are obtained by solving the BS equation in which the single-particle Green’s functions are calculated in HF approximation, and the BS kernel is obtained by summing ladder and bubble diagrams.

From theoretical point of view, the corresponding expressions for Green’s functions cannot be evaluated exactly because the interaction part of the Hubbard Hamiltonian is quartic in the fermion fields. The simplest way to solve this problem is to apply the so-called mean-field decoupling of the quartic interaction. To go beyond the mean-field approximation, we apply the idea that we can transform the quartic term into quadratic form by making the Hubbard-Stratonovich transformation for the fermion operators. In contrast to the previous approaches, such that after performing the Hubbard-Stratonovich transformation the fermion degrees of freedom are integrated out; we decouple the quartic problem by introducing a model system which consists of a multicomponent boson field interacting with fermion fields and .

The functional-integral formulation of the Hubbard model requires the representation of the Hubbard interaction of (1) in terms of squares of one-body charge and spin operators. It is known that it may be possible to resolve the Hubbard interaction into quadratic forms of spin and electron number operators in an infinite number of ways. If no approximations were made in evaluating the functional integrals, it would not matter which of the ways is chosen. When approximations are taken, the final result depends on a particular form chosen. Thus, one should check that the results obtained with the Hubbard-Stratonovich transformation are consistent with the results obtained with the canonical mean-field approximation. It can be seen that our approach to the Hubbard-Stratonovich transformation provides results consistent with the results obtained with the mean-field approximation; that is, one can derive the mean-field gap equation using the collective-mode dispersion in the limit and .

There are three advantages of keeping both the fermion and the boson degrees of freedom. First, the approximation that is used to decouple the self-consistent relation between the fermion self-energy and the two-particle Green’s function automatically leads to conserving approximations because it relies on the fact that the BS kernel can be written as functional derivatives of the Fock and the Hartree self-energy . As it is known, any self-energy approximation is conserving whenever: (i) the self-energy can be written as the derivative of a functional , that is, , and (ii) the SD equation for needs to be solved fully self-consistently for this form of the self-energy. Second, the collective excitations of the Hubbard model can be calculated in two different ways: as poles of the fermion Green’s function, , and as poles of the boson Green’s function, , or equivalently, as poles of the density and spin parts of the general response function, . Here, the boson Green’s function, , is defined by the Dyson equation where is the free boson propagator. Third, the action which describes the interactions in the Hubbard model is similar to the action in quantum electrodynamics. This allows us to apply powerful field-theoretical methods, such as the method of Legendre transforms, to derive the SD and BS equations, as well as the vertex equation for the vertex function, , and the Dyson equation for the boson Green’s function, .

The mean-field treatment of the FF and LO phases in a variety of systems shows that the FF and LO states compete with a number of other states, such as the Sarma () states and the superfluid-normal separation phase (also known as the phase separation phase). It turns out that in some regions of momentum space the FF (or LO) phase provides the minimum of the mean-field expression of the Helmholtz free energy. Phase diagrams for a mixture at zero temperature were obtained in [9], but the calculations were limited to the emergence of insulating phases during the evolution of superfluidity from the BCS to the BEC regime, and the competition between the FF and Sarma phases was ignored. The polarization versus temperature diagrams in Figure 1 show that there are three phases: the Sarma phase, the FF phase, and the normal phase in which the Helmholtz free energy is minimized for gapless phase. The zero polarization line is the conventional Bardeen-Cooper-Schrieffer state. The phase diagram contains a Lifshitz point. When the interaction strength is increased, the Lifshitz point moves toward the higher temperatures and larger polarizations. Moreover, contrary to the phase diagram of population-imbalanced Fermi gas, where the phase separation appears for low polarizations, the existence of a polarization window for the FF phase was found. This means that as soon as the system is polarized it goes into the FF phase if the temperature is low enough. This polarization window is larger for a majority of atoms compared to the majority of atoms.

In what follows, the collective-mode dispersion of a mixture loaded in a two-dimensional optical lattice is calculated numerically by solving the BS equation in the GRPA approximation.

#### 2. Collective Modes of ^{6}Li-^{40}K Mixture

The superfluid states can be described in terms of the Namby-Gor’kov single-particle Green’s function which is a thermodynamic average of the -ordered tensor product of the four-component fermion fields:Here, we introduce composite variables and , where , are the lattice site vectors, and according to imaginary-time (Matsubara) formalism the variables , range from to . Throughout this paper we have assumed , the lattice constant , and we use the summation-integration convention: that repeated variables are summed up or integrated over. The field operators allow us to define the generalized single-particle Green’s function which includes all possible thermodynamic averages:

In the GRPA, the generalized Green’s function (3) is replaced by its mean field approximation. The Fourier transform of the mean-field single-particle Green’s function is as follows:

Here, the FF vector , as well as the chemical potentials and , and the gap are defined by the solutions of following set of four mean-field equations (the number equations, the gap equation, and the -equation) [10]: As was mentioned above, as soon as the system is polarized it goes into the FF phase, and therefore, it is natural to obtain the collective-mode dispersion using one point from the FF window shown in Figure 1. This point is defined by the following system parameters: , , , and . The solution of the mean-field equations (5) provides the FF vector , the gap, and the chemical potentials. Since the formation of the FF superfluid state is driven by the mass and population imbalance, which leads to the distortion of the Fermi surfaces. There are many equivalent ways to deform the surfaces, and as a result, the direction of the single FF pairing momentum is not specified. In what follows, we shall assume that the pairing momentum is along the -axis. The corresponding mean-field solutions are as follows: , , , and . For the Sarma state the corresponding mean-field values are , , and . The FF phase is the most stable as it provides the minimum of the mean-field expression of the Helmholtz free energy (the ratio between the FF free energy and the Sarma free energy is about 0.9986).

In the case when the order parameter is assumed to vary as a single plane wave, we have a broken translational invariance, and as a result, the normal and anomalous mean-field single-particle Green’s functions have phase factors associated with the FF quasimomentum , which can be eliminated using the corresponding unitary transformation. In the appendix, we have shown that in the tight-binding approximation the BS equation (in the GRPA) can be reduced to a secular determinant, which determines the collective-mode dispersion. When the generalized Green’s function (3) is used, the BS approach provides a secular determinant , defined in the appendix.

In Figure 2, we have presented the dispersion relation calculated for the system parameters listed in Section 1. The FF vector is directed along the -axis. The speed of sound, , to the positive and negative directions of the axis is defined by at . For the dispersions presented in Figure 2, we obtain in positive direction and in the negative direction. There are two roton minima at and at .

#### 3. Discussion

We have studied the Nambu-Goldstone excitation spectrum and the corresponding speed of sound of an interacting Fermi mixture of Lithium-6 and Potassium-40 atoms Fermi gases in deep optical lattices by using the Bethe-Salpeter equation in the GRPA. The generalized single-particle Green’s function, used in our numerical calculations, takes into account all possible thermodynamic averages. The chosen temperature and polarization correspond to a point from the polarization-temperature window for the FF phase in the phase diagram. The on-site attractive interaction corresponds to the weak-coupling regime, where the GRPA is valid.

It is well known that at very low temperatures, the singularities of the integrals (A.41) in the appendix correspond physically to the possibility of depairing into two fermion excitations. At a zero temperature, the two fermion excitations have energies and . The spectrum for this kind of excitation is known as the particle-hole continuum. At a zero temperature, the lower boundary of the particle-hole continuum is defined by the condition where the minimum is to be taken over all the possible values of . The collective-mode dispersion, presented in Figure 2, lies below the particle-hole continuum, and therefore, the possibility of depairing into two fermion excitations is not important.

It is known that the superfluid ground state comes from the symmetry breaking and it is characterized by a nonzero order parameter. If the atoms have different masses and chemical potentials, the superfluid states of Cooper pairs are with nonzero momentum. A finite pairing momentum implies a position-dependent phase of the order parameter. Based on general symmetry principles, it is known that when continuous and global symmetries are spontaneously broken in Lorentz-invariant systems, the number of Nambu-Goldstone modes is always equal to the number of broken generators, and all of them have linear dispersion [11, 12]. In the FF and in the LO cases, we are dealing with a spontaneous translational symmetry breaking and with an inhomogeneous superfluid state; but in case of the FF order parameter, only one generator is spontaneously broken, because there is an unbroken combination of rotations and translations. Thus, the FF superfluid state is characterized by a single Nambu-Goldstone mode, and this statement is in accordance with the fact that the determinant in the appendix provides only one mode. It is worth mentioning that in the case of the LO order parameter, both the particle number symmetry and translations along a given direction are spontaneously broken, and therefore, there should exist two Nambu-Goldstone modes. The Bethe-Salpeter description of the collective modes in the case of LO superfluid states is expected to be more complicated because the phase factors of the nondiagonal elements of the single-particle Green’s function (see (A.32) in the appendix) one cannot eliminate as in the case of FF states.

In view of the fact that most of the previous numerical calculations are based on the Nambu-Gor’kov single-particle Green’s function,which leads to the secular determinant, one may well ask whether the collective-mode dispersion, defined by determinant, is significantly different in comparison to the dispersion obtained with the secular determinant .

To answer this question, we have calculated the collective-mode dispersion relation of a population-imbalanced atomic Fermi gas. The dispersions obtained by means of the secular determinant and are presented in Figure 3. In the range of small , we find a difference of about 15% between the speed of sound obtained by means of secular determinant and the speed of sound, calculated by secular determinant used in [13]. The figure also indicates that the difference between the two dispersion curves tends to increase with , reaching 25% at .

#### Appendix

#### Bethe-Salpeter Approach to the Collective Modes of ^{6}Li-^{40}K Mixture

Green’s functions in the functional-integral approach are defined by means of the so-called generating functional with sources for the boson and fermion fields. In our problem, the corresponding functional integrals cannot be evaluated exactly because the interaction part of the Hamiltonian (1) is quartic in the Grassmann fermion fields. However, we can transform the quartic terms to a quadratic form by introducing a model system which consists of a four-component boson field (, , ) interacting with fermion fields and . The action of this model system is assumed to be of the following form , where

The action describes the fermion part of the system. The generalized inverse Green’s function of free fermions is given by the following diagonal matrix:where and . The symbol is used to denote (for fermion fields ; ). In the case of the FF states of a population-imbalanced Fermi gas, the noninteracting Green’s function iswhere .

The action describes the boson field which mediates the fermion-fermion on-site interaction in the Hubbard Hamiltonian. The bare boson propagator in is defined as The Fourier transform of the boson propagator is given by

The interaction between the fermion and the boson fields is described by the action . The bare vertex is a matrix, whereThe Dirac matrix and the matrices are defined as (when a four-dimensional space is used, the electron spin operators have to be replaced by [14])

The relation between the Hubbard model and our model system can be demonstrated by applying the Hubbard-Stratonovich transformation for the fermion operators:The functional measure is chosen to be

According to the field-theoretical approach, the expectation value of a general operator can be expressed as a functional integral over the boson field and the Grassmann fermion fields and :where the symbol means that the thermodynamic average is made. The functional is defined by where the functional measure satisfies the condition . The quantity is the source of the boson field, while the sources of the fermion fields are included in the term:Here, we have introduced complex indices , and .

We shall now use a functional derivative ; depending on the spin degrees of freedom, there are sixteen possible derivatives. By means of definition (A.10), one can express all Green’s functions in terms of the functional derivatives with respect to the corresponding sources of the generating functional of the connected Green’s functions . Thus, we define the following Green’s and vertex functions which will be used to analyze the collective modes of our model.

The Boson Green’s function is which is a matrix defined as .

The generalized single-fermion Green’s function is matrix (3) whose elements are . Depending on the two spin degrees of freedom, and , there exist eight “normal” Green’s functions and eight “anomalous” Green’s functions. We introduce Fourier transforms of the “normal” , and “anomalous” one-particle Green’s functions, where . Here and are the creation-annihilation Heisenberg operators. The Fourier transform of the generalized single-particle Green’s function is given by Here, and are matrices whose elements are and , respectively.

The two-particle Green’s function is defined as

The vertex function for a given is a matrix whose elements are

Next, we shall obtain the corresponding equations of the boson and fermion Green’s functions. The poles of these Green’s functions provide the single-particle and the two-particle excitations.

It is well known that the fermion self-energy (fermion mass operator) can be defined by means of the SD equations. They can be derived using the fact that the measure is invariant under the translations and : where is the average boson field. The fermion self-energy is a matrix which can be written as a sum of Hartree and Fock parts. The Hartree part is a diagonal matrix whose elements are

The Fock part of the fermion self-energy is given byThe Fock part of the fermion self-energy depends on the two-particle Green’s function ; therefore the SD equations and the BS equation for have to be solved self-consistently.

Our approach to the Hubbard model allows us to obtain exact equations of Green’s functions by using the field-theoretical technique. We now wish to return to our statement that Green’s functions are the thermodynamic average of the -ordered products of field operators. The standard procedure for calculating Green’s functions is to apply Wick’s theorem. This enables us to evaluate the -ordered products of field operators as a perturbation expansion involving only wholly contracted field operators. These expansions can be summed formally to yield different equations of Green’s functions. The main disadvantage of this procedure is that the validity of the equations must be verified diagram by diagram. For this reason we will use the method of Legendre transforms of the generating functional for connected Green’s functions. By applying the same steps as in [15], we obtain the BS equation of the two-particle Green’s function, the Dyson equation of the boson Green’s function, and the vertex equation:Here, is the two-particle free propagator constructed from a pair of fully dressed generalized single-particle Green’s functions. The kernel of the BS equation can be expressed as a functional derivative of the fermion self-energy . Since , the BS kernel is a sum of functional derivatives of the Hartree and Fock contributions to the self-energy:The general response function in the Dyson equation (A.20) is defined asThe functions , , and are related by the identity:

By introducing the boson proper self-energy one can rewrite the Dyson equation (A.20) for the boson Green’s function as The proper self-energy and the vertex function are related by the following equation:It is also possible to express the proper self-energy in terms of the two-particle Green’s function which satisfies the BS equation , but its kernel includes only diagrams that represent the direct interactions:One can obtain the spectrum of the collective excitations as poles of the boson Green’s function by solving the Dyson equation (A.26), but one has first to deal with the BS equation for the function . In other words, this method involves two steps. For this reason, it is easy to obtain the collective modes by locating the poles of the two-particle Green’s function using the solutions of the corresponding BS equation.

As we have already mentioned, the BS equation and the SD equations have to be solved self-consistently. In what follows, we use an approximation which allows us to decouple the above-mentioned equations and to obtain a linearized integral equation for the Fock term. To apply this approximation we first use (A.25) to rewrite the Fock term asand after that we replace and in (A.29) by the free boson propagator and by the bare vertex , respectively. In this approximation the Fock term assumes the form:The total self-energy is , whereThe contributions to , due to the elements on the major diagonal of the above matrices, will be included into the chemical potential. To obtain an analytical expression for the generalized single-particle Green’s function, we assume two more approximations. First, since the experimentally relevant magnetic fields are not strong enough to cause spin flips, we shall neglect . Second, we neglect the frequency dependence of the Fourier transform of the Fock part of the fermion self-energy. Thus, the Dyson equation for the generalized single-particle Green’s function becomesWe can eliminate the phase factors by performing the unitary transformation between the old generalized single-particle Green’s function and the new one ; that is, , where the corresponding matrix is as follows: After performing this unitary transformation, Green’s function becomes functions of , and the corresponding Fourier transform isThe spectrum of the collective modes will be obtained by solving the BS equation in the GRPA. As we have already mentioned, the kernel of the BS equation is a sum of the direct and exchange interactions, written as derivatives of the Fock (A.30) and the Hartree (A.31) parts of the self-energy. Thus, in the GRPA the corresponding equation for the BS amplitude can be obtained by performing integration over the momentum vectors: where the two-particle propagator and the direct and exchange interactions are defined as follows:

The BS equation (A.35) written in the matrix form is , where is the unit matrix, the matrix is a matrix, and the transposed matrix of is given by

The secular determinant can be rewritten as a block diagonal determinant: where is a unit matrix. The block structure of the secular determinant allows us to separate the sixteen BS amplitudes into three independent groups related to the blocks , , and . The determinantof the block determines the amplitudes , , , , , , , and . The other four amplitudes , , , and are related to the block. The last four amplitudes are equal to zero; that is, . The collective-mode dispersion is defined by the secular determinant (A.39). At a finite temperature the elements of arewhere the symbols , , , and are defined asHere, , , , and and are one of the following form factors: