Mixture of the nearest- and next-nearest-neighbor d + id-wave pairings on the honeycomb lattice

Symmetry of superconducting pairing is one of the cardinal issues in superconductivity. In general, the nearest-neighbor and next-nearest-neighbor pairings cannot spontaneously mix with each other constrained by the distinct symmetry classification, for example the dx2−y2 - and d xy -wave pairings in the square lattice. The honeycomb lattice provides a special platform to study the mixture of the nearest-neighbor and the next-nearest-neighbor pairing components due to its unique crystal structure. Here, using the variational quantum Monte Carlo simulations, we report a robust superconducting ground state with the intrinsic mixture of the nearest-neighbor and next-nearest-neighbor chiral d + id-wave pairings on the honeycomb lattice. We numerically identify that the two pairing channels promote each other. The proportion of the next-nearest-neighbor d + id component can be readily tuned by doping density and the short-range interactions. Our results provide the efficient ways to manipulate the real-space pairings.


Introduction
Superconductivity is one of the most fascinating issues in condensed matter physics. The topologically nontrivial superconducting state [1,2], such as the spin-triplet p x ± ip y pairing and the spin-singlet d x 2 −y 2 + id xy -wave (abbreviated as p + ip and d + id hereafter, respectively), are of particular interest due to its promising applications in the quantum topological computation [3,4]. The graphene-like materials with the honeycomb lattice provide an ideal platform to realize the chiral superconductivity owing to the unique lattice symmetry and electronic structures.
In principle, if the Cooper pairs are not bound together by phonon exchange, the Coulomb interactions between electrons usually give rise to high-angular momentum Cooper pairing with a sign change in k-space to avoid the strong on-site repulsion. For example, the d x 2 −y 2 -wave pairing in the cuprates [5] and the potentially chiral p + ip-wave pairing in SrRuO 3 [6,7] in the square lattice. On the honeycomb lattice, the d x 2 −y 2 -wave, together with the d xy -wave, belongs to the same two-dimensional E 2g irreducible representation of the D 6h symmetry group [8], so the chiral d + id-wave superconductivity is generally expected. Assuming the phenomenological nearest-neighboring attractive interactions, the d + id-wave superconducting state was first proposed in the graphene at and near the half-filling in the mean-field level [9,10]. Beyond the mean-field approximation, within the standard Hubbard model, the functional-renormalization-group (FRG) [11] and the random-phase-approximation [12] studies showed that the d + id-wave pairing can dominate a certain doping range in the weak coupling limit, while the variational Monte Carlo [13] and quantum Monte Carlo simulations [14] predicted a dominating d + id instability close to the Mott insulator at half-filling in the range from the intermediate to strong coupling. Including the short-range interactions, similar results were also proposed either by the renormalization group scheme [15] or quantum Monte Carlo simulations [16]. The d + id pairing symmetry on the honeycomb lattice was also supported in the researches on the t − J-like models by using the FRG [17], renormalized mean-field [14], and tensor product [18] approaches. Near the van Hove singularity (VHS) filling, most studies suggested that the leading pairing instability is also the chiral d + id-wave pairing [19][20][21][22][23][24], though it is challenged by other competitive orders, such as the chiral spin-density-wave [21,22,24] or uniaxial spin-density-wave [25], as well as the triplet pairing states [23,26]. Besides, in the spin-singlet channel, a cos(3θ)-type pairing function in the B 1g representation of the D 6h group was suggested based on a realistic dispersion fitted from the angle-resolved photoemission spectroscopy data of graphene, where an extended van Hove singularity was observed [27].
On the other hand, the spin-triplet pairings can also naturally avoid the strong on-site Coulomb repulsion. A FRG analysis of the Hubbard model showed that a p + ip-wave pairing state may be favorable close to the half-filling [11]. This was further supported by the variational cluster approximation method, in which a nonchiral p-wave and a chiral p + ip-wave superconductivity was found in the weak and strong coupling regime, respectively [28]. The numerical studies on the Hubbard model in the extremely strong-coupling limit based on the Grassmann tensor product state approach, as well as the quantum field theory, suggested a chiral p + ip-wave superconductivity stemmed from the spin-charge separation and the charge current-current coupling [29]. Moreover, the intra-sublattice triplet f -wave instability was revealed in the extended Hubbard model [12,15,26,30]. Hence, the possible superconducting pairing symmetries of the interacting electron systems on the honeycomb lattice are still under debate.
In addition, the theoretical studies on the honeycomb lattice are mainly focussed on the nearest-neighboring pairing, which gives rise to a inter-sublattice superconducting pairing. However, in principle, the inter-and intra-sublattice pairings should coexist on the honeycomb lattice, since the energy band at any momentum point is an equal proportional mixture of the two sublattices A and B. Thus, in order to correctly describe the superconducting pairing on the honeycomb lattice, at least the nearest-neighbor (NN) pairing of the inter-sublattice component and the next-nearest-neighbor (NNN) pairing of the intra-sublattice component should be included equally. This is in sharp contrast to the situation on the square lattice, where the nearest-neighbor and the next-nearest-neighbor pairing belong to the different irreducible representation [8], the two pairings cannot coexist unless the rotational symmetry is broken. On the honeycomb lattice, both the nearest-neighbor and the next-nearest-neighbor d + id pairing belong to the same E 2g representation, and thus can naturally mix with each other.
In this paper, we study the superconducting ground state on the honeycomb lattice in the intermediate doping range. Performing the variational quantum Monte Carlo simulations, we explicitly demonstrate the positive coupling between the inter-sublattice and intra-sublattice chiral d + id-wave pairing, and thus the mixed d + id-wave pairing is the leading superconducting instability. We also show that the enlarged Fermi surface from around K point towards M point by doping, as well as the enhanced charge fluctuations contributed from the nearest-neighbor Coulomb repulsion, profits the next nearest-neighbor d + id component more. This may provide an opportunity to control the superconducting component on the honeycomb lattice.

Model and methods
The model we adopted is the t-U-V Hamiltonian defined on the hexagonal lattice in real-space as where (i, j) denotes the bonds connecting the sites i and j up to the third NN and i, j denotes the NN bond. c iσ is the electron annihilation operator at the ith site with spin index σ, and n i = σ n iσ is the electron number operator with n iσ the corresponding electron number operator for the given spin flavour. The electron-electron interactions considered here include the on-site and the NN Coulomb repulsion. Such model has already been studied by the quantum Monte Carlo near the half-filling [16] and the functional renormalization group method near the van Hove filling [21], and a superconductivity instability towards the NN d + id-wave was proposed. In this work, we fixed the parameters t = 1, t = 0.1/2.8, and t = 0.07/2.8 for the first, second, and third NN hopping as before [31]. The on-site Coulomb repulsion is set as U = 3.6t, in the intermediate coupling regime. In the pristine graphene, the Coulomb interaction U/t is estimated to be about 2 ∼ 3 [31], but may enhance to 3.3 according to the first principles simulations within the constraint random phase approximation [32]. The honeycomb structured materials can also host the strong correlations, for example In 3 Cu 2 VO 3 [33,34]. We adopt the variational Monte Carlo (VMC) method to study the superconducting instabilities on honeycomb lattice. The trial wave function, chosen to minimize the energy of the Hamiltonian, can be written in general as |Ψ VMC = P|Ψ MF with |Ψ MF the mean-field wave function and P the many-body projecting operation. Any quadratic Hamiltonian can be used to construct |Ψ MF . Moreover, with the appropriate choice of P, the local strong correlation effect such as the Mott physics can also be well simulated. These together make the VMC method not only a flexible but also a reliable method to answer whether the ground state is ordered or not and which kind of order it would develop. We will not present the details of the VMC method but refer to, for example, reference [35][36][37][38] for interested readers.
Here we consider where the operators P N and P S z =0 project the wavefunction into a subspace with the fixed number of particles N and S z total = 0. P J is the charge Jastrow factor where v i, j s are the variational parameters. The on-site v ij projections reproduce the well-known Gutzwiller projection for the Hubbard-like interaction [38], while the non-local v ij projections associate with the charge fluctuations. To well account for the Mott physics [39,40] near the half-filling, we further include the doublon and holon correlations projections P ex d−h [41][42][43][44][45] as where α 0 with = 1, 2 for the NN and NNN are the variational parameters, and ξ i(0) is given by Here τ ∈ runs over all the th NN sites, and D i = n i↑ n i↓ and H i = (1 − n i↑ )(1 − n i↓ ) are doublon and holon operators, respectively. Figure 1(c) illustrates that the holons and doublons tends to form the bonds between the NN or the NNN to avoid the energy penalty. The typical energy reduction per site after P ex projection in the present model is about 5.2 × 10 −3 and 2.2 × 10 −4 for the NN and the NNN doublon-holon pairs at 0.18 electron-doping on 2 × 12 × 12 lattice, values comparable, and even larger than the energy reduction contributed from the superconducting orders as shown below. In principle, the higher angular momentum pairings, containing more nodes, are usually energetically unfavourable [46]. We, therefore, study four most possible and controversial superconducting pairing channels in real-space, namely, the extended-s-wave(es-wave), p + ip-wave, d + id-wave and f -wave. These pairing channels belongs to four different irreducible representations as shown in table 1. The singlet and triplet pairing operators can be written as: The mean-field Hamiltonian we used for |Ψ MF is

Symmetry
Operators Irreps Here i, j denotes the NNN bond. The real-space Bogoliubov-de Gennes Hamiltonian and particle-hole transformation are used to construct |Ψ MF . The chemical potential μ and superconducting gap Δ in the mean-field Hamiltonian, as well as the parameters v ij and α l 0 in projection operation, are used as the variational parameters to minimize the ground-state energy of the t-U-V model. The finite variational parameter for the superconducting gap in the mean-field Hamiltonian optimized by the VMC method is not enough to guarantee the emergence of long-range superconducting order. We further evaluate the superconducting correlation function [13] as: where Δ l (r 0 ) and Δ l (r) are the superconducting pairing operators for the l-pairing channel at r 0 and r, respectively. The long-range superconducting order would establish in the case of the finite value of χ ll (r − r 0 ) in the long-distance limit. The VMC simulations are performed on the 2 × 12 × 12 lattice with the periodic boundary conditions to ensure that both the high symmetric K and M points, which are crucial for the fillings close to the half-filling and the van Hove filling, respectively, are included in the reciprocal space. The typical number of Monte Carlo samples is at least 2 27 and the error is about 10 −5 unless otherwise specified. To ensure the robustness of the superconducting order. We mainly focus on the intermediate doping range below the van Hove filling.

Single-component superconductivity
To begin, we study four single-component pairing channels, i.e., the es-, p + ip-, d + id-, and f -wave pairing, respectively, to find out the leading superconducting instability on the honeycomb lattice. As shown in figure 2 at the fixed electron-doping δ ≈ 0.18 (340 electrons in 2 × 12 × 12 unit cell), the chiral NN-d + id-wave pairing is the dominant instability among all the considered superconducting channels, while the triplet NN-p + ip and NN-f -wave pairing are the subdominant instability with a reduced energy reduction about ∼55% and ∼39% of that of the NN-d + id-wave pairing relative to the normal state. We find no instability towards the es-pairing at the present setting, in comparison with the subdominant tendency in quantum Monte Carlo simulations with the relative small U/t [16], but in agreement with the previous variational cluster approximation result, though a dominant p-wave pairing instability was proposed in the latter [28]. The dominating NN-d + id pairing close to the Mott insulator at intermediate doping is consistent with previous numerical studies [13,14,18]. On the other hand, the NNN-d + id pairing exhibits a subdominant instability with a relative small energy reduction (about 40%) compared with that of the NN-d + id pairing. The strong instability tendency towards the NNN-d + id pairing was recently reported at the van Hove filling in the same model by the determinantal quantum Monte Carlo calculations, in which no triplet pairing was observed in larger sized systems [24]. Previous renormalization group method also indicated a subdominant NNN-d + id-wave pairing instability at finite V [21]. We An interesting question then follows: whether these pairing instabilities can simultaneously emerge, i.e., a mixed superconducting state with multiple components. In principle, the pairings belonging to the different irreducible representation of the corresponding crystal point group cannot mix with each other unless the additional symmetry breaking. According to the classification of the irreducible representation of D 6h of the honeycomb lattice [8], no triplet and singlet mixture can be achieved. However, the NN and NNN d + id-wave pairing belong to the same irreducible representation E 2g , and therefore they can mix with each other at least from the view of group theory. It is reminiscent of the situation in the electron-doped cuprate, in which a coexistence of NN and third-NN d x 2 −y 2 -wave superconductivity was experimentally proposed to account for the non-monotonic cos 6θ-component deviating from the standard cos 2θ-component [47][48][49] (remaining controversial). Here, the honeycomb lattice provides an additional opportunity to realize the multi-component superconductivity-a mixture of the NN and NNN d + id-wave pairing. We now turn to study the potentially mixed d + id superconductivity on the honeycomb lattice.

Mixed d + id superconductivity
To explore the coexistence of the NN-and NNN-d + id-wave components, we perform the VMC calculations in the Δ NN d+id -Δ NNN d+id parameter space at electron doping ∼0.172 in a relative small size lattice of 2 × 8 × 8. As shown in figure 3(a), the resultant heatmap clearly illustrates the existence of a net energy reduction when two superconducting orders are simultaneously included. We have checked this reduction in larger lattice of 2 × 12 × 12 to rule out the size effect. In figure 3(b), we show the variational process of the respective superconducting component at the optimal condition of the mixed state on the 2 × 12 × 12 lattice, both NN-and NNN-d + id-wave superconductivity exhibit an energy minimum with distinct energy reduction. In particular, the energy reduction in the mixed state of respective d + id channel is even larger than the corresponding counterpart in the single superconducting state. This strongly suggests the positive coupling between the NN-and NNN-d + id-wave pairings, and further confirms the emergence of the mixed superconducting state. We further check the pairing correlation χ(r − r 0 ) as a function of the distance for both pairing channels as shown in figure 3(c). The pairing correlations converge to a finite value at the long distant limit, manifesting the long-range order nature of the NN-and NNN-d + id pairing for the present choice of the model parameters. The long-range order of the NN-d + id-wave pairing agrees with the conclusion in a similar VMC study [13], but conflicts with the result obtained by the previous determinant Monte Carlo simulations at relative lower doping density [16]. This enhanced tendency of the superconducting pairing may be due to the increasing density of states towards the von Hove filling. Moreover, the superconducting correlations for both channels in the mixed state are stronger than their counterparts in the single superconducting state, supporting the existence of the mixed superconducting state further. We also show the correlations between the NN and NNN-d + id components, the non-vanishing values at the long distant limit also support the coexistence of those two pairing channels, as well as their long-range order nature. Hence, the mixture of NN and NNN d + id-wave superconductivity is the ground state of the honeycomb lattice in the intermediate doping.
We further explore the degree of mixture of the NN-and NNN-d + id-wave superconductivity by defining a ratio of γ = Δ NNN d+id /Δ NN d+id . Due to the increase of the density of state towards the von Hove filling, the magnitudes of both components are expected to be enhanced. Our numerical data well support this idea (figure 4(a)), and well agree with the previous VMC data [13,21], as well as the renormalized mean-field results [14,17]. In particular, the ratio γ increases quickly upon the electron doping density δ, suggesting the NNN d + id-wave pairing benefits more from the doping. As well known, the Fermi surface on the honeycomb lattice evolves from the small pocket around the K point towards the large one around the Γ point when the doping density increases (below 1/4 doping), and, therefore, the near-Fermi-surface contributions from the K-region decreases, whereas the near M-region increases, respectively (low panels in figure 4). On the other hand, the order parameter of the NN-d + id-wave component dominates the K-region (figure 4(c)), while the NNN-d + id-wave component dominates the M-region ( figure 4(d)). With the increasing doping density, more NNN-d + id-wave proportions benefit larger gap opening around the Fermi surface, and are expected to lower the system energy. It should be noted that the pairing instabilities of d + id-wave become stronger with the increasing doping, however, the magnetic instability also increases due to the stronger Fermi surface nesting. In particular, the perfect Fermi surface nesting at the von Hove filling prefers the chiral spin-density-wave state as revealed by the various numerical simulations [20,21,24].
We also investigate the role of the long-range interactions on the two-component superconductivity by considering the NN repulsion V. The optimal superconducting variational parameter of the NN d + id-wave pairing weakens quickly upon the growing V, whereas that of the NNN d + id strengthens slightly. In particular, the NNN d + id-wave component surpasses the NN d + id-wave one above V ≈ 0.9t. This means that the long-range Coulomb interaction strongly suppresses the NN d + id-wave pairing, but enhances the NNN d + id-wave pairing. In general, the electrostatic screening is relatively weak in the strongly correlated materials due to their bad-metal nature. The virtual electron-hopping processes in doping case may also enhance the value of V in comparison with their undoped parents [50]. Moreover, V is material-dependent and can be tuned by pressure. In this sense, the long-range Coulomb interactions provide a way to manipulate the superconducting components. Similar single component NN d + id suppression via V was also observed by the determinantal quantum Monte Carlo simulations [16], whereas the single component NNN d + id enhancement was predicted by the renormalization group method but with a much smaller amplitude compared with its NN counterpart [21].
It is also interesting to investigate the NNN f -wave instability, which was proposed to dominate for the longer-range Coulomb interactions [12,15,26,30]. Our VMC data show that the NNN f -wave single component superconductivity enhances upon V, in comparison with the decreasing NN d + id-wave one ( figure 5(a)). However, as shown in figure 5(b), the single NNN d + id-wave component also enhances with V, and dominates over the NNN f -wave one in a wide range V. As discussed before, the mixed d + id-wave state has lower energy than the single d + id-wave component. Therefore, our results suggest that the mixed state is still the ground state even for strong V.

Conclusion and discussion
In summary, we employ the variational quantum Monte Carlo method to study the possible superconductivities on the honeycomb lattice based on the doped extended Hubbard model. We find that the most favorable superconductivity occurs in the singlet channel with a E 2g pairing symmetry. We explicitly demonstrate that there is a positive coupling between the inter-sublattice and intra-sublattice chiral d + id-wave pairings, and thus the mixed d + id-wave pairing with both the nearest-neighbor and next-nearest-neighbor components is the leading superconducting instability. We also show that the ratio of the nearest-neighbor and next-nearest-neighbor components can be controlled by the doping and interactions.
In this paper, we clearly demonstrate that the mixture of the NN-and NNN-d + id-wave pairing is the superconducting ground state in the intermediate doping region. Whether this mixed superconducting state could be the ground state even at the van Hove filling has not been discussed before. In principle, due to the perfect Fermi surface nesting, the chiral spin-density-wave state is believed to be the ground state at the van Hove filling [21,51]. However, the long-range hopping terms may suppress the perfect nesting, and the competition between the mixed superconducting state and the chiral spin-density wave state should be sure of interest. On the other hand, the introducing of the short-range Coulomb interaction may trigger the charge-density-wave, the competition between multiple orders, i.e., the spin-density-wave, charge-density-wave, and superconducting orders, should be also interesting. In addition, since the kagome lattice shares the similarity of the honeycomb lattice. The mixed superconducting state and its competition with other ordered states [52,53] are surely deserved to further study.