Strongly correlated two-dimensional (2D) magnets are a puzzling class of materials from a fundamental physics perspective. With the discovery of ferromagnetism (FM) in 2D systems such as CrI 3 1 and Cr 2 Ge 2 Te 6 , 2 an effort to identify and understand the underlying mechanisms of 2D magnets with a finite transition temperature (T c ) has become a highly active area of materials science. One of the most interesting and controversial 2D ferromagnetic materials is VSe 2 , which has a metallic 1T phase (octahedral (1T)-centered honeycombs) and a semiconducting 2H phase (trigonal prismatic (2H)-hexagonal honeycomb). 3 Despite discrepancies of the structural properties coupled to the energetic stability (whether or not 1T vs 2H is more favorable) [4][5][6][7] that we were successfully able to resolve in our previous work using highly accurate electronic structure methods (see ref 8), there are several remaining questions regarding 2D 1T-VSe 2 .
In a 2018 study, single layer 1T-VSe 2 was synthesized on van der Waals (vdW) substrates (graphite and MoS 2 ), where strong ferromagnetic ordering was measured above room temperature and a charge density wave (CDW) was detected with a transition temperature of 121 K. 9 This was the first study to demonstrate that the CDW transition of single layer 1T-VSe 2 is coupled with its magnetic properties. 9 Despite room temperature ferromagnetism being reproduced in some cases, such as for a chemically exfoliated monolayer (ML), 10 there have been several conflicting reports of the magnetic properties and the interplay between the CDW state and magnetism. [9][10][11][12][13][14][15][16][17][18][19][20][21][22] For example, some studies have found the nonmagnetic (NM) CDW phase to be experimentally favorable, with an inherent absence of ferromagnetism. [13][14][15][16] Several theoretical studies have attempted to explain why these discrepancies in magnetic properties might occur, citing strain, vacancies, substrate choice, doping, and chemical functionalization as possible reasons. [23][24][25][26][27][28][29][30][31] These studies also highlight the strong sensitivity of the magnetic properties of 2D 1T-VSe 2 to extrinsic factors, which can provide a viable route to tune such properties. In our previous density functional theory (DFT) work, we explored the competing magnetic and nonmagnetic states in single layer 1T-VSe 2 with and without charge density wave. 26 We found that there is strong competition between nonmagnetic and magnetic states in the CDW structures (with respect to the undistorted structure), with relative energies being on the scale of 1 meV per formula unit (fu). 26 This implies that it is possible for antiferromagnetic (AFM) ordering to compete with ferromagnetic ordering and the CDW state. 26 Although our previous semilocal DFT calculations 26 provide a good qualitative assessment of how different magnetic orderings can compete in the undistorted and CDW structures, a more quantitative answer is required to accurately understand the magnetic and CDW transitions in monolayer 1T-VSe 2 and estimate quantities such as transition temperatures. For this reason, high-fidelity many-body techniques such as fixed-node Diffusion Monte Carlo (DMC) 32 can be utilized to accurately describe the electronic and magnetic properties of 2D 1T-VSe 2 . DMC is a correlated electronic structure method that has a reduced sensitivity to approximations such as the exchange-correlation functional and the Hubbard U 33 correction. In addition, DMC has successfully been applied to several 2D and quasi-2D systems. 8,[34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] In this work, we investigate the magnetic properties of monolayer 1T-VSe 2 in various magnetic states and geometries, with and without CDW distortions, through the lens of DMC and provide a thorough benchmark using several DFT functionals. In addition to providing accurate quantitative estimates of energy differences between different phases, we coupled our DMC results to classical Monte Carlo simulations to estimate magnetic transition temperatures. Finally, we used our DMC insights to perform DFT calculations of the Raman modes and compared to our own experimental results on a synthesized 1T-VSe 2 flake.
In this work, we specifically focused on monolayer 1T-VSe 2 in its freestanding form. Motivated by the experimental discrepancies of controversial room temperature ferromagnetism, competing magnetism and CDW states in the 1T phase, we decided to perform in-depth benchmarking at the DFT level using a variety of approximations. A summary of the different structures and magnetic configurations we studied is given in Figure 1. We focused on the normal (undistorted) structures of 2D 1T-VSe 2 in its FM, AFM (spins are antialigned in a stripy pattern), and nonmagnetic (NM) orientations. In addition to the normal unperturbed crystal
structure of 1T-VSe 2 , we studied the distorted √3 × √7 × 1 supercell, which is a signature of the material being in its CDW state. This distorted √3 × √7 × 1 CDW structure has previously been studied extensively with DFT and experimentally verified. Building on the work of ref 26, we performed a thorough DFT benchmark (using several different approaches) of the various magnetic states of the distorted CDW structure. These configurations included FM, NM, and various AFM configurations. Consistent with ref 26, we studied the √3 × √7 × 1 distorted AFM-A, AFM-B1 and AFM-B2 configurations (see Figure 1). It is important to note that the distorted AFM-A configuration is equivalent to the AFM (stripy) configuration for the undistorted structure.
In our previous work (ref 8), we determined the optimal geometry (including lattice parameter and bond distance) and relative phase stability of 2D 1T-and 2H-VSe 2 in the FM orientation with DMC and benchmarked with several DFT functionals with and without the Hubbard U correction. Figure S1a depicts the DFT benchmarking calculations for the NM, FM, and AFM orientations of the undistorted structure. We performed DFT calculations with LDA, PBE, SCAN, and r 2 SCAN with U values ranging from (0 to 3) eV, where we fully relaxed each structure in its respective magnetic orientation. As seen in Figure S1a, which plots the energy differences with respect to the NM state for each functional, the discrepancies with each DFT method are enormous. Although qualitatively the results are somewhat similar (i.e., the FM state being lowest in energy across the board for most DFT functionals), quantitatively, the results vary drastically. In fact, the energy per formula unit (fu) between the NM orientation and the FM/AFM orientations can vary up to ≈0.4 eV. A quantitative energy difference between magnetic states is absolutely necessary to obtain magnetic exchange parameters (J) and therefore transition temperature (T c ) with reasonable accuracy.
In order to overcome these shortcomings of local and semilocal DFT, we performed DMC calculations to obtain accurate total energies of the NM, FM and AFM orientations of the normal undistorted structure of monolayer 1T-VSe 2 . As a starting geometry for our DMC calculations, we used the structure obtained from our previous work in ref 8. We acknowledge that the local magnetic ordering can impact the structural properties, but for our DMC calculations we used the same structure for NM, FM, and AFM calculations. In order to study the influence of using the same geometry for NM, AFM and FM on the energetics, we performed DFT calculations with the fixed FM geometry for all magnetic orientations and found a maximum energy difference of 10 meV between the fixed and relaxed structures, which is within the error bar of our DMC calculations. To reduce finite-size errors, we performed DMC calculations for two reasonably sized supercells for NM, FM, and AFM and extrapolated to the thermodynamic limit. We used the following supercell sizes: 36 and 72 atoms for NM and AFM, and 27 and 48 atoms for FM. The DMC energy differences and uncertainties are depicted in Figure S1 and Table 1 (in shorthand notation). From the figure and table, we see that SCAN and r 2 SCAN (with no Hubbard U correction) are in the closest agreement with our DMC benchmark. For the FM configuration, we also find that SCAN and r 2 SCAN successfully reproduce the DMC geometry that was previously obtained in ref 8.
Although other semilocal DFT functionals such as PBE + U (U = 1, 2 eV) can do a sufficient job at reproducing our DMC geometry, they fail in terms of correctly capturing the quantitative energy differences of various magnetic states. Meta-GGA functionals (SCAN, r 2 SCAN) can be quite successful for 2D magnetic materials such as VSe 2 . SCAN and r 2 SCAN include the kinetic energy density (in addition to the electron density and its gradient) in the exchangecorrelation functional, which allows them to capture complicated electronic interactions. These functionals can reduce the self-interaction errors of typical GGAs and improve the description of correlation effects. SCAN satisfies all 17 exact known constraints that meta-GGA can satisfy. r 2 SCAN is a regularized version of SCAN designed to improve numerical stability without hindering accuracy (r 2 SCAN is meant to obtain similar results to SCAN). 50 This does so by breaking some of the exact constraints of SCAN.
Although these functionals can simultaneously capture multiple properties such as energy differences between magnetic phases and optimal geometry, SCAN and r 2 SCAN have been widely reported to overestimate on-site magnetic moments. [51][52][53][54] In the case of 1T-VSe 2 , we found the DMC onsite magnetic moment of the V atom to be 1.06(2) μ B in ref 8. We found the on-site moment to be 1.40 μ B with SCAN and 1.43 μ B with r 2 SCAN. Nonetheless, due to the success of multiproperty predictions with SCAN and r 2 SCAN with respect to our DMC benchmark calculations, we decided to use SCAN to obtain more complex quantities (that are difficult to obtain with DMC) such as magnetic anisotropy energies. For computationally demanding simulations such as for the Raman active modes, where it is not computationally tractable to use SCAN or r 2 SCAN, we used PBE + U (U = 1 eV) due to its relative qualitative agreement with DMC results (in terms of structural properties).
Figure S1c depicts the results for the CDW-distorted structures with various magnetic orientations. The initial supercells for these distorted structures were obtained from ref 26, where the V and Se atoms were initially displaced by 0.12 and 0.18 Å respectively (similarly to ref 15) prior to atomic relaxation. For our DFT benchmarking calculations, we adapted a slightly different approach for the distorted CDW structures (in contrast with the undistorted cells). For these supercells, we modified the lattice constant to be a direct multiple of the unit cell geometry we obtained in our previous work (ref 8) and fixed the cell dimensions. We then allowed the atomic positions (with initially displaced atoms) to relax with each DFT functional. Similarly to the undistorted structures in Figure S1a, there is a massive discrepancy in relative energy between each magnetic orientation (with respect to the NM undistorted cell) of the distorted CDW structures (up to 0.5 eV). We observe that the FM orientation is the lowest energy magnetic state for all functionals depicted in Figure S1b. In contrast to the undistorted structure in Figure S1a, the energy differences between the FM and various AFM orientations is much smaller. Our results are in good qualitative agreement with the DFT calculations performed in ref 26, but are shifted slightly in energy which can be attributed to fixing the lattice constant to the DMC obtained value. Similar to the undistorted case, we also find excellent qualitative agreement between SCAN/r 2 SCAN and PBE + U (U = 1 eV) for the energy differences between the FM and AFM configurations of distorted structures. To provide more insight on these DFT calculations, we performed DMC simulations for just the FM CDW-distorted structure at supercells of 30 and 60 atoms and extrapolated to the thermodynamic limit. For the starting geometry of this structure, we optimized the initially displaced atomic coordinates with r 2 SCAN and used a direct multiple of the unit cell geometry we obtained in our previous work (ref 8). This DMC energy difference and uncertainty are shown (with respect to the undistorted NM structure) in Figure S1c and Table 1. We find that the DMC energy difference lies between those of PBE + U (U = 1 eV) and PBE + U (U = 2 eV). Additionally, we find that the SCAN/r 2 SCAN and LDA + U (U = 3 eV) energy differences are close to the DMC value. More detailed benchmarking data can be found in Tables S1-S5.
From our DMC insight, we went on to estimate the bilinear isotropic exchange (J), anisotropic exchange (λ), and easy axis single ion anisotropy (A). The incorporation of spin-orbit within DMC is a relatively new development, with successful applications in band gap calculations 41 and identifying band inversions in topological materials. 55 For 1T-VSe 2 , we were mainly focused on spin-orbit-induced magnetic anisotropy. From our DFT calculations, we know that the magnetic anisotropy energy differences are on the order of fractions of an meV/fu, while our DMC uncertainties for total energy are on the order of 10 meV/fu. Since these spin-orbit-induced magnetic anisotropy energy differences are within our DMC uncertainty, we performed collinear DMC calculations (where spins can orient either up or down) for the undistorted FM and AFM orientations (depicted in Figure S1 and Table 1).
For the undistorted structure, we obtain a J value of 71 (10) meV (no spin-orbit contribution). Taking the insight from our DMC and DFT calculations, where we observed that SCAN quantitatively reproduced our DMC benchmarking results with closer agreement than r 2 SCAN (see Table 1), we went on to perform spin-orbit SCAN calculations to explicitly obtain anisotropy parameters λ and A (see Supporting Information (SI) for more detailed information). From these calculations, we obtained a λ value of -0.57 meV and a A value of -0.77 meV. With regards to the distorted CDW structure, we performed a DMC benchmark calculation for the FM orientation (see Figure S1 and Table 1). Due to the reasonable agreement with DMC and SCAN and for consistency with the undistorted case, we went on to calculate exchange and anisotropy parameters for the distorted structure using spinorbit SCAN calculations. In the case of SCAN (since the AFM-A CDW-distorted structure relaxed to the undistorted AFM orientation), the AFM-B2 configuration is the ground state AFM configuration. Using our CDW-FM and CDW-AFM-B2 noncollinear SCAN calculations, we obtained a J value of 21 meV, a λ value of -0.77 meV, and a A value of 0.27 meV (see SI for more details). We went on to use these accurately computed magnetic exchange and anisotropy parameters for our spin Hamiltonian model shown in eq 1, where we performed classical Monte Carlo simulations to gain insight into the magnetic phase transitions of the undistorted and distorted structures.
Figure 2 illustrates the magnetic phase transitions for both (a) undistorted and (b) CDW phases. In each case, the exchange interaction J is predominant over other interactions. Notably, J is over three times higher in the undistorted phase, measuring 71 meV, compared to 21 meV in the CDW phase. This pronounced disparity in interaction strength indicates that spins in the undistorted structure align more strongly with neighbors, leading to an extended correlation length and larger clusters of aligned spins, requiring more thermal energy to disrupt these interactions and disorder the aligned spins. This implies that there exists two potential magnetic phase transitions for freestanding (no substrate, strain, or defects) 1T-VSe 2 , one for the undistorted monolayer below room temperature (≈100 K less than the transition temperature reported in ref 9) and one for the CDW-distorted (slightly below liquid N 2 ). These results indicate that the magnetic phase transition occurs below the CDW transition temperature, which has been reported to be above 121 K. This is expected since the CDW phase must exist before the CDW-FM phase can exist by further cooling. The small energy scale between the CDW-distorted phases and the undistorted phases (reported in our current work and previous work 26 ) poses a challenge in understanding which phase is energetically more favorable. This small energy scale opens the door to tune the energetic stability of various magnetic states (with and without CDW-distortion) through external mechanisms such as strain, defects, substrate engineering, or temperature.
For a deeper analysis of how structural parameters are related to the energetic favorability of different magnetic states, we performed biaxial strain calculations for the undistorted and distorted structures with the most reliable DFT functional (in our case, SCAN was the most reliable for reproducing DMC geometries and quantitative DMC energy differences). These results are depicted in Figure 3. It is important to note that for the CDW-distorted structures, we kept the displaced atomic positions fixed while applying biaxial strain and excluded the AFM-A orientation since it relaxed to the undistorted AFM (stripy) orientation with SCAN. We also kept the relative atomic positions fixed for the undistorted configurations while varying the lattice.
From Figure 3a, we observe that the CDW-distorted structure in FM and AFM configurations is more sensitive to small amounts of strain, in contrast to the undistorted structure. In other words, the energy differences between the strain curves are much smaller for the CDW-distorted structures (compared to undistorted). Notably, we observe that for positive values of strain (lattice expansion), the FM in the undistorted and CDW-distorted structures becomes more stable. For the undistorted structure, we see that the energy difference between the FM state and AFM state drastically increases with the application of positive (tensile) strain. We observe that with small amounts of compressive strain, the undistorted FM structure becomes more favorable than the CDW-distorted FM structure, which can be a promising route to stabilizing the FM in the undistorted structure. For larger amounts of compressive strain (≈-2%), AFM order begins to become more favorable for the undistorted and CDWdistorted (specifically AFM-B1) structures.
To further understand this phenomena, we computed J under strain for the energy differences depicted in Figure 3a for the undistorted and CDW-distorted structures. We then went on to estimate T c under strain using these values of J for additional Monte Carlo simulations (similar to those reported in Figure 2). It is important to note that for these Monte Carlo simulations under lattice strain we did not include A and λ in the Hamiltonian (due to increased computational demand of spin-orbit SCAN calculations). Since the driving force behind the T c is J, we do not expect major fluctuations in T c from A and λ with applied strain. These results are depicted in Figure 3b. In addition to the strong competition between the undistorted FM and the CDW-FM structures (dotted purple and orange curves in Figure 3a), there is potential straintunability within each FM ground state (with and without CDW). In Figure 3b, we see that J can be enhanced in the CDW-FM state by small amounts of biaxial expansion or compression, with the T c approaching values above 120 K. With regards to the undistorted FM structure, we find that tensile strain can enhance the T c to be above room temperature, but there is strong competition with the FM-CDW state for small values of tensile strain, with the CDW-FM becoming significantly more energetically favorable after ≈+2% strain. In addition to the tunability of magnetic states, it was revealed from anharmonic phonon calculations that the CDW order can be tuned with small amounts of lattice strain (as little as 1.5%), highlighting the strong substrate dependence. 56 For example, common substrates for VSe 2 9 such as MoS 2 (a = b = 3.16 Å) and highly oriented pyrolytic graphite (a = b = 2.46 Å) have significant lattice mismatch, which can induce strain (since the lattice constant of monolayer 1T-VSe 2 is estimated to be ≈3.41 Å 8 ). As an alternative to substrate Figure 3. An illustration of the impact of strain on undistorted and CDW-distorted 1T-VSe 2 : (a) The total energy offset (E total -E 0 ) in eV/fu (calculated with the SCAN functional) as a function of biaxial strain percentage (where E 0 , the energy of the unstrained, undistorted FM configuration, is used as the reference point). The legend specifies the magnetic orientation of each set of data points, where the solid lines represent the distorted CDW structures and the dashed lines represent the undistorted structures. The inset depicts a zoomed in section of the curves, focusing on the lowest-energy orientations. (b) Isotropic exchange (J) for the undistorted (blue) and CDW-distorted (green) structures plotted on two separate axes. These values of J were determined by the undistorted (FM, AFM) and the CDW-distorted (FM,AFMB2) structures. The dotted regions of these two curves indicate a transition to AFM favorability. The DMC expected values and uncertainties for both J and the equilibrium lattice constant (for the undistorted FM structure) are also depicted. Finally, the star points on each respective curve represent the T c value computed under strain (with classical Monte Carlo) for that particular value of J.
engineering, substitutional doping can be a promising route to inducing tensile strain in 1T-VSe 2 , which can stabilize and even enhance the magnetic properties in both the undistorted and CDW-distorted phases.
For additional benchmarking purposes, we performed density functional perturbation theory (DFPT) calculations to obtain the assignment and peak positions of the Raman modes of 1T-VSe 2 for the distorted and undistorted structures in various magnetic and nonmagnetic configurations at 300 K. Due to the higher computational cost and convergence issues of SCAN, we performed these DFPT calculations with PBE + U (U = 1), since this functional was able to match our DMC benchmark for lattice geometry. The logic behind using PBE + U (U = 1) stems from the fact that if it can correctly reproduce our lattice geometry, it can accurately capture vibrational properties. Table 2 depicts the Raman active modes calculated with DFT for the bulk and monolayer structures (with and without CDW distortions). For the undistorted structure, we performed calculations for the FM, AFM, and NM orientations, and for the CDW-distorted structure, we considered the FM and orientations. As seen in Table 2, each structure has a pronounced A 1g peak ranging from (197 to 208) cm -1 . Due to a lack of symmetry in the CDW-distorted structure, we could only identify this A 1g , but it is entirely possible that other Raman modes can appear experimentally. For the undistorted structure (bulk and monolayer), we observe an E g peak ranging from (130 to 137) cm -1 for the FM and NM orientations. The most distinguishing features are apparent for the AFM orientation of the monolayer, where we observe an additional A 1g peak at 139 cm -1 and a B g peak at 133 cm -1 .
To validate our DMC-informed DFT calculations of the Raman modes, we grew a VSe 2 crystal by Chemical Vapor Transport (CVT) (see SI for details). The 1T phase of VSe 2 was confirmed by X-ray diffraction (XRD) (see Figure S2). Thin flakes of 1T-VSe 2 were exfoliated and encapsulated in h-BN (Figure S3a) and subsequent Raman spectra measurements were performed. The Raman spectra were averaged (30 s, 20 accumulations) at four distinct positions on a 1T-VSe 2 flake capped with h-BN. Figure 4a shows an optical image with the approximate locations where the spectra were measured, and Figure 4b shows the corresponding spectra. The thinnest region measured to be about 1.7 nm by an atomic force microscope, Figure S3, position 3, was used for comparison to DFT predictions. To assign the phonon modes, polarized Raman spectra were collected from the 1.7 nm region (position 3). Data is shown in both the parallel and cross configuration, i.e., the incoming laser excitation and the Raman scatter were either parallel to one another, or at 90°(cross), confirming the respective A and E symmetry of the modes of the VSe 2 . In the 1.7 nm thick region, we identify two modes; the A 1g at 206.7 cm -1 ± 0.4 and the E g at 137.5 cm -1 ± 0.3 in good agreement with the two modes predicted by DFT. The other features correspond to the Si/SiO 2 substrate, as can be seen when compared to measurement of the bare substrate (position 4). However, in the bulk-like region on the h-BN capped flake, position 1, there are additional observable features to the right of the A 1g mode that differ from the edgelike feature observed in position 3 coming from the substrate. Previously, a broad mode at about 257 cm -1 has been attributed as an E g mode of the 1T phase in a bulk sample; 57 however, we observe in some cases two distinct modes at 235.3 and 253.2 cm -1 in the bulk-like region in the h-BN capped flake and at 229.8 and 250.3 cm -1 in the bulk crystal, shown in Figure S4. Since these are only observed in the bulk-like capped region/crystal, and are not predicted by the DFT results, these modes are likely due Se. [58][59][60] We conclude that the 1T phase of VSe 2 only has two phonon modes which we identify as the higher intensity A 1g and the lower intensity E g .
We have provided a comprehensive investigation of the magnetic properties and strain response of monolayer 1T-VSe 2 using DFT and highly accurate DMC methods. By extensively benchmarking various DFT functionals against our DMC results, we demonstrated the significant impact of the exchange-correlation functional on the accuracy of magnetic properties and energetic stability of various magnetic states of 2D 1T-VSe 2 , with and without CDW. Our high-fidelity results aim to resolve previous discrepancies that exist in the theoretical and experimental literature for 1T-VSe 2 . Through classical Monte Carlo simulations informed by DMC and DFT, we estimated the magnetic transition temperatures for the undistorted (228 K) and CDW-distorted (68 K) phases, and our strain calculations indicated that small amounts of biaxial strain can enhance the transition temperature, which provides a viable route for engineering the magnetic properties. Furthermore, our Raman spectroscopy experiments on exfoliated flakes of 1T-VSe 2 validate our predictions. Our work underscores the important role of highly accurate manybody methods such as DMC in describing the electronic and magnetic properties of 2D materials. We hope that our combined approach can pave the way for future explorations of correlated 2D magnetic materials.
Benchmarking DFT calculations for monolayer 1T-VSe 2 were performed with the Vienna Ab initio Simulation Package (VASP) code, using projector augmented wave (PAW) pseudopotentials. 61,62 For testing purposes, we employed a variety of exchange-correlation functionals including the local density approximation (LDA), 63,64 Perdew-Burke-Ernzerhof (PBE), 65 the strongly constrained and appropriately normed (SCAN) 66 meta-GGA functional and the r 2 SCAN 50 meta-GGA functional, with and without the Hubbard correction U, 33 which was used to treat the on-site Coulomb interaction of 3d orbitals of V. A k-grid of 24 × 24 × 1 was used and scaled proportionally to the supercell size, a plane-wave kinetic energy of 400 eV (increased to 800 eV for SCAN and r 2 SCAN) was used, and at least 20 Å space of vacuum was added between periodic layers. Our density functional perturbation theory (DFPT) 67 phonon calculations to obtain the peak positions of Raman modes were also performed with VASP along with the phonopy 68,69 package. Experimental peak positions and uncertainties are given in the last column.
Our quantum Monte Carlo (QMC) calculations for 1T-VSe 2 used the exact same settings as our previous work on 1T-and 2H-VSe 2 (see ref 8 for more specific details regarding kinetic energy cutoff, k-point grid, finite-size, and time step convergence). DFT-PBE was used to generate the trial wave function (nodal surface) for DMC using the Quantum Espresso (QE) 70 code with a U correction of 2 eV. This was due to the fact that PBE + U = 2 eV resulted in a trial wave function with the lowest total energy calculated with DMC. 8 We used a k-point grid of 12 × 12 × 1 and a kinetic energy cutoff of 4080 eV (300 Ry) to generate the trial wave function. In terms of pseudopotentials, RRKJ (OPT) potentials were used for V 71 and Burkatzki-Filippi-Dolg (BFD) were used for Se. 72 Variational Monte Carlo (VMC) and DMC 32 simulations were conducted using the QMCPACK 73,74 code and the Nexus 75 workflow manager was used to automate the DFT-VMC-DMC calculations. Up to twobody Jastrow parameters [76][77][78] were optimized using the linear method 79,80 to minimize the variance and energy in VMC, with the goal of modeling electron correlation and reducing uncertainty in the DMC results. 80,81 In order to compute the nonlocal part of the pseudopotentials, the locality approximation 81 was used in DMC. An optimal time step of 0.01 Ha -1 (0.27 eV -1 ) 8 was used in DMC, and to reduce finite-size errors, canonical twist averaging (with uniform weighting) 82 and extrapolating to larger supercells (up to 72 atoms) was performed. A more detailed explanation of the theory behind VMC and DMC can be found in ref 32.
Energies of different spin configurations derived from firstprinciples calculations were mapped onto a classical spin Hamiltonian, characterizing three-dimensional spins arranged on a two-dimensional lattice at sites of the V atoms. Subsequently, this Hamiltonian served as the basis for executing classical Monte Carlo simulations aimed at estimating the transition temperature (T c ) for both the undistorted and CDW phases. The Hamiltonian, which has previously been used for other 2D magnetic systems, 37,83,84 includes the bilinear isotropic exchange (J), anisotropic exchange (λ), and easy axis single ion anisotropy (A),
The indices i and j iterate over all magnetic sites and their corresponding first nearest neighbor (NN) magnetic sites, respectively. Only the first NN interactions are considered due to the strong localized magnetic moment on V. We calculated J 2 /J 1 for the CDW phase to be 0.068, supporting these assumptions (using SCAN).
According to the adopted sign convention, J > 0 favors isotropic spin alignment, leading to a preferred ferromagnetic (FM) phase over an antiferromagnetic (AFM) phase. Similarly, λ > 0 favors alignment of the spin z-components, and A > 0 favors the out-of-plane direction as the easy axis. We employ a 40a × 23 √3a supercell, equivalent to a nearly 13.6 × 13.6 nm 2 square cell, for the Monte Carlo simulations. The spin configuration space is discretized into uniformly distributed points on the surface of a unit sphere with a fine resolution of 0.5°, facilitating detailed exploration of spin orientations. The Monte Carlo simulation consists of numerous sweeps across all of the magnetic sites. During each sweep, the spin orientation at each site is updated
randomly, one site at a time, and the new configuration is then either accepted or rejected based on the Metropolis algorithm. The convergence of both energy and magnetization is meticulously monitored across varying temperatures. Upon achieving equilibration, an ensemble of configurations is collected to calculate the magnetization and specific heat. We use an ensemble size of 20,000 for the CDW simulation and 60,000 for the undistorted simulation. Additional methodology details can be found in the SI.
The authors declare no competing financial interest.
The code used to perform classical Monte Carlo simulations using the spin Hamiltonian is available at https://github.com/ UMBC-STEAM-LAB/SpinMCPack. Please note that the use of commercial software (VASP) does not imply recommendation by the National Institute of Standards and Technology. Certain commercial equipment or materials are identified in this paper to adequately specify the experimental procedures. In no case does the identification imply recommendation or endorsement by NIST, nor does it imply that the materials or equipment identified are necessarily the best available for the purpose. Notice: This manuscript has been authored by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The US government retains and the publisher, by accepting the article for publication, acknowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US government purposes. DOE will provide public access ACS Nano to these results of federally sponsored research in accordance with the DOE Public Access Plan (https://www.energy.gov/ doe-public-access-plan).
https://doi.org/10.1021/acsnano.4c15914 ACS Nano 2025, 19, 9925-9935