Ni-based transition-metal trichalcogenide monolayer: a strongly correlated quadruple-layer graphene
Yuhao Gu, Qiang Zhang, Congcong Le, Yinxiang Li, Tao Xiang, Jiangping Hu
NNi-based transition-metal trichalcogenide monolayer: a strongly correlated quadruple-layergraphene
Yuhao Gu,
1, 2
Qiang Zhang, Congcong Le,
3, 2
Yinxiang Li, Tao Xiang,
2, 3, 4 and Jiangping Hu
2, 3, 4, ∗ Beijing National Laboratory for Molecular Sciences, State Key Laboratory of Rare EarthMaterials Chemistry and Applications, Institute of Theoretical and Computational Chemistry,College of Chemistry and Molecular Engineering, Peking University, 100871 Beijing, China Beijing National Laboratory for Condensed Matter Physics,and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China Kavli Institute of Theoretical Sciences, University of Chinese Academy of Sciences, Beijing, 100190, China Collaborative Innovation Center of Quantum Matter, Beijing, China
We investigate the electronic physics of layered Ni-based trichalcogenide NiPX (X=S, Se), a member oftransition-metal trichalcogenides (TMTs) with the chemical formula, ABX . These Ni-based TMTs distinguishthemselves from other TMTs as their low energy electronic physics can be effectively described by the two e g d-orbitals. The major band kinematics is characterized by the unusal long-range effective hopping between twothird nearest-neighbor (TNN) Ni sites in the two-dimensional Ni honeycomb lattice so that the Ni lattice canbe equivalently viewed as four weakly coupled honeycomb sublattices. Within each sublattice, the electronicphysics is described by a strongly correlated two-orbital graphene-type model that results in an antiferromag-netic (AFM) ground state near half filling. We show that the low energy physics in a paramagnetic state isdetermined by the eight Dirac cones which locate at K , K (cid:48) , K and K (cid:48) points in the first Brillouin zone witha strong AFM fluctuation between two K ( K (cid:48) ) and K (cid:48) ( K ) Dirac cones and carrier doping can sufficientlysuppress the long-range AFM order and allow other competing orders, such as superconductivity, to emerge.The material can be an ideal system to study many exotic phenomena emerged from strong electron-electroncorrelation, including a potential d ± id superconducting state at high temperature. I. INTRODUCTION
Since the discovery of graphene[1] a decade ago, two-dimensional (2D) materials have been a research frontier forboth fundamental physics and practical device applications[2,3]. Transition-metal trichalcogenides (TMTs) with the chem-ical formula ABX (X=S, Se, Te), which were known morethan a century ago[4, 5], are layered van der Walls (vdW) ma-terials. Recently, this family of materials has attracted greatresearch attention as potential excellent candidates to explore2D magnetism for novel spintronics applications.All the members in the family of ABX materials are builton a common structural unit, (P X ) (X=S, Se, Te) anioncomplex. However, the cation atom A is rather flexible, rang-ing from vanadium to zinc (A=V, Cr, Mn, Fe, Co, Ni, Zn, etc.)in the row of the 3 d transition metal, partial alkaline metalin group-II, and some other metal ions. As shown in Fig.1(a,b), the cation is coordinated with six chalcogen anions toform an octehedra complex. In the two-dimensional layer,the cation forms a graphene-type honeycomb lattice. Thetransition-metal trichalcogenides exhibit a variety of intrigu-ing magnetically ordered insulating states[6]. Recently, underhigh pressure, FePSe can also become a superconductor[7].Among this family of materials, the Ni-based trichalco-genides can carry intriguing electronic physics, such as strongcharge-spin coupling[8], because of the following reasons.First, as the transition metal cation and chalcogen anions forman octahedral complex, the 3 d -orbitals of the transition metalare divided into high energy e g and low energy t g groups. ∗ [email protected] In the case for Ni which has 8 electrons in the 3 d -shell, the t g orbitals are fully occupied and the two e g orbitals arehalf-filled as shown in Fig.1 (c). The t g orbitals are inac-tive. The Ni-based trichalcogenides should be described by arelatively simpler low energy effective model than other ma-terials. Second, unlike a two dimensional square lattice, ahoneycomb lattice easily exhibits a Dirac-cone type of energydispersion. Near half filling, both one-orbital model, such asgraphene[1], and two-orbital models[9, 10] in the honeycomblattice are featured with Dirac points near Fermi energy. Withthe strong electron-electron correlation in the 3 d -orbitals, theNi-based trichalcogenide thus can be a candidate of stronglycorrelated Dirac electron systems. It is worth to mentionthat a recent major research effort has aimed to increase theelectron-electron correlation in graphene[11], in which flatbands have to be engineered to observe correlation effects be-cause of p -orbitals. Finally, both density functional theory(DFT) calculation and experimental measurements have sug-gested that the Ni honeycomb lattice forms the zigzag antifer-romagnetic insulating ground state featured as double paral-lel ferromagnetic chains being anti-ferromagnetically (AFM)coupled[6, 12]. The material offers a promising platform tostudy the interplay between the low energy Dirac electronicphysics and the magnetism. Such an interplay is believed tobe responsible for many important phenomena, for example,high temperature superconductivity in both cuprates and iron-based superconductors[13, 14].In this paper, we show that the Ni-based TMTs are Diracmaterials with strong electron-electron correlation. Their lowenergy electronic physics can be entirely attributed to the two e g d -orbitals with a band kinematics dominated by an unusual”long-range” hoppings between two third nearest-neighbor(TNN) Ni sites in the Ni honeycomb lattice. Thus, the original a r X i v : . [ c ond - m a t . s t r- e l ] N ov Ni lattice can be divided into four weakly coupled honeycombsublattices. Within each sublattice, the electronic physics isdescribed by a strongly correlated two-orbital graphene-typemodel. The couplings between four sublattices, namely, thenearest-neighbor (NN) and the second NN hoppings (SNN)in the original lattice, can be adjusted by applying externalpressure or chemical methods. In the absence of the strongelectron-electron correlation, the low energy physics is deter-mined by the eight Dirac cones which locate at K ( k (cid:48) ) andthree non-equivalent pairs of K ( K (cid:48) ) points in the first Bril-louin zone. In the presence of strong electron-electron corre-lation, strong AFM interactions arise between two NN siteswithin each honeycomb sublattice. Namely, in the view of theoriginal honeycomb lattice, the strong AFM interactions onlyexist between two TNN sites. Near half filling, Dirac conesare gapped out by the long-range AFM order. Using the stan-dard slave-boson approach, we show that the doping can suffi-ciently suppress the long range AFM order. In a wide range ofdoping, a strong AFM fluctuation can exist between the Diraccones and a d ± id superconducting state can be developed.It is interesting to make an analogy between the above re-sults and those known in high temperature superconductors,cuprates and iron-based superconductors. For the latter, itis known that the dominant AFM interactions are betweentwo NN sites in curpates and between two SNN sites in iron-based superconductors, which are believed to be responsiblefor the d -wave and the extended s -wave pairing superconduct-ing states respectively[13–15]. All these AFM interactions aregenerated through superexchange mechanism. Thus, the Ni-based TMTs, having the AFM interactions between two TNNsites, can potentially provide the ultimate piece of evidenceto settle superconducting mechanism in unconventional hightemperature superconductors.This paper is organized as follows. In Section II, webriefly review transition-metal phosphorous trichalcogenidesand specify our DFT computation methods for magnetism andband structures. In Section III, we analyze the band structureof the paramagnetic states, derive the tight-binding Hamilto-nian and discuss the low energy physics near Fermi surfaces.In Section IV, we calculate the magnetic states of the ma-terials and derive the effective magnetic exchange couplingparameters. In Section V, we use the slave-boson meanfieldmethod to derive the phase diagram upon doping and discussthe possible superconducting states. In the last section, wemake summary and discussion. II. TRANSITION-METAL PHOSPHOROUSTRICHALCOGENIDES
The MPX metal phosphorous trichalcogenides (M=Mg,Sn, Sc, Mn, Fe, Co, Ni, Cd, etc and X=S, Se) are a famousfamily of 2D van der Waals (vdW) materials[2, 3]. The bulkMPX crystal consists of AA-stacked or ABC-stacked single-layer assemblies which are held together by the vdW interac-tion. The vdW gap distance of the MPX with 3 d transitionmetal elements is about 3.2 ˚A, much wider than the well stud-ied MoS -type 2D vdW materials[16], which indicates that TABLE I. The structural parameters for monolayer MPX (spacegroup P − m ).System a ( ˚A) M-X ( ˚A) P-X ( ˚A) P-P ( ˚A) M-X-M ( ◦ )MnPS [19] 6.08 2.63 2.03 2.19 83.9FePS [19] 5.94 2.55 2.02 2.19 84.7CoPS [20] 5.91 2.51 2.04 2.17 85.5NiPS [21] 5.82 2.50 1.98 2.17 84.4NiPSe [22] 6.13 2.61 2.09 2.24 85.2 P S/SeM g e g t d (a) (c)(b) xy M FIG. 1. (a) The crystal structures of the monolayer NiPX (X=S, Se)(space group P-31m). (b) The top view of the monolayer NiPX . (c)The octahedral crystal field splitting of Ni atoms. the vdW interaction in MPX is relatively weak. The mono-layer structure of MPX is constructed by MX edge sharedoctahedral complexes. Similar to the MoS -type MX materi-als, the monolayer MPX can be considered as the monolayerMX with one third of M sites substituted by P dimers, i.e.,MPX can be considered as M (P ) X . Thus, the trianglelattice in MX transforms to the honeycomb lattice in MPX with the P X anions being located at the center of the hon-eycomb. The P-P and P-X bond lengths indicate that the P-Pand P-X bonds are covalent bonds in P X anions. As thetransition-metal atoms are in octahedral environment, the five d orbitals are split into two groups, e g and t g , as shown FIG.1(c).With the weak vdW interaction, the essential electronicphysics in MPS is determined within a monolayer. Infact, the atomically thin FePS layer has been experimentallysynthesized[17]. In this paper, without further specification,our study and calculations are performed for the monolayer asshown in Fig.1 (a) and (b). We use the monolayer structurescleaved from the experimental crystal structures in ICSD[18].The experimental structural parameters of the MPX are listedin the Table.I.Our DFT calculation is performed for the monolayer MPX -2-1 0 1 Γ K M Γ E ne r g y ( e V ) Ni d xz/yz
Ni d xy/x -y Ni d z S p -2-1 0 1 Γ K M Γ E ne r g y ( e V ) (a) Ni d xz/yz
Ni d xy/x -y Ni d z Se p E A A (b) FIG. 2. Electronic band structures of (a) NiPS and (b) NiPSe . Theorbital characters of bands are represented by different colors. structures together with built-in 20 ˚A thick vacuum lay-ers. We employ the Vienna ab initio simulation package(VASP) code[23] with the projector augmented wave (PAW)method[24] to perform the DFT calculations. The Perdew-Burke-Ernzerhof (PBE)[25] exchange-correlation functionalwas used in our calculations. Through out this work, the ki-netic energy cutoff (ENCUT) is set to be 500 eV for the ex-panding the wave functions into a plane-wave basis and the Γ -centered k-mesh is × × for the nonmagnetic unit cell.The energy convergence criterion is − eV. We perform thestatic self-consistent calculation with the monolayer structurecleaved from the experimental crystal structures in ICSD[18].Since there is no experimental data about CuPS with this lay-ered hexagonal structure, we optimize CuPS with the forceconvergence criterion of 0.01 eV/ ˚A. In the study of effec-tive Hamiltonian, we employ Wannier90[26] to calculate thehopping parameters of the tight binding model. In the studyof magnetism of MPX , the GGA plus on-site repulsion U method (GGA+ U ) in the formulation of Dudarev et al. [27] isemployed to describe the associated electron-electron corre-lation effect. The effective Hubbard U ( U eff ) is defined by U eff = U − J Hund . In order to describe different magneticorders, we build × × supercell and the k-mesh is × × ,correspondingly. III. ELECTRONIC BAND STRUCTURES AND THETIGHT BINDING MODEL FOR NI-TRICHALCOGENIDES
In Fig.2, we plot the band structures of NiPS and NiPSe ,which are very similar to each other. From Fig.2, it is clearthat the five d -orbital bands are divided into two groups sep-arated by a large crystal field splitting energy. The groups atthe high and low energy are attributed to the two e g and three t g orbitals respectively. The bands from t g orbitals are com-pletely filled while the four bands from the two e g orbitals areclose to half-filling. This is consistent with the fact that theNi cations are six-coordinated with an octahedral geometryand the d configuration in Ni contributes two electrons tothe two e g orbitals. The physics near Fermi energy is con-trolled by the two e g d xz/yz orbitals. It is also worthy notingthat the contribution of the S/Se- p orbitals is considerable, es- TABLE II. The NN, SNN and TNN hopping parameters and the bandwidths for NiPS and NiPSe ; xz and yz represent the band indexesof the hopping parameters. NiPS NiPSe t NNxzxz (eV) -0.050971 -0.059617 t NNyzyz (eV) -0.036294 -0.014879 t SNNxzxz (eV) 0.012118 0.019696 t SNNyzyz (eV) -0.015141 -0.017257 t SNNxzyz (eV) 0.003175 0.003827 t TNNxzxz (eV) -0.020218 -0.019138 t TNNyzyz (eV) 0.238574 0.256818band widths e g with different Ni-S-Ni angles; xz and yz represent the band indexesof the hopping parameters.Ni-S-Ni ( ◦ ) 84.4 81.9 87.1 t NNxzxz (eV) -0.050971 -0.038045 -0.066167 t NNyzyz (eV) -0.036294 -0.026956 -0.040663 t SNNxzxz (eV) 0.012118 0.011367 0.013838 t SNNyzyz (eV) -0.015141 -0.009442 -0.021795 t SNNxzyz (eV) 0.003175 0.001937 0.004137 t TNNxzxz (eV) -0.020218 -0.021682 -0.019388 t TNNyzyz (eV) 0.238574 0.210725 0.269932 pecially near Γ point, which indicates strong d - p hybridiza-tion.Similar to the single orbital band structure in graphene, thetwo e g orbital bands are featured by Dirac points as well.Here, there are eight Dirac points which locate at K and K (cid:48) points, as well as around K/ and K (cid:48) / points along Γ - K line as shown in Fig.2. These Dirac points are rather robust.To show the robustness, we analyze the symmetric characterof the bands, namely, their irreducible representations. Thetwo crossing bands along Γ - K line belong to A and A ir-reducible representations, and the bands at K points belongto E irreducible representation. Thus, the Dirac points areprotected by the symmetry as the bands belong to differentrepresentations.In order to capture the two-dimensional electronic physicsnear the Fermi level, we construct the tight binding Hamilto-nian based on the two e g orbitals. The Hamiltonian can bewritten as H = (cid:88) k ψ † k h k ψ k , (1)where the basis ψ † k = ( a † xk , a † yk , b † xk , b † yk ) and h k = (cid:18) ω k − µ γ k γ † k ω Tk − µ (cid:19) (2) (b)(a) FIG. 3. (a) The three NN, SNN, and TNN hopping parametersmarked by the green, brown and black dashed arrows, respectivelyand the zigzag AFM order with onsite red/blue arrows indicating spinup/down; (b) The four honeycomb sublattices. with µ being the chemical potential and ω k = (cid:88) j e − i k · a T SNNj (3) γ k = (cid:88) j e − i k · a T NNj + e − i k · a T SNNj . (4)Here a † xk ( b † yk ) is the electron annihilator operator of orbital xz ( yz ) in the usual A ( B ) sublattice of the honeycomb latticeand vectors a , a , a are the first, second and third neighborvectors. T ij = C j T i C − j is the i -th neighbor hopping matrixvia the bond along ij bond direction and C j is the threefoldrotation operation to the ij direction relative to the initial set-ting. T i ( i = N N, SN N, T N N ) is the hopping matrix withthe direction marked in Fig.3 (a): T i = (cid:18) t ixzxz t ixzyz − t i ∗ xzyz t iyzyz (cid:19) . (5)By the lattice symmetry, t NNxzyz = t T NNxzyz = 0 . We will useeV as the energy unit for all hopping parameters. The resultsof NiPS and NiPSe are similar, as shown in Table.II. Theexplicit formula of the Hamiltonian is given in Appendix A.Here we focus on the results of NiPS .It is interesting to notice that the leading term in the aboveHamiltonian is t T NNyzyz , the TNN σ -bond hoppings as shown inFig.3 (a), which is almost one order of magnitude larger thanthe other hopping parameters, namely, the NN and SNN hop-ping parameters. Thus, we can consider these TNN hoppingsas the dominant hopping parameters and treat other hoppingsas perturbations. In Fig.4 (a), we plot the band dispersion withonly the TNN hopping parameters. With only these TNN hop-pings, the original Ni honeycomb lattice is divided into fourdecoupled sublattices as shown in Fig.3 (b). Within each hon-eycomb sublattice, the model is identical to the one previouslystudied in an ultracold atomic honeycomb lattice with two de-generate p-orbitals[9, 10]. As shown in Fig.4 (a), there are twocompletely flat bands and two dispersive bands. The flat bandsstem from the localized binding and anti-binding molecularorbitals[9]. The two dispersive bands create the eight Diracpoints. With only these TNN hoppings, the second pair of Γ K M Γ K/2 Γ K M Γ K/2 E ne r g y ( e V ) E ne r g y ( e V ) (a) (b) D D D D D D D D FIG. 4. The effect of different hopping parameters on the band struc-ture and Dirac points. (a) The band dispersion with only the leadingTNN hopping t TNNxzxz ; (b) The band dispersion with all hopping pa-rameters: the green, purple and black arrows represent the motionof the band and Dirac points by increasing the NN, SNN and TNNhopping parameters, respectively. Q Q KM Γ K`/2
FIG. 5. Fermi surfaces and nesting vectors. The electron pocketsand hole pockets are marked by red and green respectively. The in-sets from top to bottom are Fermi surfaces at three doping levels , 0.1(electron), -0.1 (hole) and 0.5 (electron) per Ni atom with respect tothe half filling, corresponding to formula x = ( n − / with n thetotal electrons in each unit cell. In (a), nesting vectors Q and Q are depicted. Dirac points are exactly located at K/ and K (cid:48) / points. Thispair is simply created through the Brillouin zone folding be-cause of the sublattice structure. Thus, the presence of the twopairs of Dirac points underlines the sublattice structure.The dominance of the t T NNyzyz in NiPX can be understoodfrom the lattice chemistry. The Ni- e g orbitals are stronglycoupled with S/Se-p orbitals. These effective hoppings aremediated through the central P X anion. For the NN hop-pings, two NN Ni atoms are in two edge shared MX octa-hedral complexes. As the Ni-X-Ni angle is close to 90 ◦ , theNN indirect hopping through X is very small. The SNN ef-fective hopping is mediated by two S/Se atoms which sepa-rately locate in the top and bottom layers. The coupling be-tween these two S/Se atoms is weak due to the long distancearound 3.8 ˚A between them, which explains the weak SNNhoppings. By contrast, the TNN σ hopping parameter is me-diated through two S/Se atoms in the same layer.The effects of other hopping parameters on Dirac points andband structures are indicated in Fig.4 (b), in which the arrowsrepresent the motion of Dirac points and band structures whenthe corresponding hopping parameters increase. More specifi-cally, the weak third neighbor π -bond hoppings t T NNxzxz neitheraffect the Dirac cones at K and K/ , nor the band degeneracypoints at Γ and M . They only affect the flat bands in Fig.4 (a)far away from Fermi energy. The flat bands turn to dispersewhen t T NNxzxz increases. Therefore, in the weak hopping region,the low energy physics near Fermi surfaces are not affected bythe third neighbor π -bond hoppings. The weak second nearestneighbor hoppings, t SNNxzxz and t SNNyzyz , shift the Dirac points at K/ and K vertically. By increasing these hoppings, the twoDirac points shift in opposite directions by a shift ratio equalto as indicated by the purple arrows in Fig.4 (b). The weakNN hoppings, t NNxzxz and t NNyzyz , donot affect the Dirac cone at K because of the symmetry protection from the C , time re-versal and inversion symmetries. However, it drags the K/ Dirac cone along Γ − K line as indicated by the green arrowsin Fig.4 (b). Band crossing points D and D are also draggedalong the direction indicated by the green arrows in Fig.2 (b).The Fermi surfaces at the different doping levels are shownin Fig.5. Without the SNN hoppings, the model gives Diracsemimetals at half filling. Thus, the tiny pockets at half fillingshown in Fig.5 stem from very small SNN hoppings. Due tocharge conservation, the area of electron pockets at K/ arethree times smaller than those hole pockets at K . In principle,with very small hole doping, strong nesting can take placebetween the electron and hole pockets at K and K (cid:48) / respec-tively, but not at K and K/ , by taking into consideration ofthe shapes of Fermi pockets. By increasing hole (electron)doping, both pockets at K and K/ become hole (electron)pockets. When the doping reaches around 0.3 carriers per Niatoms, there is a Lifshitz transition of Fermi surfaces, namely,the two pockets emerge together to become one Fermi surface.From the above Fermi surface topology, we can considerthe possible Fermi surface nesting in a paramagnetic state.Near half filling, the nesting vector is given by Q = G / ,half of the reciprocal lattice vector, as highlighted in Fig.5.This vector is exactly the ordered magnetic wavevector in theAFM zigzag state. We calculate the spin susceptibility underrandom phase approximation (RPA), with the same methodand notations specified in literature[28, 29]. The result is plot-ted in Fig.6 for several different doping levels. Clearly thesusceptibility peak emerges at M ( Q ) near half filling. Be-low the critical doping at the Lifshitz transition, the peak iswell preserved, indicating the existence of strong AFM fluc-tuations. IV. MAGNETIC EXCHANGE COUPLING PARAMETERSAND THE AFM ZIGZAG STATE
Without doping, MPX are known to be magneticinsulators[2, 3, 12, 30]. As the magnetic moments are local-ized at the transition metal atoms, the magnetism can be cap- χ s RPA χ χ s RPA χ χ s RPA χ χ s RPA χ (a) (b)(c) (d) Γ K MK/2
Γ Γ
K MK/2 ΓΓ K MK/2 ΓΓ K MK/2 Γ FIG. 6. Bare (dashed blue) and RPA (solid red) approximated spinsusceptibility for different doping levels: (a) half filling, (b) 0.1, (c)-0.1 and (d) 0.3 (Liftshiz point where K and K/ pockets connect).Here the onsite energy U = 0 . eV and Hund’s coupling J h = 0 . U is adopted, similar to ref.[28]. Resonance apexes appear around thenesting vector Q (M) and Q (K/2) marked in Fig.5. tured by an effective Heisenberg model with local magneticmoments. As the effect of the spin orbital coupling is gener-ally small for e g orbitals, we expect an isotropic Heisenbergmodel. Furthermore, from the lattice structure, it is obviousthat the minimum effective model should include NN, SNN,and TNN magnetic exchange coupling parameters. Namely,the model can be written as H = J (cid:80)
Ni d xy/ -yz S p (a) (b)
Mn Fe Co Ni Cu-20020406080100120140 U eff =0 eV U eff =1 eV U eff =2 eV U eff =3 eV U eff =4 eV J ( m e V ) FIG. 7. (a) The band structure of NiPS in the AFM zigzag state.(b) J superexchange AFM interactions in MPS (M = Mn, Fe, Co,Ni, Cu), which are extracted from the GGA+ U calculations with thevalues U eff = (0 , , , , eV. values.The classical energies of the four different magnetic statesfor the effective Heisenberg model are given by E F M = S (6 J + 12 J + 6 J ) + E ,E AF M − Neel = S ( − J + 12 J − J ) + E ,E AF M − zigzag = S (2 J − J − J ) + E ,E AF M − stripy = S ( − J − J + 6 J ) + E (7)From the calculated energies of these states, we can extractthe effective magnetic exchange interactions. The results arelisted in Table.V. Some similar results have been obtainedpreviously[12, 30]. Our calculation are consistent with theseprevious calculations[12, 30].Here we pay special attention to the values in NiPX . Asshown in Table.V, for NiPX , among the three magnetic ex-change coupling parameters, J is one order of magnitudelarger than the other two parameters. Moreover, J is stronglyAFM while J and J both are weakly FM. These qualitativefeatures are independent of U eff . The dominance of J overthe other two further confirms the extracted physical pictureof weakly coupled four sublattices as shown in Fig.3 (b) basedon the hopping parameters in the electronic band structure. J stems from so-called AFM super-superexchangeinteraction[6, 30]. In Fig.7 (b), we plot the values of J asa function of M (M=Mn, Fe, Co, Ni, Cu). It is important to TABLE V. The calculated exchange interaction parameters J , J and J for monolayer MPS (M = Mn, Fe, Co, Ni) and NiPSe usingGGA+ U ( U eff = 0 or eV). U eff = 0 J (meV) J (meV) J (meV)MnPS -20.90 5.08 3.81CoPS -22.52 14.50 5.78NiPS -10.63 -2.30 131.24NiPSe -20.52 -4.38 215.55 U eff = 4 eV J (meV) J (meV) J (meV)MnPS -6.39 4.11 -1.90CoPS -1.71 -0.07 5.23NiPS -5.17 -0.78 34.93NiPSe -5.80 -0.45 53.13 note that for CuPS , the AFM Neel and AFM zigzag states arenot metastable with different U eff in our calculation, whichmeans that J ≤ . In Fig.7 (b), it is clear that J reachesthe maximum value in NiPX , which can be easily under-stood as the hall-filling of e g orbitals maximize the super-superexchange interaction. V. THE TWO-ORBITAL T-J MODEL AND DOPINGPHASE DIAGRAM FOR NIPX -20 -20-10 2010 2010-10 E(meV) E(meV) D O S ( A . U . ) |∆ | (a) Г K M M K ГK`/2 K`/2 + d ◇ s +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇◇ + d ◇ s (b)(c) (d) FIG. 8. (a) and (b) plot the amplitude distributions of the supercon-ducting gap G ( k ) of the d ± id -wave and the extended s -wave statesrespectively. (c) and (d) are the density of states in the d ± id -waveand the extended s -wave states at a doping level, n=3.8 and 4.2, re-spectively, by taking ∆ = 0 . and J = 40 meV . From the above analysis and the known experimentalfacts[6], it is clear that NiPX must belong to strongly cor-related electron systems. The band width of the two e g or-bitals is only about 1eV, much less than the band gaps in theirAFM zigzag states. Moreover, as we showed above, the ex-perimental band gaps are close to the theoretical results whenwe take U eff ∼ eV , which is much larger than the bandwidth as well. Thus, the magnetic order is caused by thestrong electron-electron correlation.Following the standard argument, NiPX , just like manyother strongly correlated electron systems, must be a Mott in-sulator. As U eff is much larger than the band width, we cantake the large U limit to derive a t-J type of model. For NiPX ,a minimum two-orbital t-J model can be written as H tJ = ˆ P H ˆ P + (cid:88) ij H J,
In summary, we have shown that the Ni-based TMTs areclose to a strongly correlated quadruple-layer graphene andare Dirac materials described by a two-orbital model withthe strong electron electron interaction. The main electronickinematics and magnetic interactions exist with unusual longrange distance between two third nearest neighbour Ni atoms,which stems from the super-super exchange mechanism. Withthis underlining electronic structure, the materials provide asimple and ideal playground to investigate strong correlationphysics.The two orbital model can be viewed as a natural ex-tension of the single orbital model in the conventional hightemperature superconductors, cuprates. Recently, materi-als with both active e g orbitals have gained much attention.Ba CuO δ [36], synthesized under high pressure, is likely anextremely heavy hole doped cuprates. As the Jahn-Teller dis-tortion of the CuO octahedron causes a shorter Cu-O bondalong c-axis than the in-plane ones, both e g orbitals become important[37]. The single CuO layer grown by MBE alsohas a similar electronic structure[38]. In La Ni Se O , a re-cent theoretically proposed candidate of Ni-based high tem-perature superconductors[39], the low energy physics is alsoattributed to the two e g orbitals of Ni atoms. If the supercon-ductivity is determined to arise in the two orbital model, theNi-based TMTs, together with these materials, can provide usmuch needed information to solve the elusive mechanism ofhigh temperature superconductivity.Comparing with recent artificial graphene systems, inwhich new bands are created to reduce the kinetic energy sothat the effect of the weak electron-electron correlation canarise in the standard graphene[11], the Ni-based TMTs aresimply in the other limit (Mott limit) with the strong electron-electron correlation. We can consider to increase the kineticenergy, namely the band width, to enter the Mott transitioncritical region. In general, the band width can be increased byapplying pressure or by atom substitutions. During these pro-cesses, the angles between Ni-S/Se-Ni can be changed greatlyas well, which can lead to different intriguing physics. Forexample, in this case, more spin-orbital superexchange inter-action terms may be important, which can lead to exciting in-terplay between orbital and spin degrees of freedom.Experimentally, electron or hole carriers have not been in-troduced to Ni-based TMTs. Recently, modern gating tech-nology can induce carriers to a variety of two dimensionalmaterials[40]. The Ni-based TMTs can be an important play-ground for this modern technology. Acknowledgement:
We thank Hong Jiang, Xinzheng Li, Xi-anxin Wu, Yuechao Wang, Shengshan Qin and Dong Luanfor useful discussion. The work is supported by the Min-istry of Science and Technology of China 973 program(Grant No. 2015CB921300, No. 2017YFA0303100, No.2017YFA0302900), National Science Foundation of China(Grant No. NSFC-11334012), and the Strategic Priority Re-search Program of CAS (Grant No. XDB07000000). Q.Z.acknowledges the support from the International Young Sci-entist Fellowship of Institute of Physics CAS (Grant No.2017002) and the Postdoctoral International Program (2017)from China Postdoctoral Science Foundation.
Appendix A: The explicit form of effective Hamiltonian
The tight-binding effective Hamiltonian H in Eq.1 is a × matrix. The explicit form of its elements is given by H = t SNNxzxz (2 cos k x + cos k x √ k y t SNNyzyz cos k x √ k y ,H = t SNNxzxz ( −√ k x √ k y t SNNyzyz ( √ k x √ k y i (4 t SNNxzyz sin k x k x − cos √ k y ,H = 12 e − iky √ (cid:18) e i √ k y cos (cid:18) k x (cid:19) ( t NNxzxz + 3 t NNyzyz ) + e i √ k y (cos( k x )( t T NNxzxz + 3 t T NNyzyz ) + 2 t NNxzxz ) + 2 t T NNxzxz (cid:19) ,H = 12 i √ e − iky √ (cid:18) sin (cid:18) k x (cid:19) ( t NNxzxz − t NNyzyz ) − e i √ k y sin( k x )( t T NNxzxz − t T NNyzyz ) (cid:19) ,H = t SNNyzyz (2 cos k x + cos k x √ k y t SNNxzxz cos k x √ k y ,H = 12 e − iky √ (cid:18) e i √ k y cos (cid:18) k x (cid:19) (3 t NNxzxz + t NNyzyz ) + e i √ k y (cos( k x )(3 t T NNxzxz + t T NNyzyz ) + 2 t NNyzyz ) + 2 t T NNyzyz (cid:19) , (A1)with H = H , H = H , H = H ∗ and H = H by symmetry. These hopping parameters are given in Table.4in the main text. Appendix B: Formulation of the slave boson meanfield
We provide the detailed procedure for the slave bosonmeanfield method on the Hamiltonian Eq. (8). In ourtwo orbital model, the two e g orbitals are degenerate sothat they have the identical occupancy. In the slave-bosonapproximation[35], the same occupancy for all the orbitalsleads to the same renormalization for all the hopping inter-action. Namely, we have ˆ P H ˆ P = | n − | H , (B1)in which n = 4 ± x , x is the doped electron (+) or hole(-) per atom. n = 4 represents the half filling, where thekinetic energy vanishes and the Hamiltonian reduces to pureHeisenberg exchange interaction. The exchange term H J, (cid:104) ij (cid:105) in Eq. (8) can be decoupled in superconducting and magneticchannels.In superconducting channel, it is − J (cid:88) ij (cid:16) (cid:104) ∆ † ij (cid:105) ∆ ij + ∆ † ij (cid:104) ∆ ij (cid:105) (cid:17) + E s , (B2)with pairing matrix ∆ † ij ≡ σa † i ¯ ij,σ a † j ¯ ij, ¯ σ with a i ¯ ij,σ =cos θ ij a ix,σ + sin θ ij a iy,σ and the constant part E s = J (cid:88) ij (cid:16) (cid:104) ∆ † ij (cid:105)(cid:104) ∆ ij (cid:105) + n i ¯ ij,σ n j ¯ ij,σ (cid:17) , (B3)with the spin indices ¯ σ = − σ = ± and link ij the TNN bond. The same spin coupling density term n i ¯ ij,σ n j ¯ ij,σ term in E s would be decoupled as bond hop-ping (cid:104) a † i ¯ ij,σ a j ¯ ij,σ (cid:105) a i ¯ ij,σ a † j ¯ ij,σ according to the approach inliterature[33]. Those bond hopping terms effectively renor-malize the third neighbor σ -bond hoppings to affect the posi-tion and shape of Fermi surface. The hopping parameters and Fermi surfaces will be renormalized as the doping varies inour model. Thus as a simple illustrating of the phase diagramunder slave boson method, those bond hopping terms wouldbe simply taken as the density correlation. Put the renormal-ized hopping and decoupled exchange terms together, in thetwo orbital momentum space, we obtain (cid:80) k Ψ † k h Bk Ψ k + E s with h Bk given in Eq. (10).In the magnetic channel, the J dominated spin exchangeinteraction prefers the zigzag AFM state as shown in Fig.3. Inthe mean field level, the magnetic order in the zigzag patternis given by < S i > = ( − i m ˆ z with i = i a + i a .It is important to point out that as only the σ bond orbitalsare considered, the spin exchange is orbital selective. That isto say, the spin operator, S i, ¯ ij = a + i ¯ ijµ σ µν a i ¯ ijν with σ thevector of three Pauli matrices. As a result (cid:104) S zi, ¯ ij (cid:105) = cos θ ¯ ij m ixx + sin θ ¯ ij m iyy (B4) + cos θ ¯ ij sin θ ¯ ij ( m ixy + m iyx ) , with m iαβ ≡ (cid:104) σ ( a † iασ a iβσ )) (cid:105) / − i m αβ . Up to a con-stant term, H J, (cid:104) ij (cid:105) is decoupled as (cid:88) H J, (cid:104) ij (cid:105) = (cid:88) k σa † kασ M αβ a k + Q βσ (B5)with Q = G / marked in Fig.5 as the ordered magneticwavevector. The scattering matrix from k to k + Q in e g orbitals space is M = − J (cid:18) m xx + m yy m xy + m yx m xy + m yx m xx + 3 m yy (cid:19) . (B6)Defining Ψ † Mk = (Ψ † k , Ψ † k + Q ) , the meanfield Hamiltoniancan be written as H = (cid:88) k ∈ rBZ (cid:16) Ψ † M A k Ψ Mk + tr ( h k + h k + Q ) − µ (cid:17) + E sm , (B7)with rBZ represent the reduced Brillouin zone due to magnetic0cell and A k is a × matrix as A k = (cid:18) h Bk I (cid:78) I (cid:78) M I (cid:78) I (cid:78) M † h Bk + Q (cid:19) (B8) E sm N = 3 J + 3 J m xx + 3 m yy + 2 m xx m yy (B9) + ( m xy + m yx ) + n ) . In A k , the first I is for particle-hole space and the secondis for A-B sublattice. The h Bk is given in Eq. (10). Theself-consistency of the chemical potential is also taken intoconsideration for a fixed doping. It is easy to show that n k + n k + Q =8 + tr (Σ (cid:104) Ψ Mk Ψ † Mk (cid:105) ) (B10) =8 − tr (Σ U k f (Λ k ) U † k ) , in which Σ is the × stagger matrix Σ ≡ − I (cid:78) σ (cid:78) I , U † k A k U k = Λ k , diagonalizes the A k and f (Λ k ) is the Fermidistribution function.It is worthy to mention that the numerical result indicatesslight difference between the intra-orbital magnetic orders, m xx and m yy , due to the rotation symmetry breaking in theAFM zigzag state. The inter-orbital magnetic orders, m xy and m yx , are very small and can be ignored in the meanfield solu-tion. [1] Novoselov, K. S. et al. Electric field effect in atomically thincarbon films.
Science , 666–669 (2004).[2] Wang, F. et al.
New frontiers on van der waals layered metalphosphorous trichalcogenides.
Adv. Funct. Mat. , 1802151(2018).[3] Susner, M. A., Chyasnavichyus, M., McGuire, M. A., Ganesh,P. & Maksymovych, P. Metal thio-and selenophosphates asmultifunctional van der waals layered materials. Adv. Mater. , 1602852 (2017).[4] Friedel, C. Compt. Rend. (1894).[5] Ferrand, L.
Bull. Soc. Chim. , 115 (1895).[6] Le Flem, G., Brec, R., Ouvard, G., Louisy, A. & Segransan, P.Magnetic interactions in the layer compounds MPX (M = Mn,Fe, Ni; X= X, Se). J. Phys. Chem. Sol. , 455–461 (1982).[7] Wang, Y. et al. Emergent superconductivity in an iron-basedhoneycomb lattice initiated by pressure-driven spin-crossover.
Nat. Commun. , 1914 (2018).[8] Kim, So Yeun et al. Charge-Spin Correlation in van der WaalsAntiferromagnet NiPS . Phys. Rev. Lett. , 136402 (2018).[9] Wu, C., Bergman, D., Balents, L. & Sarma, S. D. Flat bands andwigner crystallization in the honeycomb optical lattice.
Phys.Rev. Lett. , 070401 (2007).[10] Wu, C. & Sarma, S. D. p x,y -orbital counterpart of graphene:Cold atoms in the honeycomb optical lattice. Phys. Rev. B ,235107 (2008).[11] Cao, Y. et al. Correlated insulator behaviour at half-filling inmagic-angle graphene superlattices.
Nature , 80 (2018).[12] Chittari, B. L. et al.
Electronic and magnetic properties ofsingle-layer MPX metal phosphorous trichalcogenides. Phys.Rev. B , 184428 (2016).[13] Hu, J. & Ding, H. Local antiferromagnetic exchange and col-laborative fermi surface as key ingredients of high temperaturesuperconductors. Scientific Reports , 381 (2012).[14] Hu, J. Identifying the genes of unconventional high temperaturesuperconductors. Science Bulletin , 561–569 (2016).[15] Seo, K., Bernevig, B. A. & Hu, J. Pairing symmetry in a two-orbital exchange coupling model of oxypnictides. Phys. Rev.Lett. , 206404 (2008).[16] Leb`egue, S. & Eriksson, O. Electronic structure of two-dimensional crystals from ab initio theory.
Phys. Rev. B ,115409 (2009).[17] Lee, J.-U. et al. Ising-type magnetic ordering in atomically thin FePS . Nano letters , 7433–7438 (2016).[18] Allmann, R. & Hinek, R. The introduction of structure typesinto the inorganic crystal structure database icsd. Acta. Crystal-logr. A , 412–417 (2007).[19] Ouvrard, G., Brec, R. & Rouxel, J. Structural determination ofsome MPS layered phases (M= Mn, Fe, Co, Ni and Cd). Mat.Res. Bull. , 1181–1189 (1985).[20] Klingen, W., Eulenberger, G. & Hahn, H. ¨Uber hexathio-und hexaselenohypodiphosphate vom typ M II P X . Naturwis-senschaften , 229–230 (1968).[21] Rao, R. R. & Raychaudhuri, A. Magnetic studies of a mixedantiferromagnetic system Fe − x Ni x PS . J. Phys. Chem. Sol. , 577–583 (1992).[22] Brec, R., Ouvrard, G., Louisy, A. & Rouxel, J. Proprietes struc-turales de phases M (II)PX (X= S, Se ). In Annales de chimie–science des mat´eriaux , 499–512 (1980).[23] Kresse, G. & Furthm¨uller, J. Efficiency of ab-initio total energycalculations for metals and semiconductors using a plane-wavebasis set.
Comp. Mat. Sci. , 15–50 (1996).[24] Kresse, G. & Joubert, D. From ultrasoft pseudopotentials tothe projector augmented-wave method. Phys. Rev. B , 1758(1999).[25] Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradientapproximation made simple. Phys. Rev. Lett. , 3865–3868(1996).[26] Mostofi, A. A. et al. wannier90: A tool for obtainingmaximally-localised wannier functions. Comp. Phys. Commun. , 685–699 (2008).[27] Dudarev, S. L., Botton, G. A., Savrasov, S. Y., Humphreys, C. J.& Sutton, A. P. Electron-energy-loss spectra and the structuralstability of nickel oxide: An LSDA+U study.
Phys. Rev. B ,1505 (1998).[28] Li, Y. et al. Robust d -wave pairing symmetry in multiorbitalcobalt high-temperature superconductors. Phys. Rev. B ,024506 (2017).[29] Scalapino, D. J. A common thread: The pairing interaction forunconventional superconductors. Rev. Mod. Phys. , 1383–1417 (2012).[30] Wildes, A. et al. Magnetic structure of the quasi-two-dimensional antiferromagnet nips 3.
Phys. Rev. B , 224408(2015).[31] Castellani, C., Natoli, C. R. & Ranninger, J. Magnetic structure of V O in the insulating phase. Phys. Rev. B , 4945–4966(1978).[32] Kugel, K. & Khomskii, D. The Jahn-Teller effect and mag-netism: transition metal compounds. Physics-Uspekhi , 231(1982).[33] Wang, Q.-H., Lee, D.-H. & Lee, P. A. Doped t − J model ona triangular lattice: Possible application to Na x CoO · y H O and Na − x TiO . Phys. Rev. B , 092504 (2004).[34] Brinckmann, J. & Lee, P. A. Renormalized mean-field theoryof neutron scattering in cuprate superconductors. Phys. Rev. B , 014502 (2001).[35] R¨uegg, A., Indergand, M., Pilgram, S. & Sigrist, M. Slave-boson mean-field theory of the mott transition in the two-band hubbard model. Eur. Phys. J. B , 55–64 (2005).[36] Li, W. et al. A new superconductor of cuprates with uniquefeatures. arXiv:1808.09425 (2018).[37] Maier, T. A., Berlijn, T. & Scalapino, D. J. d -wave and s ± pairing strengths in Ba CuO δ . arXiv:1809.04156 (2018).[38] Jiang, K., Wu, X., Hu, J. & Wang, Z. Nodeless high-T c superconductivity in highly-overdoped monolayer CuO . arXiv:1804.05072 (2018).[39] Le, C., Zeng, J., Gu, Y., Cao, G.-H. & Hu, J. A possible familyof Ni-based high temperature superconductors. Science Bulletin , 957 – 963 (2018).[40] Ye, J. T. et al. Liquid-gated interface superconductivity on anatomically flat film.
Nat. Mater.9