Multiple Dirac Cones and Topological Magnetism in Honeycomb-Monolayer Transition Metal Trichalcogenides

TL;DR

This paper predicts that monolayer transition metal trichalcogenides with honeycomb structures can host multiple Dirac cones, exhibit tunable gaps, and become topological ferromagnets with high Chern numbers due to electron correlations and doping.

Contribution

It introduces a new class of 2D materials with multiple Dirac cones and topological magnetism, expanding the understanding of correlated Dirac electrons in honeycomb monolayers.

Findings

Multiple Dirac cones with tunable gaps identified.

Transition to topological ferromagnet with high Chern number.

Potential for new atomically thin heterostructures.

Abstract

The discovery of monolayer graphene has initiated two fertile fields in modern condensed matter physics, Dirac semimetals and atomically-thin layered materials. When these trends meet again in transition metal compounds, which possess spin and orbital degrees of freedom and strong electron correlations, more exotic phenomena are expected to emerge in the cross section of topological states of matter and Mott physics. Here, we show by using ab initio calculations that a monolayer form of transition metal trichalcogenides (TMTs), which has a honeycomb network of transition metal cations, may exhibit multiple Dirac cones with tunable gaps in the electronic structure. Furthermore, we elucidate that electron correlations and carrier doping turn the multiple-Dirac semimetal into a topological ferromagnet with high Chern number. Our findings raise the honeycomb-monolayer TMTs to a new paradigm to explore correlated Dirac electrons and topologically-nontrivial magnetism. In turn, the unique wide-ranging properties of the materials will deliver new building blocks for atomically thin heterostructures.

Tables

Table 1. Table 1: Transfer integrals between MLWFs. Each value in the table means ⟨ m , 𝟎 | ​ H ​ | n , 𝒓 ⟩ bra 𝑚 0 𝐻 ket 𝑛 𝒓 \bra{m,\bm{0}}H\ket{n,\bm{r}} , where H 𝐻 H is the Hamiltonian of the system and | m , 𝒓 ⟩ ket 𝑚 𝒓 \ket{m,\bm{r}} is the d m subscript 𝑑 𝑚 d_{m} -like MLWF at site 𝒓 𝒓 \bm{r} ( m 𝑚 m = 3 ​ z 2 − r 2 3 superscript 𝑧 2 superscript 𝑟 2 3z^{2}-r^{2} or x 2 − y 2 superscript 𝑥 2 superscript 𝑦 2 x^{2}-y^{2} ). We take 𝒓 = 𝑹 i 𝒓 subscript 𝑹 𝑖 \bm{r}=\bm{R}_{i} ( i 𝑖 i = 1, 2, 3, 4, or 5) illustrated in Fig.  2 (c). The unit of transfer integrals is in meV.

$(m,n)$	${𝑹}_1$	${𝑹}_2$	${𝑹}_3$	${𝑹}_4$	${𝑹}_5$
$(3z^2-r^2,3z^2-r^2)$	$-87$	$-9$	$-38$	$4$	$-12$
$(3z^2-r^2,x^2-y^2)$	$0$	$22$	$0$	$-8$	$0$
$(x^2-y^2,3z^2-r^2)$	$0$	$18$	$0$	$-8$	$0$
$(x^2-y^2,x^2-y^2)$	$-70$	$14$	$304$	$7$	$30$

Table 2. Table 2: The lattice constant and the layer thickness of optimized structures of M 𝑀 M P X 3 subscript 𝑋 3 X_{3} . Schematic pictures of a 𝑎 a and n 𝑛 n are illustrated in Fig.  5 . In the honeycomb network of the edge-sharing ideal octahedra, the ratio of the parameters should be a / n = 3 / 2 ∼ 2.12 𝑎 𝑛 3 2 similar-to 2.12 a/n=3/\sqrt{2}\sim 2.12 .

	NiPS3	PdPS3	PtPS3	NiPSe3	PdPSe3	PtPSe3
$a$ (Å)	$5.82$	$6.02$	$6.07$	$6.17$	$6.34$	$6.42$
$n$ (Å)	$2.98$	$3.17$	$3.07$	$3.07$	$3.27$	$3.12$
$a/n$	$1.95$	$1.90$	$1.97$	$2.01$	$1.94$	$2.05$

Table 3. Table 3: Transfer integrals between MLWFs. Each value in the table means ⟨ m , 𝟎 | ​ H ​ | n , 𝒓 ⟩ bra 𝑚 0 𝐻 ket 𝑛 𝒓 \bra{m,\bm{0}}H\ket{n,\bm{r}} , where H 𝐻 H is the Hamiltonian of the system and | m , 𝒓 ⟩ ket 𝑚 𝒓 \ket{m,\bm{r}} is the d m subscript 𝑑 𝑚 d_{m} -like MLWF at site 𝒓 𝒓 \bm{r} ( m 𝑚 m = 3 ​ z 2 − r 2 3 superscript 𝑧 2 superscript 𝑟 2 3z^{2}-r^{2} or x 2 − y 2 superscript 𝑥 2 superscript 𝑦 2 x^{2}-y^{2} ). See Fig.  5 (a) for the spatial positions of 𝒓 𝒓 \bm{r} . The unit of transfer integrals is in meV.

$(m,n)$	${𝑹}_1$	${𝑹}_1^1$	${𝑹}_1^2$	${𝑹}_2$	${𝑹}_2^1$	${𝑹}_2^2$	${𝑹}_2^3$	${𝑹}_2^4$	${𝑹}_2^5$
$(3z^2-r^2,3z^2-r^2)$	$-87$	$-72$	$-72$	$-9$	$-9$	$26$	$-9$	$-9$	$26$
$(3z^2-r^2,x^2-y^2)$	$0$	$9$	$-9$	$22$	$-22$	$2$	$18$	$-18$	$-2$
$(x^2-y^2,3z^2-r^2)$	$0$	$9$	$-9$	$18$	$-18$	$-2$	$22$	$-22$	$2$
$(x^2-y^2,x^2-y^2)$	$-70$	$-82$	$-82$	$14$	$14$	$-21$	$14$	$14$	$-21$

Table 4. Table 4: The total energy of each electronic state obtained by GGA calculations. The lowest energy of all electronic states is set to be zero for each compound. The blanks indicate that the corresponding state is not obtained as a stable solution. The unit of energy is in meV per rectangular unit cell including four formula units (see Fig.  8 ). For the magnetic solutions, we denote the value of the magnetic moment per transition metal cation in unit of Bohr magneton μ B subscript 𝜇 B \mu_{\rm B} .

	PM	FM	AFM	zAFM	sAFM
NiPS3	$914$	$444\left(1.24\right)$	$41\left(1.13\right)$	$0\left(1.12\right)$	$491\left(1.24\right)$
PdPS3	$5$	–	$0\left(0.26\right)$	$1\left(0.27\right)$	–
PtPS3	$0$	–	–	–	–

Equations

$$H_{\mathrm{SOC}}=-(\tilde{\lambda }/2)\sum _{\mathbf{k}}\sum _{mn}\sum _{\sigma {\sigma }^{\prime }}{c}_{\mathbf{k}m\sigma }^{†}({\hat{\tau }}_y{)}_{mn}({\hat{\sigma }}_z{)}_{\sigma {\sigma }^{\prime }}{c}_{\mathbf{k}n{\sigma }^{\prime }},$$

$$H_{\mathrm{int}}=\frac{1}{2}\sum _{mn{m}^{\prime }{n}^{\prime }}{U}_{mn{m}^{\prime }{n}^{\prime }}\sum _i\sum _{\sigma {\sigma }^{\prime }}{c}_{im\sigma }^{†}{c}_{in{\sigma }^{\prime }}^{†}{c}_{i{n}^{\prime }{\sigma }^{\prime }}{c}_{i{m}^{\prime }\sigma },$$

Full text

I Introduction

Since the success of exfoliation of a monolayer graphene Novoselov666 , atomically-thin layered materials have grown as one of the leading themes in modern condensed matter physics. In particular, van der Waals (vdW) materials, composed of atomic layers bounded via weak vdW forces, have received great attention. Electrons confined in an atomically thin layer exhibit drastically distinct behavior from the bulk form. The archetypal example is the Dirac electrons in a monolayer graphene, which show anomalous transport behavior, e.g., the anomalous integer quantum Hall effect novoselov2005two ; zhang2005experimental and the Klein tunneling katsnelson2006chiral ; young2009quantum . Another example is the valley degree of freedom in the monolayer form of transition metal dichalcogenides doi:10.1021/nl903868w , which has been intensively studied toward valleytoronics devices wang2012electronics ; xu2014spin . Furthermore, heterostructures of different vdW materials have provided a new platform for novel functionalities never seen in bulk compounds geim2013van ; novoselov20162d .

Through the intensive research in the past decade, a lot of efforts have been made to find atomically-thin magnetic materials. Among many candidates, a family of transition metal trichalcogenides (TMTs) has gained increasing interests, both from theoretical proposals of monolayer magnetism PhysRevB.91.235425 ; PhysRevB.94.184428 and experimental reports on the mono and few-layer forms doi:10.1021/acsnano.5b05927 ; C5TC03463A ; kuo2016exfoliation ; doi:10.1021/acs.nanolett.6b03052 ; 2053-1583-3-3-031009 ; Gong2017 . In addition, not only the magnetism but also anomalous electronic and transport properties are predicted in the presence of the relativistic spin-orbit coupling (SOC), e.g., the spin-valley coupling Li05032013 , the magnon spin Nernst effect PhysRevLett.117.217202 , and the gate-controllable magneto-optic Kerr effect PhysRevLett.117.267203 . Thus, the atomically-thin layered TMTs are expected to provide a unique cross section between strong electron correlations and the SOC, but their potential remains unexplored.

In this paper, we theoretically propose that monolayer TMTs with a honeycomb network of $4d$ and $5d$ transition metals would host a new playground for correlated Dirac electrons and topologically-nontrivial magnetism. By ab initio calculations, we show that the TMTs with group 10 transition elements have semimetallic band structures with multiple Dirac cones in the half-filled $e_{g}$ orbitals in the paramagnetic state. We find that the multiple Dirac cones originate in electron transfers on a hidden honeycomb superstructure emergent from spatially anisotropic $d$ orbitals with large hybridization with the neighboring ligand $p$ orbitals. We also show that the SOC gaps out these Dirac cones in the presence of the trigonal lattice distortions, and hence, the mass gap can be flexibly tuned by the tensile strain. In addition, by the mean-field analysis for an effective multi-orbital Hubbard model, we elucidate that electron correlations and chemical doping potentially change the multiple-Dirac semimetals into a topological ferromagnet with high Chern numbers.

This paper is organized as follows. In Sec. II, we introduce the method of ab initio calculations used in this paper. In Sec. III, we show the results of the ab initio calculations and the mean-field analysis of an effective multi-orbital Hubbard model for monolayer TMT. Section IV is devoted to the summary and the discussion of future issues. In Appendix A, we show the details of optimized crystalline structures, ab initio band structures, and transfer integrals of monolayer TMTs. In Appendix B, we provide the information on distortions under tensile strain. The ab initio results for magnetism of monolayer TMTs are presented in Appendix C.

II

ab initio calculation

The chemical formula for TMTs is generally given by $MBX_{3}$, where $M$ is transition metals, $B$=P, Si, or Ge, and $X$ is chalcogens. TMTs have vdW layered structures, whose stacking manner depends on the compounds ZAAC:ZAAC19733960305 ; BREC19863 ; OUVRARD198827 ; carteaux1995crystallographic . In each layer, transition metal cations $M$ comprise a honeycomb network by sharing the edges of $MX_{6}$ octahedra, and $B_{2}$ dimers locate at the centers of the hexagons of the honeycomb network [Fig. 1(a)].

In this paper, we focus on a monolayer form of TMTs with $B$=P and $X$=S and Se. In this case, the nominal valence of the transition metal cation is $M^{2+}$. We note that all 3$d$, 4$d$, and 5$d$ transition metal elements belonging to group 10 and 12 can take the stable divalent oxidation state tagkey1997ii . Indeed, $M$P$X_{3}$ with $M$=Ni, Pd, Zn, Cd, and Hg have been synthesized ZAAC:ZAAC19733960305 . In the following, we consider group 10 elements, $M$=Ni, Pd, and Pt.

We calculate the electronic band structures of monolayer $M$P$X_{3}$ by ab initio calculations based on the generalized gradient approximation (GGA). In the calculations, we used OpenMX code openmx , which is based on a linear combination of pseudoatomic orbital formalism PhysRevB.67.155108 ; PhysRevB.69.195113 . We adopted the Perdew-Burke-Ernzerhof GGA functional in density functional theory PhysRevLett.77.3865 and a $30\times 30\times 1$ $\bm{k}$-point mesh for the calculations of the self-consistent electron density and the structure relaxation. We inserted vacuum space greater than 10Å between monolayers and fully relaxed the primitive vectors and atomic positions in the unit cell with the convergence criterion 0.01 eV/Å about the inter-atomic forces.

For all combinations of $M$=Ni, Pd, and Pt and $X$=S and Se, we performed full structural optimization in the paramagnetic state without the SOC, starting from the reported crystalline data of NiPS3 or PdPS3 ZAAC:ZAAC19733960305 . We confirmed that all the cases stably converge on the similar structure, as shown in Appendix A. The SOC is incorporated by the relativistic ab initio calculations for the optimized structures. In the following sections, we will discuss the material trend of $M$P$X_{3}$ on the basis of PdPS3, as it represents the typical band structure of $M$P$X_{3}$ and locates between weakly correlated $5d$ and strongly correlated $3d$ systems. (Indeed, we will discuss that the compound might be close to the border between paramagnetic and antiferromagnetic phases.) For the microscopic analyses based on the tight-binding models, we construct the maximally localized Wannier functions (MLWFs) PhysRevB.56.12847 ; PhysRevB.65.035109 for the $e_{g}$ bands and evaluate the transfer integrals between them via a code implemented in OpenMX openmx . We also compute the magnetic solutions by using GGA scheme without the SOC under the full structural optimization, as shown in Appendix C.

III Results

III.1 Multiple Dirac cones

We show the electronic band structures in the paramagnetic state obtained by ab initio calculations without the SOC. Figures 1(b) and 1(c) show the representative results for PdPS3. The Pd $d$-orbital levels are split into two groups, $e_{g}$ and $t_{2g}$, due to the crystalline electric fields of octahedral ligands. As Pd2+ is in the $d^{8}$ electron configuration, the lower-energy $t_{2g}$ manifold is fully occupied and the higher-energy $e_{g}$ manifold is half-filled. Remarkably, the $e_{g}$ bands have two crossing points at the K point and around the midpoint on the $\Gamma$-K line in the Brillouin zone [Fig. 1(b)], and the projected density of states are almost zero at the Fermi level [Fig. 1(c)]. We find that the crossings are the Dirac cones, as shown in Fig. 1(d): two electronic bands near the Fermi level give rise to eight Dirac cones (two on the zone boundary, K and K’, and other six inside). We confirm that the multiple Dirac cones are shared by other monolayer TMTs with $M$=Ni, Pd, and Pt (see Appendix A for details).

In order to clarify the microscopic origin of the multiple Dirac cones, we show the MLWFs obtained from two initial states, $d_{3z^{2}-r^{2}}$ and $d_{x^{2}-y^{2}}$, in Figs. 2(a) and 2(b), respectively. Both MLWFs well extend over the neighboring S sites, indicating the importance of indirect hopping processes via the ligand $p$ orbitals. Table 1 shows the representative transfer integrals between the two types of MLWFs for the Pd-Pd bonds up to fifth neighbors [see Fig. 2(c)]. (See Appendix A for the extended list of transfer integrals.) We construct a tight-binding model by using these transfer integrals and confirm that the model well reproduces the ab initio band structure [see Fig. 1(b)]. Interestingly, the most dominant electron transfer is not for nearest neighbors but the third neighbors. This is understood from the fact that the indirect $d$-$p$-$d$ hoppings between nearest neighbors are almost forbidden, while the $d$-$p$-$p$-$d$ ones for third neighbors are substantial, as shown in Figs. 2(d) and 2(e). We note that a similar argument was made for magnetic exchange interactions PhysRevB.91.235425 .

The dominant third neighbor transfers explain the origin of the multiple Dirac cones. As well-known in graphene, the nearest neighbor transfers produce the Dirac cones at the zone corners, the K and K’ points. This is also the case for the $e_{g}$ electron systems xiao2011interface . On the other hand, the third neighbor transfers bring about Dirac cones at the additional six points inside the first Brillouin zone: the network of the third neighbor bonds forms honeycomb superstructures with the lattice spacing twice longer than the original honeycomb network, as exemplified in Fig. 2(c), which leads to new Dirac cones around the midpoints of the $\Gamma$-K lines (zone corners in the folded Brillouin zone). Thus, in our TMTs, the hidden honeycomb superstructures stemming from the orbital and geometric nature result in the multiple Dirac nodes.

III.2 Tunable Dirac gap

Next, we discuss the effect of SOC. Although the orbital moment is quenched in the $e_{g}$ manifold in an ideal octahedral crystal field, the SOC modifies the $e_{g}$ electronic states through $t_{2g}$-$e_{g}$ mixing in the presence of a distortion of $MX_{6}$ octahedra. Indeed, we find that the Dirac nodes are gapped out by including the relativistic effect in the ab initio calculations, as shown in Fig. 3(a). In this monolayer system, the dominant distortion is a trigonal one, which leads to an effective SOC given as xiao2011interface

[TABLE]

where $c^{\dagger}_{\bm{k}m\sigma}$($c_{\bm{k}m\sigma}$) is the creation (annihilation) operator of an electron for the wave vector $\bm{k}$, orbital $m=d_{3z^{2}-r^{2}}$ or $d_{x^{2}-y^{2}}$, spin $\sigma=\uparrow$ or $\downarrow$, and $\tau_{y}$ ($\sigma_{z}$) is the $y$ ($z$) component of the Pauli matrix for the orbital (spin) space; here, the $xyz$-axes are taken as shown in the inset of Fig. 2(a) and the quantization axis of spin is taken along the [111] direction. The coupling constant is given as $\tilde{\lambda}=\Delta_{\rm tri}\lambda^{2}/\Delta^{2}$, where $\Delta$ and $\Delta_{\rm tri}$ are the crystalline electric field from the octahedral ligands and the trigonal distortion, respectively, and $\lambda$ is the coupling constant of the atomic SOC.

Indeed, we confirm that when adding the effective SOC in Eq. (1) to the tight-binding model constructed above, its band structures reproduce the gapped Dirac nodes, as shown in Fig. 3(a). From the comparison, we obtain the rough estimate of the effective SOC $\tilde{\lambda}=15$ meV. We note that the gapped state is topologically trivial: the $Z_{2}$ topological invariant PhysRevLett.95.146802 becomes zero for all the bands in the tight-binding model.

The result indicates that the Dirac gaps can be controlled through the crystalline symmetry. Here, we demonstrate it by tensile strain, which has been commonly used for two-dimensional vdW materials Lee385 ; guinea2010energy . Starting from the fully-optimized crystalline structure at zero expansive ratio, we extend two in-plane primitive vectors, $\bm{a}$ and $\bm{b}$, while keeping the out-of-plane primitive vector $\bm{c}$. We also keep the fractional coordinates of atoms projected onto the $ab$ plane. Figure 3(b) shows the change in the Dirac gaps. In the original optimal structure (zero expansive ratio), the octahedra are slightly elongated in the out-of-plane direction. While the system is expanded in the in-plane directions, the Dirac gaps decrease and become minimal around $8$-$9$% expansive ratio, where the trigonal distortion almost vanishes (see Appendix B for the information on distortions of the octahedra). Interestingly, the valley structures of the two massive Dirac cones are shifted individually by the tensile strain, as shown in Fig. 3(c). These results indicate the flexible tunability of the massive Dirac cones.

III.3 Topological magnetism induced by electron correlations

Let us discuss the effect of electron correlations in $d$ electrons on the multiple-Dirac semimetals. Although the magnetism was studied for monolayer TMTs by ab initio calculations PhysRevB.91.235425 ; PhysRevB.94.184428 , the previous works focused on the $3d$ compounds in which the SOC is irrelevant. We here investigate the synergetic effect of electron correlations and the SOC, both of which can be relevant in $4d$ and $5d$ compounds, using an effective multi-orbital Hubbard model and the mean-filed approximation. The multi-orbital Hubbard model consists of the tight-binding model with the transfer integrals in Table 1, the effective SOC in Eq. (1) with $\tilde{\lambda}=$15 meV, and the onsite Coulomb interactions given by

[TABLE]

where $c^{\dagger}_{im\sigma}$($c_{im\sigma}$) is the Fourier transform of $c^{\dagger}_{\bm{k}m\sigma}$($c_{\bm{k}m\sigma}$). Assuming the rotational symmetry of the Coulomb interaction, we set $U_{mmmm}=U$, $U_{mnmn}=U-2J$, and $U_{mnnm}=U_{mmnn}=J$ ($m\neq n$), where $U$ is the intraorbital Coulomb interaction and $J$ is the Hund’s coupling, respectively; we take $J/U=0.2$ in the following calculations. In the mean-field calculation, we adopt the standard Hartree-Fock approximation to decouple the onsite interaction terms in Eq. (2). We take into account charge, spin, and orbital orders with the ordering vector $\bm{Q}=(0,0)$ or $(\pi,\pi)$ on the honeycomb lattice and approximate the integration in the first Brillouin zone by the summation over 128$\times$128 $\bm{k}$ points and determine the mean fields consistently within a precision of less than 10*-6*.

We focus on two commensurate fillings, half filling (two $e_{g}$ electrons per $M$) and 3/4 filling (three $e_{g}$ electrons per $M$); the former corresponds to the situation discussed above, while the latter a chemical substitution of $M$ by, e.g., Ag or Cd. We note that such substitutions were reports for the bulk form of TMTs BREC19863 . Figures 4(a) and 4(b) show the ground-state phase diagrams and the magnetic moments, obtained by the mean-field analysis in the range of Coulomb interactions including realistic values for 4$d$ and 5$d$ transition metal compounds 0953-8984-29-26-263001 .

At half filling [Fig. 4(a)], while increasing the electron interactions, the system exhibits a continuous phase transition from the paramagnetic Dirac semimetal to a Néel-type antiferromagnetic insulator (AFMI) with in-plane magnetic moments. For comparison, we also performed the ab initio calculations in the GGA scheme with allowing magnetic solutions (see Appendix C for the details). We find that, similar to the mean-field results, the lowest-energy state changes from the paramagnetic Dirac semimetal to AFMI while changing from weakly correlated $M$=Pt to strongly correlated $M$=Ni; the $M$=Pd case is close to the boarder. We note that, in general, the ab initio calculation tends to underestimate the correlation effects while the mean-field approximation tends to overestimate. From these considerations, we conclude that PdPS3 might be in the AFM phase, while PtPS3 the multiple-Dirac semimetal (see also the discussion in Appendix C). As PdPS3 appears to locate close to the border, it might also be possible to transform it to the multiple-Dirac semimetal by tuning the bandwidth by the substitution of S by Se (see Appendix A) or by tensile strain.

On the other hand, at 3/4 filling, the system shows a discontinuous phase transition from the paramagnetic metal to a ferromagnetic metal (FMM), and to a ferromagnetic insulator (FMI) with out-of-plane magnetic moments, as shown in Fig. 4(b). For comparison, we also performed the GGA calculation for AgPS3, which reproduces the 3/4-filled state, using the crystalline structure of PdPS3 ZAAC:ZAAC19733960305 . We found that the lowest-energy state is the FM state with a small magnetic moment, 0.04 $\mu_{\rm B}$, which might also reflect the general tendency of GGA calculations to underestimate electron interactions.

Interestingly, we find that the ferromagnetic states at 3/4 filling acquires nontrivial topological nature. Figure 4(c) shows the band structure of the FMI at $U=1.5$ eV. The bands are split by the exchange field into the up-spin (red) and down-spin (blue) ones, and the lower six are occupied at 3/4 filling. Note that as the mean-field Hamiltonian in the ferromagnetic state conserves the spin $z$ component, we can separate the mean-field Hamiltonian into up-spin and down-spin sectors and distinguish the spin state of each electronic band in Fig. 4(c). We calculate the Berry curvature and the Chern number for each band of the mean-field solution by using the standard Kubo formula PhysRevLett.49.405 . Summing the Chern number of the occupied bands [see Fig. 4(c)], we find that the FMI is a topologically-nontrivial ferromagnet with rather high Chern number $C=4$. Figure 4(d) displays the wave-number dependence of the Berry curvature of the highest-occupied band with $C=6$. The Berry curvatures shows spikes at the K and K’ points and around the midpoints of the $\Gamma$-K lines. These anomalous contributions can be traced back to the Dirac cones in the original semimetallic state. Thus, our results suggest that the multiple-Dirac semimetal can be turned into an unconventional topological ferromagnet with high Chern number by electron correlations and carrier doping.

IV Summary and concluding remarks

To summarize, we have theoretically uncovered two potential electronic properties of TMTs with 4$d$ and 5$d$ transition metals in the monolayer form. One is the highly-tunable multiple Dirac cones. This will bring about new transport phenomena, such as the unconventional Hall responses and the multiple valley operations. The other is the topological ferromagnetism with high Chern number driven by electron correlations and chemical doping. This will provide new candidates for quantized anomalous Hall insulators, whose multiple chiral edge modes might be used for a thin-film transmitter with high efficiency. We believe that the two features will stimulate further material exploration in 4$d$ and 5$d$ TMTs for delivering missing pieces in material science of atomically-thin films and the heterostructures.

The recipe for multiple Dirac cones found here is generic and simple: transition metal cations with $e_{g}$ orbitals and octahedral ligands forming edge-sharing honeycomb structure with extended $p$ orbitals. As such a crystalline and electronic structure is widely seen in layered transition metal compounds, e.g., transition metal chalcogenides and halides Hulliger1976 , it is of great interest to search other candidates of multiple-Dirac semimetals. Another intriguing issue would be the effect of electron correlations on the Dirac cones, which has been intensively studied in graphene and related compounds. The additional valley degrees of freedom in our multiple Dirac nodes may give rise to richer physics between the Dirac semimetal and the Mott insulator. Therefore, we believe that our findings stimulate the further exploration of TMTs and their relatives, which would pave a way to post-graphene nanotechnology.

Acknowledgements.

Y.S. is supported by the Japan Society for the Promotion of Science through a research fellowship for young scientists and the Program for Leading Graduate Schools (MERIT). This research was supported by Grants-in-Aid for Scientific Research under Grants No. JP15K05176. The crystal structures and MLWFs are visualized by using VESTA 3 Momma:db5098 .

Appendix A Ab initio results for monolayer $M$P$X_{3}$ in the paramagnetic state

In this Appendix, we show the details of ab initio results for $M$P$X_{3}$ ($M$=Ni, Pd, and Pt, and $X$=S and Se). We confirm that the honeycomb monolayer form (Fig. 5) is structurally stable for all $M$P$X_{3}$, though the lattice constant $a$ is different; see Table 2. We note that the optimized structures are elongated along the out-of-plane direction compared to the ideal octahedra: in all cases, the ratio $a/n$ (see Fig. 5) is smaller than the ideal value $3/\sqrt{2}\sim 2.12$.

Figure 6 shows the electronic band structures of monolayer $M$P$X_{3}$ obtained by ab initio calculations with and without the SOC. In all cases, the band structures show multiple Dirac dispersions at the middle points on the $\Gamma$-K lines as well as the K and K’ points near the Fermi level, which are gapped out by the SOC. This indicates that the multiple Dirac cones originating from the orbital and geometric nature discussed in Sec. III.1 are ubiquitous in the monolayer TMTs. We note that similar band dispersions were already seen in the previous ab initio studies PhysRevB.94.184428 although less attention has been paid. We also note that the previous experimental work reported insulating behavior even above the Néel temperature for NiPS3 in the bulk form ic50197a018 , while our result for the monolayer form is semimetallic in the paramagnetic state. The discrepancy would be ascribed to the difference between bulk and monolayer or the importance of electron interactions in 3$d$ electrons, which is in general underestimated in ab initio calculations.

In addition, we provide the detailed information on the transfer integrals estimated by MLWFs for PdPS3. While the representative values are shown in Table 2, we enlist all the values in Table 3 up to fifth neighbors. The spatial positions of the neighbors are illustrated in Fig. 5(a). As discussed in Sec. III.1, the transfer integrals for third neighbors are the most dominant among them on average. Note that the transfer integrals for different directions in the same distance are related with each other via the point-group operations; for instance, when the $MX_{6}$ octahedra have $D_{3h}$ symmetry, the hopping matrices should follow $\hat{h}(\hat{C_{3}}\bm{R})=\hat{R}_{C_{3}}\hat{h}(\bm{R})\hat{R}^{-1}_{C_{3}}$, where $(\hat{h}(\bm{R}))_{mn}=\bra{m,\bm{0}}H\ket{n,\bm{R}}$ and $\hat{R}_{C_{3}}$ is the threefold-rotation operator on the $e_{g}$-orbital basis. Indeed, Table 3 indicates that such relations hold approximately for the optimized structure in the ab initio calculation.

Appendix B Distortion of octahedra by the tensile strain

In this Appendix, we show how PdS6 octahedra are distorted for the tensile strain discussed in Sec. III.2. We adopt two conventional measures for the distortion robinson1971quadratic : one is the bond angle variance $\sum^{12}_{i=1}(\theta_{i}-90^{\circ})^{2}/11$, where $\theta_{i}$ is the angle of a S-Pd-S bond for neighboring S, and the other is the quadratic elongation $\sum^{6}_{i=1}(l_{i}/l_{0})^{2}/6$, where $l_{i}$ is the length of a Pd-S bond and $l_{0}$ is the bond length of an ideal octahedron with the same volume. Figure 7 shows the the bond angle variance and the quadratic elongation of PdS6 octahedra as functions of the expansive ratio obtained by our ab initio calculations. The result indicates that PdS6 octahedra is originally elongated along the $c$ axis at the zero expansive ratio (see also Table 2) and approaches to ideal octahedra while increasing the expansive ratio; the distortion is minimized around 10%. This is consistent with the behavior of the Dirac gaps plotted in Fig. 3(b), which is predominantly opened by the trigonal distortion of the octahedra.

Appendix C Total energy comparison by GGA calculation

In this Appendix, we show the ab initio results while allowing magnetic solutions. Using GGA calculations without the SOC under the full structural optimization, we compute the total energy of paramagnetic (PM), ferromagnetic (FM), Néel antiferromagnetic (AFM), zigzag antiferromagnetic (zAFM), and stripy antiferromagnetic (sAFM) states. We adopted a rectangular unit cell (see Fig. 8) and $16\times 16\times 1$ $\bm{k}$-point mesh.

Table 4 shows the total energy and the magnetic moment in each state. The lowest energy solution is for the zAFM for the $3d$ electron system with $M$=Ni, while it is PM for the $5d$ $M$=Pt. We note that the zAFM was indeed found for the bulk NiPS3 BREC19863 . In the case of the $4d$ $M$=Pd, the lowest-energy state is the AFM, while the obtained solutions are almost degenerate in energy among PM, AFM, and zAFM. The results indicate the trend from strongly correlated $3d$ to weakly correlated $5d$; the $4d$ case is in the boarder between the magnetic insulator and the Dirac paramagnetic semimetal. Thus, PtPS3 would be a prime candidate for the multiple Dirac semimetal rather than PdPS3. In addition, from the comparison of the magnetic moments between the GGA calculation for PdPS3 and the mean-field analysis at half filling [see Fig. 4(a)], the corresponding $U$ value in the mean-field analysis is about 0.5 eV, which appears to be weak for 4$d$ transition compounds. This may be ascribed to the generic tendency that the GGA calculation underestimates the electron interactions of $d$ orbitals.