Multi-orbital physics of edge-magnetism in a Hubbard model of transition-metal dichalcogenide nanoribbons: Comparing Mean Field Theory and Determinant Quantum Monte Carlo

We study the emergence of edge-magnetism in zigzag nanoribbons from electron-electron interactions in a minimal 3-band tight-binding description of transition-metal dichalcogenide monolayers with an intra-orbital Hubbard term. We explain the influence of each orbital in the magnetic ordering, comparing the results of Mean Field Theory and determinant Quantum Monte Carlo. We focus on a gapped edge-antiferromagnetic phase appearing at threequarter edge-filling for realistic values of the interaction.


Introduction
Transition-metal dichalcogenides (TMDs) are part of the landscape of 2D materials beyond graphene [1]. These materials have a number of possible technological applications [2][3][4]. They are advantageous mainly because, unlike graphene, they are semiconductors in the monolayer form, with an intrinsic gap that exceeds 1 eV [5,6]. Additionaly, TMDs are good candidates for use in spin and valleytronics, where it is crucial to control the electronic spin and the valley degree of freedom.
We focus on molybdenum disulfide, MoS 2 , a part of the family of group 6 TMDs, denoted MX 2 , with M = Mo, W and X = S, Se, Te. We consider an atomically thin (∼ 6.5 Å), trigonal prismatic (2H) structure with ABA stacking. The structure consists of a layer of Mo-atoms in a triangular lattice, sandwiched between two layers of S-atoms, also organized in triangular lattices, as shown in the a) and b) panels of figure 1. We assume that the layers are much longer on one direction than on the other. This type of geometry is known as a nanoribbon since the length along the narrow direction is typically of the order of the nanometer.
The typically strong electron-electron interactions in transition-metal atoms potentially lead to the emergence of magnetism in TMDs. In fact, a ferromagnetic spin-valley polarized phase was recently predicted in a minimal model of a hole-doped TMD monolayer [7]. In this work, we probe TMD nanoribbons for magnetic order emerging from intra-orbital Hubbard interactions. Using Mean Field Theory (MFT) and determinant Quantum Monte Carlo (detQMC), we predict and characterize an edge-antiferromagnetic phase, appearing in the minimal model used here when the edges are at three-quarter filling, that is when 3/4 of the electronic edge-states are occupied.

Model
The model we use is based on the 3-band tight-binding (3BTB) model introduced in [8] with an additional Hubbard term that penalizes double orbital occupation. The Hamiltonian reads where i, j are sites on the triangular lattice, σ is the electron spin, and α, β are orbitals from the basis set are the hopping integrals (which depend on the orbitals, spin, and which of the 6 neighbors the electron hops to, with on-site terms also being allowed) [8]. Here, R i, j are lattice vectors and φ α,β (R i, j ) are Wannier states corresponding to orbitals α, β localized at R i, j . U is the Hubbard interaction. The 3BTB model captures edge-physics well. It is based on the symmetries of the monolayer and the fact that at low energies, both band edges have main contributions from the d z 2 ,d xy , and d x 2 −y 2 orbitals of the M-atoms. Spin-orbit coupling (SOC) leads to a splitting of the energy bands, with ε k,σ = ε −k,−σ (as imposed by time-reversal symmetry and shown in the panel d) of figure  1). However, this energy scale is much smaller than the typical Hubbard interaction required for magnetic ordering. Thus, we neglect SOC in what follows. In MFT, the Hamiltonian of equation (1) is replaced by the Hamiltonian that best describes the physics of the system, while being quadratic in the operators c ( †) i,α,σ , H MF . We minimize the variational grand potential Ω V , an upper bound on the real grand potential Ω: where T is the temperature of the system, S MF and H MF are, respectively, the entropy and energy computed with the mean field Hamiltonian, i.e. using the distribution ρ MF = e −βH MF , and µ is the chemical potential. The total particle number operator N commutes with the Hamiltonian, [H, N] = 0. Applying the variational principle to equation (1) and assuming spin is still a good quantum number, we obtain

Methods
Finding a self-consistent solution n i,α,σ corresponds to finding a minimum of Ω V . The global minimum corresponds to the mean field solution. We wish to reach the mean field solution starting with unbiased random or paramagnetic initial mean fields. To do so, we successively diagonalize the mean field Hamiltonian and use the obtained spectrum and eigenfunctions to recompute the mean fields, repeating the process until the mean fields vary only negligibly from iteration to iteration [9]. We use the grand canonical ensemble, fixing µ and obtaining the edge filling, which we fix using a bissection method. Starting from 8 random and 8 paramagnetic uncorrelated initial conditions, we find self-consistent solutions and choose the one that reaches the lowest variational grand potential. To improve convergence, we use annealing, starting from a high temperature and decreasing it smoothly until the desired temperature is reached. The dimensionality of our minimization problem is reduced by using the longitudinal translation symmetry of the ribbon, while considering 2 sublattices to allow antiferromagnetic solutions (see panel c) of figure 1). The determinant Quantum Monte Carlo method is based on a mapping to a higherdimensional non-interacting system using an auxiliary field. The Hubbard interaction is mediated by this discrete, and binary field h. The mapping leads to a problem of independent fermions coupled to h, and the fermionic part of the partition function may be traced out explicitly [10,11]: where S is a fermion-field action, p({h}) is the statistical weight of a configuration of the field, I is the identity matrix, imaginary-time is discretized in L = β∆τ slices labeled by index l, and B l,σ ( h l ) are square matrices with the dimension of the size of the system that depend on the parameters of the model and on the field h at each imaginary-time slice l. We use the grand canonical ensemble: H K → H K − µN. Using Wick's theorem, observables may be written in terms of the Green's matrices for a fixed configuration of the field h, The spin-spin correlator is measured by sampling configurations of h using Monte Carlo and then averaging over the field to obtain an estimate of S z i,α S z j,β ≡ (n i,α,↑ − n i,α,↓ )(n j,β,↑ − n j,β,↓ ) . Here, we compare the results obtained with MFT and detQMC at three-quarter edge filling for a MoS 2 nanoribbon, focusing on the role played by each orbital in the emergence of edge-antiferromagnetic order. In the left panel of figure 2, we show part of the mean field band structure, and we show how the staggered magnetization depends on U at T = 0 (see right panel of figure 2). Although this is not shown, in MFT, we have first order phase transitions, occurring, for the most part, above or near room temperature. The transitions to antiferromagnetic phases occur for different critical interactions on the two edges: U X c ≈ 1.84 eV, U M c ≈ 2.76 eV. The X-edge is more susceptible to antiferromagnetic order since the transition occurs sooner: U X c < U M c . Different orbitals are responsible for the magnetism of the two nonequivalent edges. In panels a) (bottom) and b) of figure 3, we see that the imbalance in the number of electrons occupying the d x 2 −y 2 orbital is the main cause of magnetic order on the X-edge, and similarly for the d xy orbital for order on the M-edge. The d z 2 orbital contributes significantly for both edges, but more on the X-edge. Thus, the staggered magnetization is higher on the X-edge. This could justify the fact that U X c < U M c (see panel c) of figure 3 from an orbital perspective). From the point of view of the mean-field bands, it is obvious that the AF ordering of the Xedge implies the opening of a gap at the Fermi level and, consequently the minimization of the variational potential of the system. As can be seen in panel a) of figure 3 (bottom), the orbitals have equal occupations on the bulk, inhibiting magnetic order. However, at the edges, sublattices A and B have different occupations, leading to antiferromagnetic order. Our detQMC results confirm the predictions of MFT. In figure 4, we can see that while antiferromagnetic spin-spin correlations decay in detQMC, unlike in MFT (in which case they are constant), they decay slower for the d z 2 and d x 2 −y 2 orbitals at the X-edge, and for the d z 2 and d xy orbitals at the M-edge, decaying quickly in the bulk in both cases.

Conclusion
We studied the emergence of edge-antiferromagnetism from the Hubbard interaction in zigzag TMD nanoribbons. Using Mean Field Theory and determinant Quantum Monte Carlo, we focused on the multi-orbital aspects of the magnetic ordering, explaining how antiferromagnetism sets in on the transition-metal and chalcogen-terminated edges at three-quarter edge filling, comparing the results obtained with the two methods. We found a greater tendency towards magnetic ordering on the X-edge compared with the M-edge, with a lower critical interaction U X c < U M c required for order to set in at T = 0.