The doped Mott insulators^[1–6] have been subject to extensive investigations because they contain plentiful phases, such as spin-density wave (SDW), pseudogaps, and superconductivity. In parent compounds, the charge insulation of most Mott insulators raises solely from the Coulomb interaction. However, recently a new family of Mott insulator, in which the charge insulation is resulted from the cooperative interplay between spin–orbit coupling (SOC) and electron correlation, has been proposed, such as Sr₂IrO₄,^[7,8] Na₂IrO₃,^[9–11] and α-RuCl₃.^[12–14] In these Mott insulators, the electron configuration for the active ions is d⁵, such as 5d⁵ of Ir⁴⁺ in iridates and 4d⁵ of Ru³⁺ in ruthenium. Within the octahedral crystal field, the d orbitals are split into a two-fold e_g manifold and a three-fold t_2g manifold with an effective l_eff = 1 orbital moment. The latter are further divided into a fully filled j_eff = 3/2 band and a half-filed j_eff = 1/2 narrow band by strong SOC. Then the j_eff = 1/2 band only needs a moderate electronic correlation to open a Mott gap. Thus these Mott insulators are usually referred to as spin–orbit assisted Mott insulators.

The novel physics of these j_eff = 1/2 Mott materials in different lattices has been reported both experimentally and theoretically. For the square lattice, a canted antiferromagnetic structure has been found in Sr₂IrO₄,^[15–17] which is isostructural to La₂CuO₄, a parent compound for cuprate superconductors. Similar to the doped cuprate Mott insulator, pseudogap, Fermi arc,^[18–21] and superconductivity^[22–26] are also addressed in doped Sr₂IrO₄. For the honeycomb lattice, an intriguing physics is that the Kitaev interactions underlying the celebrated Kitaev model,^[27] whose ground state is exactly a Z₂ spin liquid with Majorana fermion excitations, can be realized.^[28] In this aspect, the honeycomb-lattice material α-RuCl₃ has recently attracted much attention.^{[12–14,29,30]} Although the ground state of α-RuCl₃ has a long-range zigzag magnetic order, a broad continuous spectrum in inelastic neutron scattering has been observed^[12–14] and the magnetic order can be fully suppressed by either an in-plane magnetic field^[31–33] or pressure.^[34–36] There is accumulating evidence that the field-induced disordered state is a promising candidate as a quantum spin liquid.^[31–33] For doped systems, the j_eff = 1/2 Mott materials with honeycomb lattice also exhibit many exotic physics, such as unconventional superconductivity^[37–41] and unique charge dynamics.^[42–49]

In this paper, based on the multi-band Hubbard model including all t_2g orbitals, we investigate various particle–hole (PH) excitations and possible superconducting pairings mediated by these fluctuations in the doped systems with the random-phase-approximation (RPA) method. We find that the j_eff = 1/2 pseudospin fluctuations dominate over the fluctuations from other particle–hole excitations, suggesting the robustness of the j_eff = 1/2 picture even in the doped systems, and the most favorable superconducting state mediated by the pseudospin fluctuations has a d-wave pairing symmetry.

The low-energy physics of α-RuCl₃ is dominated by the t_2g orbitals due to the 4d⁵ configuration of Ru³⁺ and the octahedral crystal field, so our calculations are based on the t_2g three-orbital model,^[30] which is derived through fitting the band structure from the density-functional theory (DFT). The three-orbital Hamiltonian is given as 1 where the parameters in the three-orbital tight-binding Hamiltonian and SOC are from Ref. [30]. The band structure of this three-orbital model with SOC is shown in Fig. 1. The interactions between electrons are included in H_int as follows: 2 where U ( ) is the intra-orbital (inter-orbital) Coulomb interaction, and are the Hund’s coupling and pair hopping, respectively. In this paper, we employ and .

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (a) Brillouin zone of the honeycomb lattice. (b) Band structure of the t2g three-orbital model with SOC for α-RuCl3. The horizontal red dashed line indicates the Fermi level for the undoped system.

Based on the scenario that the pairing interaction arises from the PH fluctuations, we can calculate the effective electron–electron interactions using the RPA. However, as a result of the effect of the SOC, the effective particle–particle interaction cannot be simply divided into bubble and ladder diagrams as is the case with SU(2) symmetry.^[50] Thus, to include the SOC in the RPA approach, we use the Hugenholtz diagrams in our calculations.^[20] In this method, the susceptibilities from the PH excitations are given by 3 where is the Hugenholtz bare vertex^[20] and is the identity matrix. The non-interacting susceptibility is given by 4 with the number of unit cell N and the temperature T. The bare Green’s function is written as 5 Here, the labels i, j, μ, and ν involve all of the degrees of freedom (sublattice, orbital, and spin) in an unit cell. Thus, for the three-orbital model on the honeycomb lattice, is a 12 × 12 matrix, while , , and are 144 × 144 matrices. In the above, ( ) with ( ) for fermions (bosons).

The effective interaction in the particle–particle channel is given by^[20] 6 The gap function is derived from the linearized superconducting gap equation, which is given by 7 Considering that the dominant scatterings occur near the Fermi surface (FS), we can project the effective pairing interaction (Eq. (6)) and the gap equation (Eq. (7)) on the FS. In this way, the scattering amplitude of a Cooper pair from the state on the FS to the state on the FS is calculated by^[51–53] 8 where is the -th eigenvector of the matrix. We then solve the following linearized gap equation: 9 where is the Fermi velocity for band and represents the normalized gap function along the FS . The integration and summation in Eq. (9) are along various FS patches. Here, the temperature is set at T=0.001 eV and the calculation of the susceptibility is done with uniform 64 × 64 meshes in the momentum space. This equation can be solved as an eigenvalue problem, where the eigenvector corresponding to the largest eigenvalue determines the favourable pairing symmetry.

Let us first study the main properties of the PH excitations in the hole doped system. To find the dominant PH excitations, we classify various PH channels according to the total angular momentum of a PH pair. Generally, the susceptibility can be written in terms of i.e., . To show the labels of the effective angular momentums explicitly, we can express the annihilation operator of electrons as , where j and m indicate the quantum numbers of the total effective angular momentum and its z component, respectively, ζ contains other quantum numbers except the angular momentum. Then, the PH excitation operators are written as , which are tensor operators. Since the representations of and are in the sector of effective angular momentum, the operator has the following decomposition in the sector of effective angular momentum: 10 Thus, the operator could be divided into irreducible tensor operators as given in Appendix A. Here, j₁ and j₂ are the angular momenta of the particle and hole. Therefore, including the sublattice indices on the honeycomb lattice, the susceptibilities χ are split into 144 × 144 channels which include j_eff = 1/2 excitations, j_eff = 3/2 excitations, and various spin–orbit excitations between the j_eff = 1/2 and j_eff = 3/2 states. Among these channels, we find numerically that the j_eff = 1/2 pseudospin fluctuations are dominant in the doping regime that we consider in the following discussions. The j_eff = 1/2 pseudospin susceptibilities are defined as 11 where and are the sublattice indices, and the pseudospin operators are given by 12 When the interactions are increased, the magnitudes of pseudospin fluctuations are drastically enhanced in contrast to other types of fluctuations.

In Figs. 2(a)–2(c), we show the diagonal ( , , ) and off-diagonal ( ) channels of the pseudospin susceptibilities together with the largest eigenvalue of along the –M–K– path for the number of electrons in primitive cell n_e=9.8,9.7, and 9.6 (in the parent compound, n_e=10). It can be seen that the pseudospin J_eff = 1/2 susceptibilities capture nearly all features of the maximal eigenvalue of and also show comparable intensities. Thus, the j_eff = 1/2 picture is robust in the doping range we consider here. The differences of , , and shown in Figs. 2(a)–2(c) are the reflections of the asymmetry of the Kitaev exchange interaction in the weak-coupling limit. Moreover, when rotating 120° clockwise along k_c-axis, the transformation of pseudospin fluctuations has a cycle rule: , namely, , and , which is also coincident with the symmetry of the Kitaev interactions.^[27] In the undoped α-RuCl₃, the material shows a zigzag magnetic order, which gives a peak of the susceptibilities at M point in the Brillouin zone.^[14] Our results of the doped systems also show obvious peaks around the M point, which result from the zigzag magnetic fluctuations, though the static zigzag order for the undoped system (n_e=10) can not be obtained in our calculations due to the divergence of the RPA susceptibility. For n_e=9.8 (Fig. 2(a)), the peaks of the susceptibilities deviate from the M point. When the electron concentration is decreased, the peaks move to the M point again, meanwhile the intensities of the peaks are enhanced. In Figs. 2(d)–2(f), we present the pseudospin susceptibility , from which we can see more clearly the distribution of the pseudospin susceptibilities in the whole Brillouin zone.

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. Static susceptibilities and superconducting gap functions for different electron concentrations: (a), (d), (g) ne=9.8, (b), (e), (h) ne=9.7, and (c), (f), (i) ne=9.6. The up panel shows the jeff = 1/2 pseudospin susceptibilities (blue), (green), (red), (magenta), and the largest eigenvalue of defined in Eq. (3) (black) along the –M–K– path. Middle panel: the distribution of the sum of the pseudospin susceptibility in the first Brillouin zone. Panels (g), (h), and (i) are the pairing functions for ne=9.8, ne=9.7, and ne=9.6, respectively.

We then study the properties of superconducting pairing by solving Eq. (9). We find that the maximal values of are in the pseudospin-singlet pairing channel and are far greater than those of the pseduospin-triplet pairing states. So, we will only discuss the superconducting pairing in the singlet channel. In Figs. 2(g)–2(i), we present the real parts of the leading pairing gap functions for three different electron concentrations, which show basically a d-wave symmetry. To understand the detail forms of the pairing gap functions, we try to fit the gap functions by including the electron pairings up to the third nearest-neighbor sites in real space. The fitting functions are given by 13 with being the first, second, and third nearest-neighbor pairing functions, respectively, and , , and being the corresponding pairing intensities. Here, n denotes the electron concentration. Using the adjustable parameters , , and , we find the following optimal fitting functions for the three electron concentrations: 14

Thus, the pairing gap functions for n_e=9.7 and n_e=9.6 have a dominant d-wave superconducting pairing along the nearest-neighbor sites, which can also been seen from the sign and nodes of the gap functions shown in Figs. 2(h) and 2(i). However, the gap function for n_e=9.8 has a dominant d-wave superconducting pairing along the third-nearest-neighbor sites, so additional gap nodes appear in comparison with those for n_e=9.7 and n_e=9.6.

Next, we turn to the electron doped systems. In Fig. 3, we show the pseudospin susceptibility and superconducting gap function for n_e=10.1. From Fig. 3(a), we can find that the peaks of the susceptibility are also located near M and its symmetric points. However, the intensities of the peaks are much larger than those for the hole-doped systems shown in Figs. 2(d)–2(f), which might be the origin of the persistence of the charge gap in the electron-doped systems found in the experiment,^[44] i.e., the strong pseudospin fluctuations hinder the movements of the doped electrons. The superconducting gap function in Fig. 3(b) also shows a d-wave symmetry and the optimal fitting function is given by 15 We find that the electron pairing is dominated by both the nearest-neighbor and next-nearest-neighbor parings, and they have nearly the same intensity.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (a) Static susceptibility and (b) superconducting gap function for electron doped system with ne=10.1.

Thus, we find that the d-wave pairing mediated by the j_eff = 1/2 pseudospin fluctuations is the most favorable superconducting pairing symmetry for both hole and electron dopings. The explicit forms of the d-wave gap functions we obtained here are expected to be helpful for identifying the pairing symmetry if superconductivity in α-RuCl₃ is realized experimentally in the future.

In summary, based on the multi-band Hubbard model with all t_2g orbitals and by use of the random-phase approximation, we have investigated various particle–hole excitations and possible superconducting pairing symmetry in doped α-RuCl₃, in which the strong spin–orbit coupling and electron correlations cause an effective total angular momentum j_eff = 1/2 Mott insulator. We find that the j_eff = 1/2 pseudospin fluctuation dominates over other fluctuations in a large doping range, which suggests that the j_eff = 1/2 picture is still robust in the doped systems. Through projecting particle–particle interaction onto Fermi surfaces, we find that the d-wave pairing is the most favorable superconducting pairing symmetry in both the hole- and electron-doped systems.