Phase Diagram of the Attractive Kane-Mele-Hubbard Model at Half Filling

: Motivated by recent developments in the experimental study of ultracold atoms in graphene-like honeycomb optical lattices, we investigate superconductivity of the attractive Kane-Mele-Habbard (KMH) model with the next-nearest-neighbor (NNN) hoping at half ﬁlling. The mean-ﬁeld approximation is used to study the phase diagram which interpolates the trivial and the non-trivial topological states. It is shown that: (a) when the NNN hoping is taken into account, one has to introduce two mean-ﬁeld gap equations for the two sublattices, instead of a single gap when the NNN hopping is neglected, and (b) in the non-trivial topological region the phase diagram with the NNN hopping is signiﬁcantly different compared to the phase diagram calculated previously, but without the NNN term. We also discuss the superconducting instability of the attractive KMH model that is driven by condensation of Cooperons.


Introduction
The present-day experiments with ultracold atoms in optical lattices allow for us to simulate both the Haldane's model [1] and the situation [2,3] considered by Kane and Mele (KM). The Haldane's model [4,5] is a tight-binding representation of electron motion on a honeycomb lattice in the presence of a magnetic field, where vector potential has the full symmetry of the lattice and generates a magnetic field with zero total flux through the unit cell. It was pointed out by Haldane that in two-dimensional (2D) honeycomb lattice the topological ordering requires time reversal symmetry breaking. Because of the zero magnetic flux through each unit cell, the phase accumulated through a nearest neighbor hopping vanishes, whereas the phase accumulated through NNN hopping is nonzero. This extra phase breaks the time-reversal symmetry. The Haldane model has been experimentally realized with ultracold atoms [6,7]. However, the electron spin is not included in the Haldane's model.
It is known that the spin-orbit coupling preserves the time-reversal symmetry, but the spin-orbit effects can be used to obtain topological insulators. The first example was the KM Hamiltonian for the electrons in a graphene [8,9], which consists of two copies of the Haldane's model, one for spin-up electrons and one for spin-down electrons. In the KM model, each spin component breaks time-reversal symmetry, but the time reversal symmetry is restored when taking two copies with different signs for the spin together.
In a recent paper [20], the phase diagram of the attractive Kane-Mele-Hubbard (KMH) model at half filling has been obtained as a function of a tuning parameter x = 3 √ 3λ/(m AB + 3 √ 3λ). Here, λ is the strength of intrinsic spin-orbit (ISO) coupling and m AB is sublattice potential. However, the NNN hopping term has been neglected, although it is several orders of magnitude stronger than the ISO coupling.
In this paper, we have presented the phase diagram of the attractive KMH model with NNN hoping at half filling as a function of a tuning parameter where t is the NNN hoping amplitude. It is shown that, in the mean-field approximation, we have two gap equations for ∆ A and ∆ B , instead of a single gap ∆ when NNN hopping is neglected.
We shall discuss the case of Fermi (spin-1/2) atoms loaded into honeycomb optical lattice, but the results are also valid for the tight-binding description of electrons in a graphene as well. Our tight-binding Hamiltonian H = H KM + H NNN + H U includes the KM terms, as well as the NNN hopping and the onsite attractive Hubbard interaction, where The first term in (1) takes into account the possibility for nearest-neighbor hopping. H ISO represents the ISO interaction, which originates from the hybridization of the higher angular momentum orbit and exerts opposite magnetic fields upon electrons with opposite spin polarizations. The third term in (1) describes the possible energy offset between sites of A and B sublattices. t is the nearest-neighbor hopping amplitudes, m AB is the energy offset parameter, and is the creation operator for spin-up and spin-down fermions at site i. Throughout this paper, we have assumed the lattice constant a = 1. The hopping to the nearest neighbor sites the vectors j are: Figure 1a). The NNN hopping term is . For hopping to the next nearest neighbor sites the vectors j are d 1,2 = ±(3/2, √ 3/2), d 3,4 = ±(3/2, − √ 3/2), and d 5,6 = ±(0, √ 3). The attractive Hubbard interaction is described by H U = −U ∑ i n i,↑ n i,↓ , where U > 0. The KM Hamiltonian along with the NNN hopping term are represented by the following 4 × 4 matrix in the momentum space on the basis of the four-component wave function where , and λ is the strength of the ISO interaction. The eigenvalues of the Hamiltonian (2) are The authors of Ref. [20] pointed out that, depending on the value of the parameter 0 ≤ x ≤ 1, the system can be taken across the topological phase transition. For x < 1/2, we have a topological trivial insulator, while, for x > 1/2, the system is in topological non-trivial insulator state. The phase diagram of the KMH model, as reported in Ref. [20], was calculated within the self-consistent Bogoliubov-de Gennes theory while using a supercell with six sites. The corresponding system parameters are In what follows, we apply the mean-field approximation to obtain the matrix elements of the single-particle Green's functions of the KMH model with NNN hopping term at half filling. The above-mentioned supercell approach is more complicated than the mean-field approximation, but, as can be seen from Figure 2a, in this paper, our approach reproduces the corresponding phase diagram (Figure 2a in Ref. [20]). If the NNN term is taken into account, then the matrix elements of the mean-field single-particle Green's function in the momentum space, proportional to < Ψ Ak↓ Ψ A−k↑ > and < Ψ Bk↓ Ψ B−k↑ >, are different ands, therefore, we have to introduce two gaps, ∆ A , and ∆ B . If the NNN hopping is neglected, we have ∆ A = ∆ B = ∆. When the NNN hopping is included, we have m AB = −3t + (1 − x)E g , and 3 √ 3λ = xE g . At 3t + m AB = 3 √ 3λ the the single-particle gap does , while the mass of the other bands at K = 2π 3a , − 2π remains constant throughout the transition for all values of x, and viceversa. As in Ref. [20], the ground state of the KM Hamiltonian with NNN hopping is topological nontrivial (or topological trivial), when the parameter . Solid lines mark topological phase transitions between the topological non-trivial and trivial insulating states, and the s-wave spin-singlet superfluid state: (a) without the NNN term (t = 0), and (b) with (t = −0.15t).

Single-Particle Dispersion in the Mean-Field Approximation
In the presence of an onsite attractive interaction between the fermions, the fermion atoms form bound (Cooper) pairs. As a result, the system becomes unstable against the formation of a s-wave spin-singlet superfluid ground state. At low energies, the system admits an effective description in terms of massless Dirac fermions; therefore, in a honeycomb optical lattice we have a possibility to observe a superfluidity of Fermi atoms with the Dirac spectrum. We restricted our calculations to half-filling (the chemical potential µ = 0), where the particle-hole symmetry takes place. We further assume that the BCS mean-field order parameters ∆ A(B) = U < ψ A(B),−k,↓ ψ A(B),k,↑ > are real constants. When the attractive Hubbard interaction is taken into account, the KM basis of the four-component wave function where the corresponding 4 × 4 blocks are defined in the Appendix A.
We define the Matsubara single-particle Green's functions: Here, the Matsubara fermion energies are ω m = (2π/β)(m + 1/2), m = 0, 1, 2, . . . , β = 1/(k B T), and k B is the Boltzmann constant (throughout this paper we have assumedh = k B = 1). The single-particle excitations ω i (k), i = 1, 2, 3, 4 in the mean-field approximation are defined in the Appendix A. They manifest themselves as poles of the Green's function G(k, ıω m ). The corresponding zero-temperature Green's function G(k, ω) is an 8 × 8 matrix with elements G (1,2) n 1 ,n 2 (k, ω) {n 1 , n 2 } = 1, 2, . . . , 8 written in the following forms: The functions A The momentum distribution for the spin components n ↑(↓) (k) can be evaluated using the corresponding elements of the 8 × 8 Green's function matrix: Very similarly, one can derive a set of two gap equations 67 (k, ıω m ), which at a zero temperature assume the form: Our next step is to solve the gap Equations (5) assuming t = −0.15t (the same value has been used previously in [21]). At x = 1/2, i.e., 3t + m AB = 3 √ 3λ, the band structure is at a topological phase transition, because the Dirac bands at points K = (2π/(3a), 2π/(3 √ 3a) or K = (2π/(3a), −2π/(3 √ 3a)) being massless (at x = 1/2, we have ω 4 (k K ) = 0 and ω 2 (k K ) = 0. Let us assume that ∆ B = α∆ A , so we use Equations (5) to obtain an equation, α = f (α, ∆ A ), for α and Next, we fixed the value of ∆ A , and solve iteratively our equation for α. Having α and ∆ A , we used Equation (5) to obtain the corresponding value of U. The results of our numerical calculations for x = 1/2 are presented in Figure 3. The minimum value of the Hubbard interaction that creates non-zero superfluid gaps is almost the same with and without the NNN hopping term, as can be seen. The difference between the two gaps becomes important for higher values of U. The fact that the minimum value U c ≈ 2.6 that can produce non-zero gap (or a Cooperon bound state, particularly in the vicinity of the Gamma point at which the dispersion has its minimum) has been observed previously in Ref. [22]. Our numerical results that U c does not change significantly with and without NNN hopping term tells us that on-site attractive interaction, rather than the NNN term, boosts the formation of the a Cooperon bound state.

Discussion
To summarize, we use the mean-field approximation to numerically calculate the phase diagram of the attractive KMH model with NNN hoping at half filling, which interpolates the trivial and the non-trivial topological states. It is shown that, as soon as the NNN hoping is included, we have to solve two mean-field gap equations for ∆ A and ∆ B , instead of a single gap equation for ∆ in the case when NNN hopping is neglected.
Recently, the possibility to have superconducting instability in the attractive KMH model has been analyzed within the T-matrix approximation [22]. The question that naturally arises here is about the contributions due to the bubble diagrams, which are included in the Bethe-Salpeter (BS) equation, but neglected by the T-matrix approximation. In Ref. [23], we apply the BS formalism to calculate the slope of the Goldstone mode and the corresponding sound velocity. We found 8% difference between the values of the sound velocity provided by the T-matrix approximation and the BS equation.
It is known that the Gaussian approximation also neglects the bubble diagrams, but, in a square lattice, the difference between the speeds of the sound calculated in the Gaussian and in the BS approximations is about 25% (see Figure 10 in Ref. [24]). To explain the difference of 8% in our numerical calculations, we refer to the system parameters: x = 1, λ = 0.1t, U = 2.69t, and ∆ = 0.151t. From the value of λ follows that E g ≈ 0.52t, and, therefore, the phase diagram is very close to that presented in Figure 2. From another point of view, the value of U = 2.69t tells us that the system is very close to the topological phase transition line at x = 1. The fact that close to the phase transition boundary the speed of sound calculated with in the Gaussian and the BS approaches is essentially the same as has been previously found in a square lattice [25]. Thus, it is natural to expect that away from the phase transition boundary contributions due to the bubble diagrams will be more important.