  This study investigates the spin-orbital states of FePS3 and FePSe3
monolayers and the origin of their Ising zigzag AFM, using DFT, crystal field
level diagrams, superexchange analyses, and parallel tempering MC simulations.
Our calculations show that under the trigonal elongation of the FeS6 (FeSe6)
octahedra, the $e_g^\pi$ doublet of the Fe 3d crystal field levels lies lower
than the $a_{1g}$ singlet by about 108 meV (123 meV), which is much larger than
the strength of Fe 3d SOC. Then, the half-filled minority-spin $e_g^\pi$
doublet of the high-spin Fe$^{2+}$ ions ($d^{5\uparrow,1\downarrow}$) splits by
the SOC into the lower $L_{z+}$ and higher $L_{z-}$ states. The spin-orbital
ground state $d^{5\uparrow}$$L_{z+}^{1\downarrow}$ formally with $S_z$ = 2 and
$L_z$ = 1 gives the large z-axis spin/orbital moments of 3.51/0.76 $\mu_{B}$
(3.41/0.67 $\mu_{B}$) for FePS$_3$ (FePSe$_3$) monolayer, and both the moments
are reduced by the strong (stronger) Fe 3d hybridizations with S 3p (Se 4p)
states. As a result, FePS3 (FePSe3) monolayer has a huge perpendicular
single-ion anisotropy energy of 19.4 meV (14.9 meV), giving an Ising-type
magnetism. Moreover, via the maximally localized Wannier functions, we find
that the first nearest neighboring (1NN) Fe-Fe pair has large hopping
parameters in between some specific orbitals, and so does the 3NN Fe-Fe pair.
In contrast, the 2NN Fe-Fe pair has much smaller hopping parameters and the 4NN
Fe-Fe pair has negligibly small ones. Then, a combination of those hopping
parameters and the superexchange picture can readily explain the computed
strong 1NN ferromagnetic coupling and the strong 3NN antiferromagnetic one but
the relatively much smaller 2NN antiferromagnetic coupling. Furthermore, our
PTMC simulations give TN of 119 K for FePS3 monolayer and also predict for
FePSe3 monolayer the same magnetic structure with a close or even higher TN.
