#### Abstract

The speed of the Goldstone sound mode of a spin-orbit-coupled atomic Fermi gas loaded in a square optical lattice with a non-Abelian gauge field in the presence of a Zeeman field is calculated within the Gaussian approximation and from the Bethe-Salpeter equation in the generalized random phase approximation. It is found that (i) there is no sharp change of the slope of the Goldstone sound mode across the topological quantum phase transition point and (ii) the Gaussian approximation significantly overestimates the speed of sound of the Goldstone mode compared to the value provided by the BS formalism.

#### 1. Introduction

It is well known that topological superfluids are new states of matter that can be observed in two-dimensional atomic Fermi gases with strong Rashba spin-orbit coupling (SOC), conventional s-wave pairing, and out-of-plane magnetic field (Zeeman field) which breaks time-reversal symmetry. In what follows, we study the speed of sound of a two-component pseudospin- fermionic gas loaded on a square optical lattice with a non-Abelian gauge field in the presence of a Zeeman field, where is independently tunable parameter, and are the Pauli matrices. The external non-Abelian gauge field effectively creates a SOC in the gas. The pseudospin of atoms can couple with the Zeeman field and the orbital degrees of freedom of atoms, which gives rise to a topological quantum phase transition (TQPT) between gapped (topological trivial) and gapless (topological nontrivial) superfluid states. In other words, the spectrum of the Hamiltonian becomes gapless at some set of points in the momentum space. Due to absence of symmetry breaking across the TQPTs, it remains experimentally challenging to detect these phase transitions. It is experimentally possible to detect TQPTs by measuring the dynamic response of the bulk under an external force [1], by observing the edge modes [2, 3], or the momentum distributions of atoms across the phase transition boundary using a noise-correlation imaging [4] and momentum-resolved spectroscopy [5, 6].

Recent numerical calculations of the sound speed of atomic Fermi gases with s-wave attraction, synthetic Rashba SOC, and out-of-plane Zeeman field in two- and three-dimensional* free space* [7–9] show sharp changes of the slope of the Goldstone sound mode across the phase transition points, and this can be used to detect TQPTs. Our goal is to study whether a similar behavior of the sound speed occurs in the lattice case if the SOC is created by the non-Abelian gauge field. We found that the sharp change in the slope does not exist in the presence of lattice geometry.

It is widely accepted among the cold-atom community that the system of a two-component Fermi gas with s-wave attraction in square optical lattice with an external non-Abelian gauge field (or a synthetic Rashba SOC) and out-of-plane Zeeman field can be described by a negative-U Hubbard model [10–26]. We restrict the discussion to the case of atoms confined to the lowest-energy band (single-band model). We assume that there are atoms distributed along sites, and the filling factor is smaller than unity. The corresponding polarization is . The Hamiltonian for an uniform system is , whereThe out-of-plane Zeeman field is described by the term :The Hubbard part of the Hamiltonian isHere, is the chemical potential, and is the density operator on site . The Fermi operator creates (destroys) a fermion on the lattice site with pseudospin projection . The sum means sum over the nearest-neighbor sites of the 2D lattice. describes the hoping of the atoms between site and site , and is accumulated via the Peierls substitution phase factor. The strength of the on-site Hubbard interaction is , which corresponds to attractive interaction.

Hamiltonian (1) in a momentum representation assumes the form . Here, , , is the tight-binding energy, is the identity matrix, and characterize the SOC.

In the case of a square lattice the vectors (in units ) are restricted to be within the Brillouin zone, i.e., a square defined by the following four points . The tight-binding energy and are given by where is the strength of the nearest-neighbor hopping (in what follows, the lattice spacing and spacing are set to unity).

The opening or closing of the gap can be achieved by varying the tunable parameter along with the Zeeman field. We first assume that the opening or closing of the gap is achieved by varying the parameter , while by changing the out-of-plane Zeeman field the polarization is kept fixed. We also discuss the second possibility, when one can close or open the gap by changing the Zeeman field while keeping the parameter fixed. In this case the polarization changes during the transition from gapless to gapped superfluid states.

The single-particle and collective excitations of the system under consideration manifest themselves as poles of the single-particle and two-particle Green’s functions, but the corresponding expressions for the 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. According to the standard mean field theory on the Hubbard Hamiltonian, the on-site interaction is decoupled (up to Hartree-Fock correction terms) aswhere the order parameter is a real number.

Instead of applying the mean field decoupling, we shall transform the quartic terms to a quadratic form by introducing a four-component boson field which mediates the interaction of fermions in the same manner as in quantum electrodynamics, where the photons mediate the interaction of electric charges. Green’s functions are thermodynamic averages 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 perturbation expansions 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 shall use the method of Legendre transforms [27] of the generating functional for connected Green’s functions to derive the Schwinger-Dyson [28, 29] (SD) equation for the poles of the single-particle Green’s function, , and the Bethe-Salpeter [30, 31] (BS) equation for the poles of the poles of the two-particle Green’s function, . 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. The kernel of the BS equation is defined as a sum of the direct interaction, , and the exchange interaction , where and are the Fock and the Hartree parts of the fermion self-energy . Since the fermion self-energy depends on the two-particle Green’s function, the positions of both poles must to be obtained by solving the SD and BS equations self-consistently. Instead of solving self-consistently the SD and BS equations, we shall employ the generalized random phase approximation (GRPA). In this approximation, the single-particle excitations are obtained in the mean field approximation, while the collective modes are obtained by solving the BS equation in which single-particle Green’s functions are calculated in Hartree-Fock approximation, and the BS kernel is obtained by summing ladder and bubble diagrams.

#### 2. Functional-Integral Formalism

##### 2.1. Field-Theoretical Approach to the Hubbard Model

The functional-integral formulation of the Hubbard model requires the representation of the Hubbard interaction in terms of squares of fermion operators. This can be done by employing a certain Hubbard-Stratonovich transformation. This field-theoretical approach to the Hubbard Hamiltonian has already been used to describe the collective modes of ultracold - mixture in a square optical lattice (see [32]). Hamiltonian (1) in the present study has two more interactions. It is possible to include both the SOC and the Zeeman field into initial definition of the free fermion Green’s function, and therefore most of the general equations in [32] remain unchanged. In particular, the general expressions for the SD and BS equations are the same as those without the SOC and the Zeeman field, but with more complicated fermion self-energy and BS kernel.

As it is well known, 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, but the corresponding functional integrals cannot be evaluated exactly because the interaction part of the Hubbard Hamiltonian 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 , whereHere , , , and are composite variables, where , and are the lattice site vectors. According to imaginary-time (Matsubara) formalism the variables ; is the Boltzmann constant ().

The field operators (6) allow us to define generalized single-particle Green’s function by using a tensor product of these two matrices. The corresponding Green function, represented by a matrix, includes all possible thermodynamic averages:

The action of the above-mentioned model system is assumed to be of the following form , where , , . The fermion action corresponds to the Hamiltonian . The corresponding inverse Green function of free fermions is given by the following matrix:The symbol is used to denote (for fermion fields ; ).

The Fourier transform of noninteracting Green’s function is After including the Zeeman field and the SOC interaction into free fermion Green’s function (9), one can perform the same steps as in [32] to come up with the following single-particle Green function:where , and , where is a constant in space. Physically, it describes a superfluid state of Cooper pairs with zero momentum when the pairing is only between atoms with equal and opposite momenta.

##### 2.2. Mean Field Approximation

The Fourier transform of zero-temperature single-particle Green’s function (10) in the mean field approximation is

The matrix elements () of are given bywhere the corresponding expressions for , , , and can be obtained by inverting matrix (11) and setting .

For a fixed filing factor and at a given Zeeman field , the chemical potential and the gap have to be determined by solving the mean field number and gap equations:Here, , is the Fermi-Dirac distribution function, and the following notations have been introduced:The number equation (13) follows from the definitions of the filing factors through corresponding mean field Green’s functions. It is known that the gap equation can be obtained by minimizing the mean field Helmholtz free energy with respect to .

As we mentioned before, in order to keep the polarization fixed, one has to adjust the corresponding Zeeman field. This means that , , and have to be calculated by solving the number and the gap equations, along with the following equation:

Another interesting feature is that because of the SOC term in Hamiltonian (1) the pairing field contains both singlet and triplet component. The singlet and triplet , amplitudes, obtained by means of Green’s function elements , , and , are

With the help of the pairing amplitudes, one can calculate the total condensate fraction , where and are the singlet and the triplet contributions, correspondingly. At zero temperature, we obtain

##### 2.3. The Bethe-Salpeter Equation for the Collective Excitations in the Generalized Random Phase Approximation

The 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 and the Hartree parts of the self-energy. Thus, in the GRPA the corresponding equation for the BS amplitude is given by (see [32]) where and are the direct and the exchange interactions, correspondingly. The two-particle propagator in the GRPA is defined as follows:The BS equation (21) can be written in the matrix form as , where is the unit matrix, the matrix is a matrix, and the transposed matrix of is given by

At a given vector from the Brillouin zone, the collective excitation spectrum is obtained by the condition to have a nontrivial solution of the BS equation. By applying simple matrix algebra, the above secular determinant can be simplified to a BS determinant , where the BS matrix is given byThe elements , and of the blocks , , , and are given in the Appendix in terms of the propagators .

In contrast to our functional-integral formalism, one can use the Hubbard-Stratonovich transformation to introduce the energy gap as an order parameter field. In this approach, one can integrate out the fermion fields and to arrive at an effective action. The next steps are to consider the state which corresponds to the saddle point of the effective action and to write the effective action as a series in powers of the fluctuations and their derivatives. The exact result can be obtained by explicitly calculating the terms up to second order in the fluctuations and their derivatives. This approximation is called the Gaussian approximation. Within this approximation, the collective excitation spectrum is defined by the Gaussian secular determinantwhile within the BS formalism the secular determinant is more complicated ( dependence is understood):

#### 3. Numerical Results

There are various mean field quantities of physical interest, such as the chemical potential, the pairing gap, the singlet and triplet pairing amplitudes, and the singlet and triplet condensate fractions. We focus on the zero-temperature case assuming two different filling factors of and . The strength of the attractive interaction is , and for both filling factors the polarization is fixed to be .

In Figure 1, we plot the chemical potential , the gap , and the corresponding Zeeman field , as functions of the SOC parameter , obtained by solving the mean field equations (13)-(17). In Figure 2 we have shown the singlet and the triplet condensate fractions, as functions of the SOC parameter , obtained by solving the mean field equations along with definitions (19) and (20). It can be seen that at the value , where and , the chemical potential has a minimum, while the gap and the Zeeman field reach their maximum values. As expected, the singlet and triplet condensate fraction-graphs have maximum at the same as the gap.

Next, we have calculated the speed of sound as a function of the SOC parameter . The sound speed is proportional to the slope of the Goldstone sound mode . The slope has been numerically obtained using the collective-mode dispersion for small wave vectors in -direction.

It is clear that the single-particle dispersion, given by (16), becomes precisely zero at some special points in momentum space satisfying the condition when . The condition is met at the following points: , , , and . For , , and a filling factor , we have two points at which the single-particle dispersion . The first one is at the center of the Brillouin zone , and for between and . At the second point, , we found also the existence of opening and closing of a gap in momentum space when . Therefore, at both points one should have a TQPT between superfluid states with different topologies. During the TQPT the different phases of matter are not characterized by an order parameter, but rather an integer number, the Chern number, describing the system as a whole. This is because the different phases across the TQPT are characterized by different topologies of Fermi surface instead of being classified by different symmetries.

In Figure 3, we have shown the slope calculated by using the Gaussian and the BS secular determinants (25) and (26), correspondingly. It is worth mentioning the significant difference between the slope obtained within the Gaussian approximation and the slope obtained by the BS secular determinant. The arrows show the points where the opening and closing of a gap take place. As can be seen, the sharp change in the slope does not exist in the presence of a non-Abelian gauge field.

The boundaries for the TQPT in [7–9] were accessed by varying the Zeeman filed at a constant strength of the synthetic Rashba SOC. It was found that on a topological trivial side, where , the slope is suppressed by the Zeeman field, while in the nontrivial side, where , the slope is enhanced by the Zeeman field. As mentioned in [7, 8], this is quite unusual since we generally expect that the superfluidity should be suppressed by the Zeeman field.

We have fixed the non-Abelian parameter , and we have accessed the opening and the closing of a gap in momentum space at the point by varying the Zeeman filed. To do this, first we used (13)-(14) to calculate numerically the chemical potential and the gap as functions of the Zeeman field. The results are shown in Figure 4. The TQPT takes place for and . The insert shows the behavior of the slope close to . We found that there is no sharp change in the slope across the two TQPT points.

#### 4. Discussion

In this paper, we are concerned with the question which naturally arises here as to whether there is a possibility of detecting TQPTs by monitoring the behavior of slope of the sound mode if a non-Abelian gauge field is used instead of a synthetic Rashba SOC. The answer is not trivial because the non-Abelian gauge field is taken into account via the Peierls substitution. As a result, the tight-binding energy does depend on the parameter . The second difference is that the nondiagonal part of , i.e., , describes the spin-flipped hopping between site and site . In the case of a synthetic Rashba SOC we do not have spin-flipped hopping, and the hopping to the nearest-neighbor site does not depend on the strength of Rashba SOC.

In summary, we have derived the BS equation in the GRPA for the collective excitation energy of a Fermi gas loaded on a square optical lattice with a non-Abelian gauge field in the presence of a Zeeman field. To the best of our knowledge, there is no other calculation of the speed of sound in a lattice geometry as a function of the non-Abelian parameter with which we can make comparisons. According to our numerical calculations, there is no sharp change of the slope of the Goldstone sound mode across the phase transition point. It is found that the Gaussian approximation significantly overestimates the speed of sound of the Goldstone mode compared to the value obtained within the BS formalism.

#### Appendix

#### Blocks of the BS Matrix

The two-particle propagator in BS equation (21) is defined as

The elements , , , and of the blocks , , , and are given in terms of the propagators as follows:

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

Israel Chávez Villalpando gratefully acknowledges funding from CONACyT Grant no. 291001 and UNAM-DGAPA-PAPIIT IN102417.