Mott metal-insulator transitions in pressurized layered trichalcogenides
MMott metal-insulator transitions in pressurized layered trichalcogenides
Heung-Sik Kim,
1, 2
Kristjan Haule, and David Vanderbilt Department of Physics and Astronomy, Rutgers University, Piscataway, New Jersey 08854-8019, USA Department of Physics, Kangwon National University, Chuncheon 24341, Korea
Transition metal phosphorous trichalcogenides, M P X ( M and X being transition metal and chalcogen el-ements respectively), have been the focus of substantial interest recently because they are unusual candidatesundergoing Mott transition in the two-dimensional limit. Here we investigate material properties of the com-pounds with M = Mn and Ni employing ab-initio density functional and dynamical mean-field calculations,especially their electronic behavior under external pressure in the paramagnetic phase. Mott metal-insulatortransitions (MIT) are found to be a common feature for both compounds, but their lattice structures show dras-tically different behaviors depending on the relevant orbital degrees of freedom, i.e. t or e g . Under pressureMnPS undergoes an isosymmetric structural transition in monoclinic space group by forming Mn-Mn dimersdue to the strong direct overlap between the neighboring t orbitals, accompanied by a significant volume col-lapse and a spin-state transition. In contrast, NiPS and NiPSe , with their active e g orbital degrees of freedom,do not show a structural change at the MIT pressure or deep in the metallic phase within the monoclinic sym-metry. Hence NiPS and NiPSe become rare examples of materials hosting electronic bandwidth-controlledMott MITs, thus showing promise for ultrafast resistivity switching behavior. Since the first identification of the Mott metal-insulatortransition (MIT) by Mott and Peierls in 1937 [1] and the sug-gestion of the canonical Hubbard model in 1963 [2], manysystems showing the Mott MIT have been found. They can bebroadly classified into two categories, a ) the filling-controlledMITs, such as in the doped cuprates [3], or b ) the bandwidth-controlled MITs, such as in the rare-earth nickelates R NiO ( R being a rare-earth element) [4, 5] or vanadium oxidesV O [6–8] and VO [6, 9]. Considering applications to elec-tronic resistive switching devices, the filling-controlled MITof type ( a ) is not favorable due to the inevitable strong inho-mogeneity at the atomic scale introduced by the chemical dop-ing. The bandwidth-controlled MITs of type (b), on the otherhand, are typically coupled strongly to the structural degreesof freedom, as for examples in the bond disproportionation be-tween the short and long Ni-O bonds in R NiO [10–13] andthe dimerization of vanadium atoms in VO [14, 15]. Suchinvolvement of slow lattice dynamics in the MIT is also notfavorable for fast switching. Hence, systems with electronic bandwidth-controlled MITs ( i.e. weak or no lattice distortionsinvolved) are desirable for fast resistive switching [16, 17].Surprisingly, there are very few solids that are known toundergo purely electronic and bandwidth-controlled MITs, aswas originally envisioned by Hubbard. From the theoreticalside, there is growing evidence that starting from the metal-lic side, the MIT would not have occurred in any of abovethree systems ( R NiO , V O , and VO ) in the absence of a si-multaneous structural distortion. Even less common are suchtransitions in two-dimensional materials, which might be use-ful for ultrathin electronic and spintronic applications, and toour knowledge there is no known example of an electronicallydriven MIT among the van der Waals (vdW) materials exceptthe recently discovered Mott phase and superconductivity intwisted bilayer graphene [18, 19].Here we propose new candidates for the electronicbandwidth-controlled MIT without significant structural dis-tortion among the emerging class of two-dimensional vdW materials. Our target is a series of transition metal phospho-rous trichalcogenides M P X ( M = Mn, Ni, X = S, Se) [20–22]. To incorporate the electronic and structural degrees offreedom on an equal footing, we employ the state-of-the-artembedded dynamical mean-field theory combined with den-sity functional theory (eDMFT), which implements forceson atoms and allows relaxation of internal atomic coordi-nates [23]. For the optimization of the size and shape ofunit cells we use density functional theory (DFT) augmentedby the on-site Coulomb repulsion U (DFT+ U ), after whichoptimizations of internal atomic coordinates are performedboth in eDMFT and DFT+ U yielding consistent results[24].We mainly focus on paramagnetic phases of MnPS andNiP { S,Se } above their N´eel temperatures ( T N = 78 and 154K for MnPS and NiPS respectively [25–27]), with disor-dered local Mn ( d ) S = 5 / and Ni ( d ) S = 1 mo-ments, although the behavior of their MITs in the magneticphases is discussed in the Supplementary Material (SM). Wewill show that the recently discovered MIT in MnPS falls un-der the family of transitions coupled to structural changes, inwhich the dimerization plays a crucial role, therefore bear-ing a resemblance to the MIT in VO [20]. On the otherhand, theoretical simulations in NiPS and NiPSe suggestthat the MIT in these two vdW compounds occurs at evenlower pressure, and does not involve a simultaneous structuraltransition. Therefore they become rare examples of electronicbandwidth-controlled transitions with a potential for very fastresistive switching. Crystal structures versus pressure.
Fig. 1(a) and (b) showDFT+ U results on the pressure-induced change of the threelattice parameters ( a/a , b/b , and c/c , where { a, b, c } de-note their zero-pressure values) for NiPS and MnPS , respec-tively. Note that here we focus on the monoclinic C /m struc-ture as shown in Fig. 1(c). Because both of the compoundsare vdW-type layered systems, the inter-plane lattice param-eter c shows a steeper decrease compared to the in-plane a and b , and the three-fold symmetry within each layer forces a r X i v : . [ c ond - m a t . s t r- e l ] D ec b/b c/c
0 20 40 60 80 0.75 0.8 0.85 0.9 0.95 1a/a b/b c/c (a) NiPS Pressure (GPa) (b)
MnPS V/V
0 20 40 60 80 0.75 0.8 0.85 0.9 0.95 1a/a b/b c/c Insulator Metal P c in exp. Mn dimer
Mn (S~5/2) / Ni
Mn (S~1/2) P S (c) (d) xyz d xy d xy abc ac FIG. 1. (a,b) Evolution of DFT+ U -optimized lattice parameters ( a , b , and c as depicted in (c,d)) as a function of pressure, where (a) and (b)panels show results from NiPS and MnPS respectively. Inset in (b) shows a volume versus pressure plot for MnPS , where black and graycurves represent ground-state and metastable structures respectively. Thick vertical dashed lines in both plots indicate the values of criticalpressure where the MIT happens. The thin dotted line in (a) shows the pressure where the a/a and b/b begin to branch in NiPS ( i.e. b (cid:54) = √ a ) due to the enhanced monoclinicity by pressure. Note that, the critical pressure for MnPS reported in Ref. 20 is around 30 GPa, asdepicted in the figure. (c) Crystal structures for NiPS and MnPS at the ambient pressure. (d) MnPS structure when P >
64 GPa, where theMn-dimer is formed parallel to b . b (cid:39) √ a in the low- P regime. The resulting volume de-crease under pressure is substantial: 40% of volume reduc-tion at ∼ GPa compared to the ambient pressure volume, asshown in the inset of Fig. 1(b).Both compounds show MIT and structural phase transitionsunder pressure, but the nature of their transition is drasticallydifferent. As shown in Fig. 1(a), the MIT and the structuraltransition in NiPS occur at very different pressures, around31 and 57 GPa respectively, while they coincide in MnPS .Remarkably, theoretical simulations suggest that the MIT inNiPS accompanies no significant structural distortion (dis-continuous structural changes, for example), and is thus arare example of an electronically driven bandwidth-controlledMIT. On the other hand, in MnPS the isosymmetric struc-tural transition ( i.e. , structural transition within the same spacegroup symmetry) with a volume collapse at GPa is crucialfor the occurrence of the MIT, hence the transition is betterclassified as the structurally assisted MIT (see Fig. 1(b)). Wenote that the theoretical critical pressure of GPa is some-what overestimated compared to the experimentally reportedvalue of ∼ GPa [20]. However, we show in the SM thatwithin eDMFT, spinodal lines extend down to a much lowerpressure of GPa with a much reduced energy barrier be-tween the metallic and insulating solutions compared to theDFT+ U results. Inclusion of the phonon free energy and thelattice zero-point energy, which is neglected here, could thenmove the position of the transition significantly (see the SMfor further details).In addition to the volume collapse, mostly from the dis-continuous change of a , DFT+ U simulations of MnPS showa Mn-Mn dimerization along the b -direction with the tiltingof the P dimer as shown in Fig. 1(d). The Mn-Mn bondlengths between the dimer and non-dimer bonds are 2.42 and 3.10 ˚A at 63 GPa, respectively, which is a rather large differ-ence. This Mn dimer formation is attributed to the direct d - d overlap between the Mn t orbitals, pointing directly towardsthe nearest-neighbor Mn as shown in the inset of Fig. 1(d).Note that the previous experimental study suggested the for-mation of Mn zigzag chains in the high-P phase [20], in con-trast to our DFT+ U and eDMFT results.NiPS , on the other hand, shows no such intermetallicdimerization or chain formation at the MIT or beyond thestructural transition pressure, because the partially-filled Ni e g orbitals point towards the S atoms. This makes NiPS more sensitive to the p - d hybridization, yielding a smallerMIT pressure in NiPS compared to MnPS . Such a starkcontrast between NiPS and MnPS , originating from the dif-ference in their orbital physics, affects the nature of the struc-tural behavior of their MIT as shown below. Note that thestructural transition in NiPS at 57 GPa is not orbital in natureand comes purely from the reduced interlayer distance and thelarge overlap between the layers.We note that the two compounds show markedly differentpressure dependence of the lattice parameters even before thestructural transition. By comparing Fig. 1(a) and (b) it is ev-ident that the compression of c under pressure is stronger inNiPS than in MnPS ; while a/a − c/c in MnPS at 60GPa is about 0.03 (See Fig. 1(b)), in NiPS it is about 0.15(Fig. 1(a)) despite the similar volume change. In other words,it is much easier to compress NiPS along the layer-normal di-rection compared to MnPS . Because the kinetic energy scaleset by the hopping integrals between the t (for MnPS ) and e g (for NiPS ) orbitals show different anisotropy, the t and e g yield strong in-plane d - d and inter-plane d - p - p - d overlapsrespectively (See SM). As a result, while the t orbitals fa-vor in-plane compression for the larger in-plane kinetic energy g ’a e g -4 -2 0 2 4 P D O S ( s t a t e s / e V / N i ) Energy (eV) (a)
P = 0 GPa30.4 GPa31.8 GPa88.0 GPa | S z | ( μ B ) (c) P r o b a b ili t y (b) Insulator Metal
Pressure (GPa) |2S z | (DMFT) P(d , |S z |=0)P(d , |S z |=1/2)P(d , |S z |=1) FIG. 2. (a) Projected density of states (PDOS) in the paramagneticphase of NiPS at T = 232K, calculated by eDMFT with the in-creasing pressure from the ambient condition (top panel) to 88 GPa(lowest). (b) Monte Carlo probabilities for the d | S z | = 0 (purpledashed line), d , | S z | = 1/2 (red dash-dotted), and d | S z | = 1(blacksolid) as a function of pressure. (c) Pressure dependence of the sizeof PM Ni spin moment | S z | from PM eDMFT results at T = 232K.Note a cusp at P = 30.4 GPa where the MIT happens. gain, the e g orbitals prefer to reduce the inter-plane distance,yielding the tendency shown in Fig. 1(a) and (b). Electronic MIT in NiPS . Below we take a closer look intothe nature of the MIT in NiPS . Note that all the spectrapresented hereafter are eDMFT results, where the DFT+ U -optimized cell parameters and estimated pressure values areemployed. Fig. 2(a) shows projected densities of states(PDOS) of NiPS with varying pressure from 0 to 88 GPa. It is clear that the t states ( a g and e (cid:48) g ) are mostly occupied,while the e g states are partially filled and show a narrow dipat the Fermi level at 30.4 GPa. The self-energies of the e g orbitals show poles at the Fermi level (see SM), confirmingthe presence of the paramagnetic Mott phase. Previously, itwas suggested that NiPS is a negative charge-transfer (NCT)insulator with a d L configuration ( L denoting a S p -ligandhole) [21]. However, our eDMFT results show that when theNi occupation is close to n d ≈ , where the d L configu-ration is dominant, the Mott insulating state cannot be stabi-lized, i.e., the material is metallic. The experimentally ob-served Mott insulating behavior can only be achieved with theNi occupancy of n d ≈ , where the high-spin S = 1 con-figuration is dominant, i.e., corresponding to approximatelyhalf-filled e g states (see Fig. 2(b) for the probability distribu-tion in the insulating and metallic states). This observationis corroborated by X-ray absorption spectroscopy, indicatingthat NiPS is close to the NCT regime, but is still dominatedby the d S = 1 configuration, consistent with our eDMFTresults [28].Fig. 2(b) shows the valence histogram for the few most im-portant Ni d configurations versus pressure at T = 232K. TheMott insulating state is stable as long as the high-spin state( | S z | = 1) of the Ni- d configuration is dominant. Note that wereport S z values rather than S values, because of our choice ofan Ising-type approximation of the Coulomb interaction in theeDMFT impurity solver [29]. Around 31 GPa the | S z | = 1/2states (of d and d configurations) become equally probable,at which point the Mott state collapses and a narrow metallicquasiparticle peak appears (see the third panel in Fig. 2(a)).Despite the enhanced charge fluctuation, the change of Ni d -orbital occupation ( n d ) across the transition is negligible: n d = 8.15 and 8.19 at P = 0 and 88 GPa, respectively. The in-crease of charge fluctuations with increasing pressure has anadditional effect of unlocking the | S z | = 0 sector of the d con-figuration, which is favored in the itinerant low spin regime atlarge pressure. Note that the increase of the probability for | S z | = 0 at the expense of the | S z | = 1 state has a large effecton the size of the fluctuating moment | S z | , which is plottedin Fig. 2(c). Its zero-pressure value is around . µ B , whichis quite reduced from the maximum atomic value of µ B ,and once it is reduced below . µ B it drops very suddenlyand takes values of ≤ . µ B in the metallic state. We notethat the change between 30.4 and 31.8 GPa is abrupt, which islikely associated with a first-order transition, for which a co-existence of both solutions is expected. However the hystere-sis was not observed, probably because it is too narrow at thetemperature studied ( T = 232 K). We mention that the MITwas carefully checked by employing beyond-Ising Coulombinteraction terms (spin-flip and pair-hopping types) at a lowertemperature of T = 116 K, but neither a hysteretic behaviornor a discontinuity in the energy-volume curve were found,signifying very weak coupling between the lattice and chargedegrees of freedom in NiPS (see Sec. III. A in the SM andFig. S3 therein for more details). MIT driven by uniaxial pressure in NiPSe . While the crit- c=1.00c (P = 0 GPa)c=0.88c (P = 8.0 GPa)c=0.84c (P = 12.9 GPa)c=0.82c (P = 15.4 GPa)c=0.81c (P = 16.6 GPa)c=0.80c (P = 17.8 GPa) c =1.00 c ( σ z = 0 GPa) c =0.88 c ( σ z = 8.0 GPa) c =0.86 c ( σ z = 10.5 GPa) c =0.84 c ( σ z = 12.9 GPa) c =0.82 c ( σ z = 15.4 GPa) c =0.81 c ( σ z = 16.6 GPa) e g P D O S ( s t a t e s / e V / N i ) Energy (eV)
FIG. 3. NiPSe e g -PDOS in the presence of an uniaxial stress σ z along the layer-normal ˆ z -direction from PM eDMFT results at T =58 K (with varying σ z from 0 to 16.6 GPa). Blue and red curves arePDOS for insulating and metallic phases, respectively, when the MIThappens between σ z = 15.4 and 16.6 GPa. ical pressure for the MIT in NiPS can be reached in modernhigh-pressure experimental setups, a substitution of S by themore polarizable Se is expected to further reduce the criti-cal pressure. Therefore the recently synthesized NiPSe [30]can be a better candidate for realizing the pressure-driven MITcompared to NiPS . Moreover, the collapse of the interlayerdistance is expected to be sufficient to induce the MIT, whichcan even be achieved by the tip of an atomic-force micro-scope [31]. As a zeroth-order approximation, we simulatesuch layer-normal strain by varying the interlayer distancewith fixed in-plane lattice parameters, and allowing the in-ternal coordinates to relax within eDMFT. In Fig. 3 we showthe e g -PDOS of NiPSe where the MIT happens at the mod-est stress of ≤ σ z ≤ GPa at T = 58 K, suggestingNiPSe as another promising electronic bandwidth-controlledMott transition system among these layered vdW materials. Volume collapse and MIT in MnPS . We now address thevolume-collapse transition in MnPS . We first checked thatonce the optimized lattice parameters from DFT+ U are em-ployed, both DFT+ U and eDMFT optimizations for the in-ternal coordinates yield practically the same result. Fig. 4(a)shows the evolution of the fluctuating Mn moment | S z | within the PM eDMFT at T = 580 K. We also show the or-dered magnetic moment M (per Mn) within DFT+ U in theN´eel-type antiferromagnetic ordered state [20, 25, 32] (greyline in the plot). Perhaps not surprisingly, the two methodsshow very similar behavior: an insulating state of almost max-imum spin S = 5/2 configuration in the low-to-intermediate-Pregime, and a metallic state with strongly reduced Mn mo-ments above P = 64 GPa. The inset of Fig. 4(a) schemati-cally depicts the spin-orbital configuration below and abovethe transition. In the presence of strong external pressure, theorbitally inert high-spin S = 5/2 configuration becomes ener-getically unstable, and the low-spin S = 1/2 with the partiallyfilled t orbital is stabilized with an octahedral volume col- M , | S z | ( μ B ) Pressure (GPa) (a)
P <
64 GPa P > 64 GPaMn t e g Mn 1 Mn 2 σσ * π / δ * π / δ e g P D O S ( s t a t e s / e V / N i ) (b) Energy (eV)
M (DFT+U)|2S| (DMFT) 1.2 1.4 25 30 35 d xy σ σ * π / δ * π / δ P = 0 GPa62 GPa 70 GPa g ’a e g -4 -2 0 2 4 xz/yzxye g FIG. 4. (a) Pressure dependence of the size of Mn spin moment Mfrom DFT+ U (gray line) and | S z | from paramagnetic eDMFT at T = 580K (black), where both results show the spin-state transition at63 GPa. The schematic spin-orbital configurations below and abovethe transition are shown in the inset. (b) PDOS from eDMFT, at P =0 GPa (upper panel), 62 GPa (center), and 70 GPa (lower). Note thatthe choice of d -orbital onto which the DOS are projected is differ-ent below and above the transition ( e (cid:48) g and a below, and d xz,yz,xy above 63 GPa). lapse [33]. As a result, the t open shell in this edge-sharinggeometry leads to a strong σ -like direct d - d overlap betweenthe nearest-neighboring (NN) Mn sites, resulting in a strongtendency toward Mn dimerization [34]. At ambient pressure,because of the weak inter-layer coupling, the three in-planeNN bonds are essentially equivalent to each other. In the high-pressure regime, however, the monoclinicity originating fromthe layer stacking is no longer negligible. Therefore, the NNbond parallel to the b direction, which becomes nonequivalentto the other bonds, dimerizes as shown in Fig. 1(d) (see Sec.III. B in the SM for more details on inclusion of the beyond-Ising Coulomb terms). Summary.
We report a theoretical study of the Mott MIT in-duced by external pressure in strongly correlated layered vdWmaterials. We comment that the Mott phases in other metaltrisulfides, such as FePS [35] or CoPS [36], are also of greatinterest because of their partially filled t shells even underambient conditions. Overall, this family of vdW-layered tran-sition metal trichalcogenides can be an excellent platform forthe study of strong electron correlations and their cooperationwith spin and lattice degrees of freedom. Acknowledgments:
We thank Matthew J. Coak, Michael O.Yokosuk, Nathan C. Harms, Kevin A. Smith, Sabine N. Neal,Janice L. Musfeldt and Sang-Wook Cheong for helpful discus-sions. The work was supported by NSF DMREF grant DMR-1629059. HSK thanks the support of the 2019 research grantfor new faculty members from Kangwon National University,and also the support of supercomputing resources includingtechnical assistances from the National Supercomputing Cen-ter of Korea (Grant No. KSC-2019-CRE-0036) [1] N. F. Mott and R. Peierls, Proc. Phys. Soc. , 72 (1937).[2] J. Hubbard, Proc. Royal Soc. Lond. A , 238 (1963).[3] J. G. Bednorz and K. A. M¨uller, Z. Phys. B , 189 (1986).[4] J. B. Torrance, P. Lacorre, A. I. Nazzal, E. J. Ansaldo, andC. Niedermayer, Phys. Rev. B , 8209 (1992).[5] M. L. Medarde, J. Phys. Condens. Matter , 1679 (1997).[6] F. J. Morin, Phys. Rev. Lett. , 34 (1959).[7] D. B. McWhan, J. P. Remeika, T. M. Rice, W. F. Brinkman, J. P.Maita, and A. Menth, Phys. Rev. Lett. , 941 (1971).[8] S. A. Carter, T. F. Rosenbaum, P. Metcalf, J. M. Honig, andJ. Spalek, Phys. Rev. B , 16841 (1993).[9] M. M. Qazilbash, M. Brehm, B.-G. Chae, P.-C. Ho, G. O.Andreev, B.-J. Kim, S. J. Yun, A. V. Balatsky, M. B. Maple,F. Keilmann, H.-T. Kim, and D. N. Basov, Science , 1750(2007).[10] K. Haule and G. L. Pascut, Sci. Rep. , 10375 (2017).[11] S. Johnston, A. Mukherjee, I. Elfimov, M. Berciu, and G. A.Sawatzky, Phys. Rev. Lett. , 106404 (2014).[12] H. Park, A. J. Millis, and C. A. Marianetti, Phys. Rev. Lett. , 156402 (2012).[13] T. Mizokawa, D. I. Khomskii, and G. A. Sawatzky, Phys. Rev.B , 11263 (2000).[14] W. H. Brito, M. C. O. Aguiar, K. Haule, and G. Kotliar, Phys.Rev. Lett. , 056402 (2016).[15] S. Biermann, A. Poteryaev, A. I. Lichtenstein, and A. Georges,Phys. Rev. Lett. , 026404 (2005).[16] D. Ruzmetov, G. Gopalakrishnan, C. Ko, V. Narayanamurti,and S. Ramanathan, J. Appl. Phys. , 114516 (2010).[17] Z. Yang, C. Ko, and S. Ramanathan, Annu. Rev. Mater. Res. , 337 (2011).[18] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo,J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxi-ras, R. C. Ashoori, and P. Jarillo-Herrero, Nature , 80 EP(2018).[19] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxi-ras, and P. Jarillo-Herrero, Nature , 43 EP (2018).[20] Y. Wang, Z. Zhou, T. Wen, Y. Zhou, N. Li, F. Han, Y. Xiao,P. Chow, J. Sun, M. Pravica, A. L. Cornelius, W. Yang, andY. Zhao, J. Am. Chem. Soc. , 15751 (2016).[21] S. Y. Kim, T. Y. Kim, L. J. Sandilands, S. Sinn, M.-C. Lee,J. Son, S. Lee, K.-Y. Choi, W. Kim, B.-G. Park, C. Jeon, H.-D.Kim, C.-H. Park, J.-G. Park, S. J. Moon, and T. W. Noh, Phys.Rev. Lett. , 136402 (2018). [22] J.-G. Park, J. Phys. Condens. Matter , 301001 (2016).[23] K. Haule, J. Phys. Soc. Jpn. , 041005 (2018).[24] For DFT+DMFT calculations we use the RutgersDFT+Embedded DMFT code [37], while for DFT andDFT+U calculation we used the Vienna Ab-initio
Simula-tion Package [38, 39]. Dependency to different choices ofexchange-correlation functionals, including van der Waalsfunctionals, was also checked. For computational details referto the Supplementary Material (SM), where we discussed thethe difference between the U -values employed in our eDMFTcalculations and other studies[40, 41].[25] R. Brec, Solid State Ion. , 3 (1986).[26] A. R. Wildes, H. M. Rønnow, B. Roessli, M. J. Harris, andK. W. Godfrey, Phys. Rev. B , 094422 (2006).[27] A. R. Wildes, V. Simonet, E. Ressouche, G. J. McIntyre,M. Avdeev, E. Suard, S. A. J. Kimber, D. Lanc¸on, G. Pepe,B. Moubaraki, and T. J. Hicks, Phys. Rev. B , 224408 (2015).[28] B. Pal, H.-S. Kim, K. Haule, D. Vanderbilt, J. Chakhalian, andJ. Freeland, (2018), to be submitted.[29] This approximation leads to some mixing between S = 0 and 1states, but is not expected to change qualitative aspects of theresults.[30] J.-Q. Yan, B. C. Sales, M. A. Susner, and M. A. McGuire, Phys.Rev. Mater. , 023402 (2017).[31] H. Lu, C.-W. Bark, D. Esque de los Ojos, J. Alcala, C. B. Eom,G. Catalan, and A. Gruverman, Science , 59 (2012).[32] K. Okuda, K. Kurosawa, S. Saito, M. Honda, Z. Yu, andM. Date, J. Phy. Soc. Jpn. , 4456 (1986).[33] J. Kuneˇs, A. V. Lukoyanov, V. I. Anisimov, R. T. Scalettar, andW. E. Pickett, Nat. Mater. , 198 EP (2008).[34] S. V. Streltsov and D. I. Khomskii, Proc. Natl. Acad. Sci. U.S.A. , 10491 (2016).[35] J.-U. Lee, S. Lee, J. H. Ryoo, S. Kang, T. Y. Kim, P. Kim, C.-H.Park, J.-G. Park, and H. Cheong, Nano Letters , Nano Letters , 7433 (2016).[36] A. R. Wildes, V. Simonet, E. Ressouche, R. Ballou, and G. J.McIntyre, J. Phys. Condens. Matter , 455801 (2017).[37] K. Haule, C.-H. Yee, and K. Kim, Phys. Rev. B , 195107(2010).[38] G. Kresse and J. Hafner, Phys. Rev. B , 558 (1993).[39] G. Kresse and J. Furthm¨uller, Phys. Rev. B , 11169 (1996).[40] K. Haule, T. Birol, and G. Kotliar, Phys. Rev. B , 075136(2014).[41] S. Mandal, K. Haule, K. M. Rabe, and D. Vanderbilt, arXivpreprint arXiv:1907.10498 (2019).[42] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E.Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Phys. Rev.Lett. , 136406 (2008).[43] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys,and A. P. Sutton, Phys. Rev. B , 1505 (1998).[44] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. , 566 (1980).[45] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. ,3865 (1996).[46] P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka, andJ. Luitz, WIEN2k, An Augmented Plane Wave + Local Or-bitals Program for Calculating Crystal Properties (KarlheinzSchwarz, Techn. Universit¨at Wien, Austria, 2001).[47] K. Haule and G. L. Pascut, Phys. Rev. B , 195146 (2016).[48] K. Haule and T. Birol, Phys. Rev. Lett. , 256402 (2015).[49] K. Haule, Phys. Rev. B , 155113 (2007).[50] V. Grasso, F. Neri, P. Perillo, L. Silipigni, and M. Piacentini,Phys. Rev. B , 11060 (1991).[51] M. Dion, H. Rydberg, E. Schr¨oder, D. C. Langreth, and B. I.Lundqvist, Phys. Rev. Lett. , 246401 (2004). [52] K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist, and D. C.Langreth, Phys. Rev. B , 081101 (2010).[53] J. Klimeˇs, D. R. Bowler, and A. Michaelides, J. Phys. Condens.Matter , 022201 (2009).[54] J. Chen, A. J. Millis, and C. A. Marianetti, Phys. Rev. B ,241111 (2015).[55] A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Van-derbilt, and N. Marzari, Comput. Phys. Commun. , 2309(2014).[56] Q. Han, T. Birol, and K. Haule, Phys. Rev. Lett. , 187203(2018). COMPUTATIONAL DETAILSDensity functional theory calculations
For unit cell optimizations and relaxations of initial internalcoordinates, the Vienna ab-initio
Simulation Package (
VASP ),which employs the projector-augmented wave (PAW) basisset [38, 39], was used for density functional theory (DFT) cal-culations in this work. 340 eV of plane-wave energy cutoffand 8 × × k -grid sampling were employed.For the treatment of electron correlations within DFT, arevised Perdew-Burke-Ernzerhof exchange-correlation func-tional for crystalline solid (PBEsol) was employed[42], in ad-dition augmented by on-site Coulomb interactions for transi-tion metal d -orbitals within a simplified rotationally-invariantform of DFT+ U eff formalism[43]. − eV/ ˚A of force crite-rion was employed for structural optimizations. For test pur-pose, Ceperley-Alder local density approximation[44] and theoriginal PBE[45] functionals were also used.Structural relaxations for all compounds were performedin the presence of the DFT+ U eff (4 eV) on-site Coulomb in-teraction and a N´eel-type antiferromagnetic order[25], whichgives reasonable agreements of lattice parameters and gapsizes with experimentally observed values[20, 21]. It shouldbe mentioned that, without incorporating magnetism and U eff to open the gap, the volume is severely underestimated forboth compounds, especially ∼
20% smaller in MnPS . Thisobservation signifies the role of electron correlations in struc-tural properties of these compounds. Dynamical mean-field theory calculations
A fully charge-self-consistent dynamical mean-field method[37], implemented in Rutgers DFT+ Embedded DMFT (eDMFT) Functional code(http://hauleweb.rutgers.edu/tutorials/) which is com-bined with
WIEN K code[46], is employed for computationsof electronic properties and optimizations of internalcoordinates[47]. In DFT level the Perdew-Wang local densityapproximation (LDA) is employed, which was argued toyield the best agreement of lattice properties combinedwith DMFT[48]. 500 k -points were used to sample thefirst Brillouin zone with RK max = 7.0. A force criterionof 10 − Ry/Bohr was adopted for optimizations of internalcoordinates. A continuous-time quantum Monte Carlomethod in the hybridization-expansion limit (CT-HYB) wasused to solve the auxiliary quantum impurity problem[49],where the full 5 d -orbitals of Ni and Mn were chosen as ourcorrelated subspaces in a single-site DMFT approximation.For the CT-HYB calculations, up to 10 Monte Carlo stepswere employed for each Monte Carlo run. In most runstemperature was set to be 232K, but in some calculationswith high pressure it was increased up to 580K because of theincreased hybridization between the impurity and bath. -10to +10 eV of hybridization window (with respect to the Fermi level) was chosen, and U = 10 eV and J H = 1 eV of on-siteCoulomb interaction parameters were used for both Mn andNi d -orbitals. A simplified Ising-type (density-density termsonly) Coulomb interaction was employed in this work, and itwas tested that the use of full Coulomb interaction yields onlyquantitatively different results in terms of pressure-inducedevolution of electronic structures; see Sec. III for moredetails. A nominal double counting scheme was used, withthe d -orbital occupations for double counting corrections forNi and Mn were chosen to be 8 and 5, respectively.We comment that the choice of optimal values of theCoulomb interaction U is method-dependent. Other than theeDMFT approach chosen in this study, there are two widelyemployed first-principles methods using U to incorporateelectron correlations; (a) DFT+ U , and (b) DFT+DMFT withWannierized correlated orbitals. Both methods use smallervalues of U ( (cid:39) d orbitals in transition-metal compounds compared to eDMFT (10 eV)[40, 41]. First,unlike in DFT+ U , in DFT+DMFT formalisms all (local) dy-namic screening processes are included via exactly solvingthe many-body impurity problem. Being such screening pro-cesses explicitly treated within DFT+DMFT means that, theinput Coulomb interaction U should be closer to the bare one(only screened by the core and semi-core states). In DFT+ U ,on the contrary, one should use the screened U (whose valuesmaller than the DMFT U ) to compensate the missing screen-ing processes therein. Hence the U values employed in oureDMFT results are larger than the values used in DFT+ U eff calculations.Secondly, for DFT+DMFT with Wannierized correlated or-bitals, it is well known that the Wannier functions containa substantial amount of p character from anions (oxygen orchalcogen ions) if a narrow Wannierization energy windowthat contains only the correlated subspace is chosen. Conse-quently the critical U becomes much smaller because of themixing of p character, to be indeed of the order of the band-width W . Note that this approach is equivalent to solvingHubbard-type models, where only correlated orbital degreesof freedom are considered.Contrary to the aforementioned approaches, in the eDMFTformalism we are solving a generalized Anderson-lattice typeHamiltonian (actually the p - d type Hamiltonian), where theeffective U eff that could be compared with the Hubbard- U inthe Hubbard-type model is actually the p - d splitting. The ad-vantage of using such p - d type Hamiltonian in DFT+DMFT isevident; the U values in such models are much more system-independent for many transition-metal compounds, as demon-strated recently[40, 41]. Therein we established that a reason-able U for a wide range of transition-metal oxides within theeDMFT ( i.e. the p - d model with very localized d orbitals de-scribed above) is around 10 eV and J = 1 eV, and is muchmore universal than the U values in downfolded Hubbard-likemodels. CA-NMPBE-NMPS-NM CA-AFPBE-AFPS-AF
CA-NMPBE-NMPS-NM CA-AFPBE-AFPS-AF
CA-AFPBE-AFPS-AF
CA-AFPBE-AFPS-AF M n PS V / V ( d i m e n s i o n l e ss ) M n PS G a p ( e V ) N i PS V / V ( d i m e n s i o n l e ss ) N i PS G a p ( e V ) U e ff (eV) V/V (dimensionless) (a) Gap-Exp.Gap-Exp. (b)(c) (d) U e ff = - > e V U e ff = - > e V FIG. S5. Calculated unit cell volumes and band gap sizes of(a,b) MnPS and (c,d) NiPS under different choices of exchange-correlation (XC) functionals, the presence of antiferromagnetic or-der, and U eff values at the ambient pressure condition ( i.e. optimiz-ing cell volume without any volume/shape constraint). (a,c) showhow the cell volume depends on the choice of XC functionals, thepresence of antiferromagnetic order, and U eff values for each com-pound, and (b,d) show the size of gap as a function of cell volumewith respect to the choice of XC functionals. Therein 4 data pointsfor each XC functional, from bottom to top, represent results with U eff = 0, 2, 4, and 6 eV respectively. In the legend PS, NM, AFdenote PBEsol, nonmagnetic, and antiferromagnetic orders respec-tively. N´eel- and zigzag-type AF order were employed for MnPS and NiPS , respectively. V denotes experimental cell volume foreach compound[25], and horizontal black dashed lines in (b,d) showexperimentally measured gap sizes. EXCHANGE-CORRELATION AND VAN DER WAALSFUNCTIONAL DEPENDENCE WITHIN DFT+ U eff Dependence on exchange-correlation functionals and U eff -valuein DFT+ U eff results Fig. S5 shows how MnPS and NiPS behave under thechoice of different exchange-correlation functionals and thevalue of U eff , where experimental cell volumes and gap sizesare from Ref. 21, 25, and 50. Fig. S5(a) and (c) show howthe cell volume depends on the choice of exchange-correlation(XC) functionals, the presence of antiferromagnetic (AF) or-der, and the U eff values for each compound. We notice that i) the absence of AF order, which prevents the formation ofhigh-spin configurations in both compounds, yields signifi-cantly underestimated cell volumes in all cases. Such behav-ior is more evident in MnPS , where the absence of mag-netism leads to the low-spin configuration that favors inter-metallic bonding. ii) the PBEsol XC functional gives better PBEsol+U eff vdW-DFvdW-DF2vdW-optPBEvdW-optB88vdW-optB86b
PBEsol+U eff vdW-DFvdW-DF2vdW-optPBEvdW-optB88vdW-optB86b
PBEsol+U eff vdW-DFvdW-DF2vdW-optPBEvdW-optB88vdW-optB86b
PBEsol+U eff vdW-DFvdW-DF2vdW-optPBEvdW-optB88vdW-optB86b P r e ss u r e ( G P a ) P r e ss u r e ( G P a ) B a n d g a p ( e V ) B a n d g a p ( e V ) (a) (b)(c) (d) MnPS V / V PS vs. PNiPS V / V PS vs. P MnPS V / V PS vs. gapNiPS V / V PS vs. gap V / V PS V / V PS FIG. S6. Dependence of crystal and electronic structures on fivevan der Waals (vdW) functionals (see the text) for (a,b) MnPS and(c,d) NiPS . Colored symbols represent vdW functional results, andPBEsol+ U eff (4 eV) results (dark gray curves) are shown as a refer-ence. In panel (a) and (c), the calculated pressure P is shown as afunction of dimensionless volume ( V /V PS , where V PS is the com-puted ambient-condition volume with PBEsol+ U eff = 4 eV). In panel(b) and (d), band gap is shown as a function of V /V PS . agreement with experimental volume than CA or PBE, andusing U eff = 4 ∼ U eff setup producesthe best fit. The same conclusion can be also made from thegap-volume dependence shown in Fig. S5(b) and (d), wherethe use of PBEsol+ U eff (4 eV) yields the best fit of gap sizeand cell volume for both compounds. Overall, employingPBEsol+ U eff (4 eV) does seem reasonable for studying pres-surized MnPS and NiPS . van der Waals functional dependence Fig. S6 presents how the estimated pressure and the sizeof Kohn-Sham energy gap at a given volume depend onthe choice of van der Waals (vdW) functionals in MnPS and NiPS . The range of unit cell volume is chosen tobe . V PS ≤ V ≤ . V PS , where V PS is the ambient-pressure cell volume optimized with PBEsol+ U eff (4 eV).We employed 5 different vdW functionals implemented in VASP ; vdW-DF[51], vdW-DF2[52], optPBE, optB88, andoptB86b[53]. For both compounds, different functionals tendto give similar results, while NiPS shows more noticeablefunctional dependence compared to MnPS . It can be spec-ulated that, in NiPS with an open-shell Ni e g orbitals, theunquenched orbital degree of freedom makes the system abit more sensitive to the treatment of correlations. Neverthe-less the qualitative features we address in this manuscript, thepressure-driven insulator-to-metal transitions and their orbitaldependence, remain basically the same as shown in Fig. S6(b)and (d).Specifically we notice that vdW functionals, except vdW-DF and vdW-DF2 functionals, show reasonable agreementswith PBEsol results. vdW-DF and vdW-DF2 functionalstends to prefer larger volume ( i.e. larger P estimated at thesame volume compared to other functionals). This is becausethese functionals favor larger interlayer distances in the high-pressure regime than conventional XC functionals. Overall,even though use of different vdW functional induces somequantitative differences, it does not seem to change our mainconclusions in this work.We comment that, since vdW functionals favor larger cellvolume, the inclusion of them should enhance the magnitudeof critical pressures for structural/electronic transitions in bothcompounds, which is actually making the discrepancy be-tween the theoretical prediction and experimental observationslightly worse in MnPS . Hence we argue that PBEsol+ U eff can be a more reasonable choice in this pressurized setupwhere the direct orbital overlap between layers becomes sig-nificant. COMPARISON WITH DFT+ U eff AND DMFT
Table I presents the comparison between PBEsol+ U eff -and eDMFT-optimized atomic coordinates of MnPS andNiPS , both at ambient and high-pressure regimes. Hereambient and high-pressure results represent Mott-insulatingand weakly correlated metallic phases, respectively, for bothcompounds. In eDMFT calculations, as commented in themanuscript, optimized cell parameters a , b , c , and monoclinicangle β from PBEsol+ U eff were employed. This is due tothe absence of stress tensor formalism implemented in anyof full-potential linearized augmented plane wave codes, andin DFT+DMFT formalisms as well. Under this constraint,eDMFT-optimized atomic coordinates show very similar re-sults with PBEsol+ U eff ones despite different magnetizationconditions; paramagnetic order for eDMFT, and antiferro-magnetism (N´eel order for MnPS , zigzag order for NiPS )in PBEsol+ U eff .However, the validity of employing DFT+ U eff with a mag-netic order in optimizing crystal structures of paramagneticsystems, especially the cell parameters, may need to bechecked. This is because, like in MnPS as shown in themanuscript, some structural phase transitions are stronglycoupled to elastic deformations of the unit cell. One may evensuspect that the discrepancy between the predicted and ex-perimentally reported[20] critical pressures of the structuraltransition in MnPS might originate from the use of magneticDFT+ U eff in optimizing the unit cell size and shape.To resolve the issue mentioned above, energy landscapesfrom DFT+ U eff and eDMFT in the cell parameter space needto be compared with each other. Even though full structuralrelaxations may not be possible within eDMFT formalism,several trials to compare DFT+ U eff - and DMFT-optimizedstructures were performed, where DMFT calculations were Heung-Sik Kim - Mott MIT in Pressurized Layered Trichalcogenides Cell volume (a.u. ) To t a l e n e r g y ( e V / u . c . ) V exp (P=0) • Paramagnetic
DMFT results at
T=116K . • Accurate total-E vs. cell volume within DMFT . E ( e V / u . c . ) Cell volume (a.u. ) Energy (eV) D O S ( s t a t e s / e V / u . c . ) • No volume collapse/discontinuity found. • Crossover-like behavior due to temperature effects. (a) (b)(c)
FIG. S7. (a) Calculated energy versus volume plot for NiPS . Brightgreen and red symbols represent data points from PBEsol+ U eff andeDMFT results respectively. Sizes of error bars in eDMFT points( (cid:46) done in the paramagnetic configuration. NiPS According to our result presented in the manuscript, NiPS does not have a noticeable structural change at the MIT pres-sure of P (cid:39)
30 GPa, which is somewhat unusual. Hence, fora closer look on the structure-free MIT point, we computed atotal energy curve versus cell volume for paramagnetic (PM)eDMFT near the MIT. For more accurate results, the rotation-ally invariant form of the Coulomb interaction (spin-flip andpair-hopping included) was employed at a lower temperatureof T = 116 K. PM MIT is usually known to accompany a sud-den volume change (as reported in Ref. 35 in the manuscript,for example), which should be captured as a discontinuity inthe energy-volume curve at the MIT point. Figure S7 shows asummary of the results; note that all the data points were ob-tained from calculations started from scratch to capture boththe metallic and insulating phases (with optimized internalparameters). It can be seen that both the PBEsol+ U eff andDFT+DMFT data points remarkably collapse onto a singleBirch-Murnaghan energy-volume curve (Fig. S7(a)), and thatno discontinuity can be noticed near the MIT point (Fig. S7(b)and (c)). Note that the use of the additional Coulomb interac-tion with the spin-flip and pair-hopping terms (hereafter de-noted as ‘beyond-Ising’ terms) lowered the MIT critical pres-0 MnPS NiPS V /V PS P DFT (GPa) 0.0 49.8 0.0 88.0 a b c ( ˚A) 6.870 5.542 6.736 4.938 β (degree) 106.67 108.72 106.64 110.17DFT+ U eff eDMFT DFT+ U eff eDMFT DFT+ U eff eDMFT DFT+ U eff eDMFTMn ( g ) y g ) y i ) x i ) x z z i ) x i ) x z z j ) x j ) x y y z z and NiPS from DFT+ U eff and eDMFT results, both at ambient and high-pressureregimes. Ambient and high-pressure results represent Mott-insulating and weakly correlated metallic phases, respectively, for both com-pounds. PBEsol+ U eff = 4 eV was adopted for DFT+ U eff . Cell parameters ( a , b , c , and β ) optimized in DFT+ U eff calculations were employedin eDMFT ones. Nonzero components of Wyckoff positions of the C /m space group are shown. All eDMFT calculations were done at T = 232K except the high-pressure ( V = 0 . V PS ) MnPS one, where T = 580K was used for computational issues. V PS denotes the ambientpressure cell volume for both compounds, obtained with PBEsol+ U eff = 4 eV sure from 31 to 24 GPa. The crossover-like behavior can beattributed to the temperature effect, but due to the computa-tional cost issue the temperature could not be lowered below T = 116 K. MnPS As marked in Fig. 1(a) in the main text, the value of crit-ical pressure predicted by PBEsol+ U eff calculations is twicebigger than the one reported in Ref. 20 (63 vs. 30 GPa). Fora better understanding of this discrepancy, we perform calcu-lations with interpolating structures between the honeycombhigh-spin and the dimerized low-spin structures as shown inFig. S8. Because the computation of stress tensor is notyet available in the current eDMFT formalism, two constant-volume cuts at V = 0.57 and 0.54 V are taken for a total en-ergy comparison as shown in Fig. S8(a). For PBEsol+ U eff re-sults, high-spin and low-spin states are first converged in theirhoneycomb and dimerized structures respectively, and thenthe crystal structures are slowly distorted towards the otherside while maintaining the local minima spin states.Fig. S8(b) and (c) are relative total energies and size ofspin moments from the high-, low-spin PBEsol+ U eff , andparamagnetic eDMFT calculations at V = 0.57 V . A re-markable feature is, while the energy difference between thehigh- and low-spin ground states is 1.72 eV in PBEsol+ U eff ,it is 0.58 eV in eDMFT, which is almost one third of thePBEsol+ U eff value. Furthermore, while the high- and low- spin local minima states remain (meta)stable even after thestructural changes, as shown in Fig. S8(c), in eDMFT we havea spin-state crossover as the structure evolve from one limit toanother. This features persist at V = 0.54 V (Fig. S8(d,e)),where the height of the energy barrier from the high-spin tothe low-spin state (0.3 eV) is substantially suppressed (60meV) with the same spin-state crossover. These observationsshow that, the dynamical fluctuation effect inherent in eDMFTcauses mixing between different spin configurations, henceintroducing the crossover behavior shown in Fig. S8(c) and(e) and suppressing the energy differences. This observation isconsistent with a previous DFT+DMFT study on a spin-state-crossover molecule[54]. Note that, our analysis here does notexplicitly predict that the low-spin state is stabilized at lowerpressure in eDMFT results shown in Fig. 1(a) in the main text.However, since the energy difference between different statesis reduced to a fraction compared to PBEsol+ U eff results, thevalue of critical pressure might be reduced after the lattice freeenergy contribution and the zero-point fluctuation ignored inthis work are included.We comment that, due to the increased computational costof the CT-HYB impurity solver by the enhanced intermetallichybridization in the pressurized setup, all eDMFT data pointspresented in Fig. S8 were obtained at T = 580 K. Even with-out the beyond-Ising Coulomb terms, which significantly in-creases sign problems in the impurity solver stage, the temper-ature could not be lowered due to the computational cost issue.Therefore we could not check all the results with employing1 Dimerized Non-dimer V V PS V PS P PS a / a b / b c / c β ◦ ◦ Coulomb Ising Beyond-Ising Ising Beyond-IsingMn (4g) (0.0000, 0.1415, 0.0000) (0.0000, 0.1407, 0.0000) (0.0000, 0.1664, 0.0000) (0.0000, 0.1661, 0.0000)P (4i) (0.5256, 0.0000, 0.8230) (0.5278, 0.0000, 0.8238) (0.4358, 0.0000, 0.8064) (0.4358, 0.0000, 0.8063)S1 (4i) (0.8170, 0.0000, 0.6461) (0.8176, 0.0000, 0.6451) (0.2598, 0.0000, 0.6964) (0.2603, 0.0000, 0.6961)S2 (8j) (0.7538, 0.3145, 0.7219) (0.2542, 0.3142, 0.7227) (0.2231, 0.3205, 0.6980) (0.2228, 0.3201, 0.6984)
TABLE II. Optimized crystal structures of MnPS under different volume constraints ( V = 0 . V PS and . V PS , where V PS and P PS areoptimized ambient-pressure volume and exerted pressure at a fixed volume from PBEsol+ U eff (4 eV) results. Equilibrium lattice parametersfrom PBEsol+ U eff calculations are given as ( a , b , c , β ) = (6 . ˚ A , . ˚ A , . ˚ A , . ◦ ) . Internal coordinates optimized fromeDMFT calculations, with different choices of Coulomb interactions (Ising vs. beyond-Ising), are listed below. the beyond-Ising Coulomb terms. As a partial check, we havedone two eDMFT calculations for MnPS , employing twolattice parameter sets ( a , b , c , and the monoclinic angle β )corresponding to non-dimerized metallic and dimerized Mott-insulating states around GPa. Initial internal coordinates,optimized within eDMFT afterwards, were adopted from theambient pressure structure. The purpose of this comparisonis to see whether the choice of different Coulomb interac-tions yields noticeable difference. Table II summarizes theoptimized structures with Ising and beyond-Ising Coulomb in-teractions, which shows negligible difference with respect toeach other. Hence we are safe to use Ising form in this case.
IN-PLANE AND OUT-OF-PLANE HOPPING AMPLITUDES
Fig. S9(a) and (b) show root mean square (RMS) ampli-tudes of nearest-neighbor (n.n.) in-plane t and out-of-plane e g hopping integrals, respectively, for MnPS and NiPS .The d -orbital hopping integrals were computed using WAN - NIER
90 package[55], employing optimized crystal structuresin the presence of external pressure, without including U eff and magnetism.As shown in Fig. 1(a) and (b) in the main text, the in-plane lattice parameters (with respect to their ambient pres-sure value) a/a and b/b for MnPS around 30 GPa aresmaller by ∼
2% compared to those of NiPS , while the out-of-plane c/c of MnPS is larger than that of NiPS . In ac-cordance with the tendency of lattice parameter changes, theenhancement of RMS in-plane t hopping integrals is morepronounced in MnPS , which drives the formation of in-planeMn dimer formation after the transition to the low-spin statewith the open t shell. Note that, the in-plane t hoppingintegrals for Ni is also enhanced as the pressure is increased,but its effect is not significant due to the closed t shell in the Ni d configuration. Similarly, the enhanced out-of-planekinetic energy between the e g orbitals in NiPS , depicted inFig. S9(b), induces more pronounced reduction of the c pa-rameter in NiPS compared to that of MnPS . It should bementioned that, while the out-of-plane e g hopping terms arealso enhanced in MnPS , their role in structural response topressure is less significant both in the low- and high-pressureregimes; in the low-pressure regime the hybridization betweenthe d high-spin Mn ion and anions is weak, and so is theelectron-lattice coupling, while in the high-pressure regimethe high-spin Mn has empty e g shell. ELECTRONIC STRUCTURES WITH A N ´EEL ORDERAmbient pressure cases for both compounds
Fig. S10 shows the spectral functions of NiPS (a-c) andMnPS (d-f) with a N´eel-type antiferromagnetic (AF) order incomparison with the paramagnetic (PM) phases. In both com-pounds, the presence of magnetism does not alter the size ofcharge gap, consistent with the Mott character of the insulat-ing phases in both compounds. As the temperature is loweredand magnetism sets in, the broadening of spectral functionsdue to the imaginary part of self-energy is weakened. Indeed,DFT+ U eff PDOS with the magnetic order shows very similarqualitative features with DMFT PDOS at T =58K (not shown).Note that, the use of the Ising-type Coulomb repulsion givesrise to the stabilization of magnetic order well above the N´eeltemperatures of both compounds, T N = 154 and 78 K forNiPS and MnPS , respectively[25–27], as reported in pre-vious DFT+DMFT studies[56]. It is argued that a larger in-plane kinetic energy scale originating from the e g orbital inNi yields the higher T N in NiPS compared to MnPS [21].Despite the overestimated T N , such tendency can be noticed2 Pressure (GPa) V / V LS HS (a) DFT+U LSDFT+U HSPM-DMFT
LS HS
Str. interpolation (dimensionless)
LS HS
Str. interpolation (dimensionless) ∆ E ( e V ) M , | S | ( μ B ) (b)(c) ∆ E ( e V ) M , | S | ( μ B ) (d)(e) FIG. S8. (a) A magnified view of the MnPS volume vs. pres-sure plot, where the black solid and gray dashed lines are curvesfor ground and metastable states respectively. Two dash-dotted linesare at V = 0.57 and 0.54 V on which structural interpolations be-tween the high-spin honeycomb and low-spin dimerized structuresare made. (b,c) Total energy differences (b) and size of Mn mo-ments (c) as a function of structural interpolation at V = 0.57 V ,where red (blue) curve depicting total energy and Mn magnetizationstarting from high-spin (low-spin) structure and approaching to thelow-spin (high-spin) side, and purple symbols depicting same resultsfrom paramagnetic eDMFT calculations at T = 580K. (d,e) Sameplots at V = 0.54 V . Size of QMC error bars for eDMFT results are4 meV at most, smaller then the symbol size. in our results by comparing Fig. S10(b) and (e). While the sizeof the Ni magnetization in NiPS (2 S Ni = 1.45 µ B ) is almostsaturated to the value of PM moment ( | S Ni | = 1.65 µ B ), evenat T = 232K, the Mn magnetization in MnPS at the same T is 2 S Mn = 0.87 µ B , just a fraction of the S = 5/2 moment size(4.80 µ B ) of the high-spin Mn. As T is lowered to 58K, mag-netizations in both compounds saturate to the local momentsize, as shown in Fig. S10(c) and (f). Near the MIT critical pressure
MnPS As discussed in the main text, the pressure-induced MIT inMnPS in the paramagnetic phase is a discontinuous transi-tion accompanied by a spin-state transition and isosymmetricstructural distortion with a volume collapse. Such discontinu- ,in-plane t MnPS ,in-plane t ,out-of-plane e g MnPS ,out-of-plane e g R M S n . n . h o pp i n g a m p li t u d e ( e V ) Pressure (GPa) (a)(b)
FIG. S9. Root mean square (RMS) of nearest-neighbor (n.n.) Wan-nier hopping amplitudes for Ni and Mn d -orbitals as a function ofpressure, where the RMS of (a) in-plane t - t and (b) out-of-plane e g - e g are shown. ous character does not change in the presence of magnetism,as shown in Fig. S12, where the change of the Mn magnetiza-tion M (per Mn) from DFT+ U eff and eDMFT results (at T =232 K) are shown as a function of pressure. Note that, the up-turn of M in the small-pressure regime ( <
10 GPa) in eDMFTis due to the enhancement of magnetic exchange interactionsoriginating from increased kinetic energy scale under the pres-sure.
NiPS In the PM phase of NiPS , the MIT can be indicated notonly by the gap opening in the PDOS, but also by the changeof electron self-energies Σ( ω ) . In Fig. S11(a) and (c), the MITcan be noticed by the presence and absence of the dip at theFermi level in their PDOS, but it is slightly unclear whetherthe phase at P = 30.4 GPa is an insulator due to the smallbut finite e g DOS at the Fermi level due to the broadening.However, (cid:61) Σ( ω ) , plotted in Fig. S11(b) and (d), show a cleardifference between the two phases, because the presence (ab-sence) of a pole at the Fermi level in the e g - (cid:61) Σ( ω ) signifiesthe presence (absence) of the Mott physics.In the AF-ordered phases, the gap opening is induced bythe exchange splitting between the spin up and down com-ponents, i.e. (cid:60) (Σ up − Σ down ) ( ω ) . In cases where (cid:61) Σ( ω ) is weak compared to (cid:60) Σ( ω ) and the frequency dependenceof (cid:60) (Σ up − Σ down ) ( ω ) is small, then the eDMFT results be-come equivalent to the DFT+ U eff results. Fig. S11(e-j) presentsuch situation, where the PDOSs shown in Fig. S11(e) and3 P D O S ( s t a t e s / e V / N i ) PM, T=232K AF, T=232K AF, T=58K
Energy (eV) Mn PS (d)(e)(f) |2S Mn | = 4.80 μ Β |2S Mn | = 4.80 μ Β Mn = 0.87 μ Β |2S Mn | = 4.80 μ Β Mn = 4.80 μ Β g ’a e g -1-0.5 0 0.5 1-2-1 0 1 2 -4 -2 0 2 4 AF, T=232K AF, T=58K
Energy (eV)
PM, T=232K Ni PS (a)(b)(c) g ’a e g -4-2 0 2 4-4-2 0 2 4 -2 -1 0 1 2 |2S Ni | = 1.65 μ Β |2S Ni | = 1.65 μ Β Ni = 1.45 μ Β |2S Ni | = 1.64 μ Β Ni = 1.62 μ Β Fig. 2
FIG. S10. (a-c) PDOS of NiPS and (d-f) MnPS with a N´eel type antiferromagnetic (AF) order (b,c,e,f) in comparison with paramagnetic(PM) PDOS (a,d). The second and third rows show AF PDOS with T = 232K and 58K, respectively. (h) become qualitatively equivalent to DFT+ U eff PDOS (notshown), with the exchange splitting of ∼ e g bands. Hence, in AF phases theMIT critical pressure is mainly determined by the e g band-width and the exchange splitting (cid:60) (Σ up − Σ down ) ( ω ) . Be-cause the above quantities change continuously as the pres-sure is increased, it is not easy to point out at which pressurethe MIT happens from the PDOS plots due to the presence ofsmall broadening from (cid:61) Σ( ω ) . At T = 232K, the MIT seemsto happen around 24 GPa, and this pressure does not changeas T is lowered to 116K.Fig. S13(a) depicts the change of Ni magnetizations fromDFT+ U eff and AF-eDMFT as a function of pressure. Notethat, the MIT critical pressures are 36 and 24 ∼
30 GPa forDFT+ U eff and AF-eDMFT results, as shown in the figure,but the magnetization persists within the metallic phase. Thepressure where the magnetism disappears in AF-eDMFT re-sults increases slightly from 36 to 40 GPa as the T is reducedfrom 232K to 116K, but it does not reach 58 GPa where themagnetization disappears in DFT+ U eff results. Note that, thehigh-pressure anisotropic structural distortion in NiPS hap-pens at the pressure where the DFT+ U eff magnetization be-comes zero.Fig. S13(b) shows the change of Ni local spin moment size | S | as a function of pressure. Unlike the | S | in the PMphase, which shows a cusp at the MIT critical pressure, | S | in the AF phases does not show such behavior at the MIT pres-sure. As the pressure is increased, the AF | S | is suppresseduntil the magnetization disappears and a PM metallic phasehappens. Note that, the pressures that AF | S | joins the PM | S | curve are the points where the magnetization disappears.4 g ’a e g -4 -2 0 2 4 -4 -2 0 2 4 - I m Σ ( ω ) ( e V ) - I m Σ ( ω ) ( e V ) P D O S ( / e V ) P D O S ( / e V ) Energy (eV)(a)(b)(c)(d)
P = 30.4 GPa30.4 GPaP = 31.8 GPa31.8 GPa -1 0 1 e g ’a e g -202048 -4 -2 0 2 4 -4 -2 0 2 4 NiPS , PM , T = 232 K P = 30.4 vs. 31.8 GPa NiPS , AF , T = 232 K P = 24.6 vs. 31.8 GPa P D O S ( / e V ) - I m Σ u p / d o w n ( ω ) ( e V ) R e ( Σ u p - Σ d o w n )( ω ) ( e V ) Energy (eV) Energy (eV)(e) (f)(g) (h) (i)(j)P = 24.6 GPa P = 31.8 GPa
Spin UpSpin Down -Im Σ up ( ω )-Im Σ down ( ω ) FIG. S11. PDOS and electron self-energies Σ( ω ) of (a-d) PM and (e-j) AF NiPS at T = 232K. The imaginary part of PM self-energiesin (b) and (d) show a clear e g peak and its absence in the Mott-insulating and metallic phases, respectively. (e) and (h) show spin- andorbital-projected PDOS at P = 24.6 and 31.8 GPa, respectively (spin up and down PDOS plotted in positive and negative values respectively).Imaginary part of Σ( ω ) , plotted in (f) and (i), show the absence of peak near the Fermi level. (e) and (h) show the difference between the realpart of spin up and down self-energies (cid:60) (Σ up − Σ down ) ( ω ) , which is the on-site exchange splitting opening the gap in the magnetic phases.Note that, the MIT point is not clear due to the smearing of spin up and down e g bands in the PDOS. M , | S z | ( μ Β ) Pressure (GPa)
DMFT |2S z |DMFT MDFT+U M FIG. S12. Evolution of sizes of DFT+ U eff and eDMFT local mo-ments | S z | and Mn magnetization M (per Mn) in MnPS as a func-tion of pressure. eDMFT results are obtained at T = 232 K with theN´eel order. PM-DMFT MIT P c DFT+ U MIT P c A F - D M FT M I T P c Pressure (GPa) M ( μ Β ) | S z | ( μ Β ) PM-DMFT
MIT P c A F - D M FT M I T P c (a)(b) M (DFT+U)M (AF-DMFT,T=232K)M (AF-DMFT,T=116K) |2S z | (PM-DMFT,T=232K)|2S z | (AF-DMFT,T=232K)|2S z | (AF-DMFT,T=116K) FIG. S13. (a) Evolution of Ni magnetization M (per Ni) fromDFT+ U eff (gray curve), AF-eDMFT at T = 232 (red) and 58K (blue)as a function of pressure. The MIT critical pressures from DFT+ U eff ,PM- and AF-eDMFT are depicted. (b) Sizes of Ni spin moments ||
MIT P c A F - D M FT M I T P c (a)(b) M (DFT+U)M (AF-DMFT,T=232K)M (AF-DMFT,T=116K) |2S z | (PM-DMFT,T=232K)|2S z | (AF-DMFT,T=232K)|2S z | (AF-DMFT,T=116K) FIG. S13. (a) Evolution of Ni magnetization M (per Ni) fromDFT+ U eff (gray curve), AF-eDMFT at T = 232 (red) and 58K (blue)as a function of pressure. The MIT critical pressures from DFT+ U eff ,PM- and AF-eDMFT are depicted. (b) Sizes of Ni spin moments || S ||