Speed of the Goldstone Sound Mode of an Atomic Fermi Gas Loaded on a Square Optical Lattice with a Non-Abelian Gauge Field in the Presence of a Zeeman Field

The speed of the Goldstone sound mode of a spin-orbit-coupled atomic Fermi gas loaded in a square optical lattice with a nonAbelian 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 soundmode 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.


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 twocomponent pseudospin-↑, ↓ fermionic gas loaded on a square optical lattice with a non-Abelian gauge field A = ( 푦 , − 푥 ) 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 threedimensional free space [7][8][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 swave 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][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26].We restrict the discussion to the case of atoms confined to the lowest-energy band 2 Advances in Condensed Matter Physics (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 Ĥ = Ĥ0 + Ĥ푈 , where The out-of-plane Zeeman field is described by the term Ĥ푍 : ) . ( The Hubbard part of the Hamiltonian is Here,  is the chemical potential, and n푖 ,휎 =  † 푖,휎  푖,휎 is the density operator on site .The Fermi operator is the tight-binding energy,  0 is the identity matrix, and  푥,푦 (k) characterize the SOC.
In the case of a square lattice the vectors k = ( 푥 ,  푦 ) (in units  −1 ) 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) as where 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 T푢 -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 T푢 -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  −1 =  (0)−1 − Σ for the poles of the single-particle Green's function, , and the Bethe-Salpeter [30,31] (BS) equation [ (0)−1 − ]Ψ = 0 for the poles of the poles of the two-particle Green's function, .Here,  (0)  is the free single-particle propagator, Σ is the fermion selfenergy,  is the BS kernel, and the two-particle free propagator  (0) =  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.

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 − ∑ 푖 n푖 ,↑ n푖 ,↓ 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 6 -40  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 functionalintegral 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 fourcomponent boson field  훼 () interacting with fermion fields ψ() = Ψ † ()/ √ 2 and ψ() = Ψ()/ √ 2, where Here  = 1, 2, 3, 4,  = (r 푖 , ),  = (r 푗 ,  耠 ), and  = (r 푘 , V) are composite variables, where r 푖 , r 푗 , and r 푘 are the lattice site vectors.According to imaginary-time (Matsubara) formalism the variables 0 ≤ V, ,  耠 ≤  = ( 퐵 ) −1 ;  퐵 is the Boltzmann constant (ℏ =  퐵 = 1).
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 singleparticle Green function: ) , (10) where , 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.

Mean Field Approximation.
The Fourier transform of zero-temperature single-particle Green's function (10) in the mean field approximation is The matrix elements Ĝ푀퐹 where the corresponding expressions for , and  푛 1 푛 2 (k) can be obtained by inverting matrix (11) and setting  → 0.
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, (k) =  2 푥 (k) +  2 푦 (k), () 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  ↑,(↓) = ∑ k ∑ 횤휔   푀퐹  11,(22) (k,  푚 ) 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 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

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 selfenergy.Thus, in the GRPA the corresponding equation for the BS amplitude Ψ Q 푛 2 ,푛 1 is given by (see [32]) where ) are the direct and the exchange interactions, correspondingly.The two-particle propagator  (0) in the GRPA is defined as follows: The BS equation ( 21) can be written in the matrix form as ( Î +  Ẑ) Ψ = 0, where Î is the unit matrix, the matrix Ẑ is a 16×16 matrix, and the transposed matrix of Ψ is given by At a given vector Q from the Brillouin zone, the collective excitation spectrum (Q) is obtained by the condition det | Î+  M((Q), Q))| = 0 to have a nontrivial solution of the BS equation.By applying simple matrix algebra, the above secular 16 × 16 determinant can be simplified to a 10 × 10 BS determinant det |((Q), Q)| = 0, where the BS matrix is given by The elements  푖,푗 ,  푖,푗 ,  푖,푗 , and  푖,푗 of the blocks  2×2 (, Q),  2×8 (, Q),  8×2 (, Q), and  8×8 (, Q) 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 (Q) is defined by the 2 × 2 Gaussian secular determinant while within the BS formalism the secular determinant is more complicated (((Q), Q) dependence is understood):

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  = 0.4 and  = 0.6.
The strength of the attractive interaction is / = 4, and for both filling factors the polarization is fixed to be  = 0.48.
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  = /2, where (k) = − and (k) = 4 2 [sin 2 ( 푥 ) + sin 2 ( 푦 )], 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  13)-( 17) using two different filling factors:  = 0.4 and  = 0.6.The on-site attractive strength is  = 4, and the polarization is  = 0.48. has been numerically obtained using the collective-mode dispersion (Q) for small wave vectors in  푥 -direction.
(k 0,1 ) = 0.The first one is at the center of the Brillouin zone k 0 = (0, 0), and (k 0 ) = 0 for / between 0.39 and 0.40.At the second point, k 1 = (, ), we found also the existence of opening and closing of a gap in momentum space when 0.61 < / < 0.62.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][8][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 ℎ < ℎ 푐 = √ ((k 0 ) − ) 2 + Δ 2 , the slope is suppressed by the Zeeman field, while in the nontrivial side, where ℎ > ℎ 푐 = √ ((k 0 ) − ) 2 + Δ 2 , 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  = 0.38, and we have accessed the opening and the closing of a gap in momentum space at the point k 0 = (0, 0) 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 ℎ 푐 / = 1.39 and ℎ 푐 / = 1.49.The insert shows the behavior of the slope close to ℎ 푐 / = 1.49.We found that there is no sharp change in the slope across the two TQPT points.

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 (k) = −2 cos()[cos( 푥 ) + cos( 푦 )] 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 nearestneighbor 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 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.

Figure 1 :
Figure 1: Chemical potential , pairing gap Δ, and the Zeeman field ℎ of a Fermi gas in a square optical lattice subject to a non-Abelian gauge field A = ( 푦 , − 푥 ).The graphs are obtained by solving the mean field equations (13)-(17) using two different filling factors:  = 0.4 and  = 0.6.The on-site attractive strength is  = 4, and the polarization is  = 0.48.
Figure2: The singlet  푠 and the triplet  푡푟 contributions to the total condensed fraction  =  푠 +  푡푟 of a Fermi gas in a square optical lattice subject to a non-Abelian gauge field A = ( 푦 , − 푥 ).The graphs are obtained by solving the mean field equations along with definitions(19) and(20).The system parameters are the same as in Figure1.

Figure 3 :
Figure3: The slope of the Goldstone sound mode of a Fermi gas in a square optical lattice subject to a non-Abelian gauge field A = ( 푦 , − 푥 ).The graphs are obtained by applying the Gaussian and the BS secular determinants.The system parameters are the same as in Figure1.The arrows show the values of  for which the opening and closing of a gap take place.

Figure 4 :
Figure 4: Chemical potential , pairing gap Δ, and the slope of the Goldstone sound mode of a Fermi gas in a square optical lattice subject to a non-Abelian gauge field with  = 0.38 as a function of the Zeeman field ℎ.The chemical potential and the gap are obtained by solving the mean field number and gap equations (13)-(14).The arrows show the values of the Zeeman field for which the opening and closing of a gap take place.The insert shows the slope close to ℎ 푐 / = 1.49.