Large Orbital Magnetic Moment in VI3

The existence of the V3+-ion orbital moment is an open issue of the nature of magnetism in the van der Waals ferromagnet VI3. The huge magnetocrystalline anisotropy in conjunction with the significantly reduced ordered magnetic moment compared to the spin-only value provides strong but indirect evidence of a large V orbital moment. We used the unique capability of X-ray magnetic circular dichroism to determine the orbital component of the total magnetic moment and provide a direct proof of an exceptionally sizable orbital moment of the V3+ ion in VI3. Our ligand field multiplet simulations of the XMCD spectra in synergy with the results of DFT calculations agree with the existence of two V sites with different orbital occupations and OM magnitudes in the ground state.


MAGNETIZATION MEASUREMENTS
We performed detailed magnetization measurements of VI3 bulk samples at 2 K to check the total value of V moment. The magnetization data was measured in magnetic field up to 7 T using a SQUID magnetometer MPMS7T (Quantum Design Inc.). To probe the angular dependence of magnetization in ac R -plane (at 2 K in 5 T) we used a homemade rotator with rotation axis orthogonal to the applied magnetic field. The angular-depended magnetization was measured in relative units and then scaled to the data measured along the c R -axis and ab-plane.
The results are shown in Figure S1. The maximum magnetization value of 1.23 B was found along a direction tilted by ~ 40º from the c R -axis.
Figure S1 -top: magnetization isotherms measured at 2 K. Bottom: angular dependent magnetization. The angular dependence was scaled according to the absolute values of magnetization measured along the c R -axis and ab-plane (see blue and red dots). The small discontinuities at around 90º and 270º are of an instrumental origin, most likely caused by the change of magnetic torque.

XMCD sum rules and Oxygen contribution
We used magneto-optical sum rules 2,3 to analyze the direction and value of the projection of an orbital magnetic moment (morb). The spin sum rule requires separating the XMCD contribution from the L3 and L2 edges. Correction factors of the spin sum rule have been calculated for the late transition metals 4,5 . However, for Ti, V, and Cr the overlap is large and no correction factor can be properly calculated. Therefore, we focus only on the orbital sum rules given by the expression below: where ̂ is the x-ray propagation unit vector, q and r are defined in figure S2, ℎ stands for the number of holes per atom in the 3d orbital, is Bohr magneton. In expression (1) the quantity r is calculated from the XAS defined as the sum of individual spectra measured with left and right circularly polarized X-rays. Figure S2top: XAS and bottom: XMCD measured at 90 K normal incidence, 5 T applied field parallel to the X-ray direction. Figure S3 shows the XMCD and corresponding integral measured at 2 K, 5 T, normal incidence. The integral value indicated is calculated as the average in the energy region 540-580eV and the standard deviation error bar. This data was used for the OM calculation presented in the manuscript.

Figure S3
-XMCD measured at 2 K, 5 T normal incidence and corresponding integral. The integral value indicated is the average within the energy range 540-580 eV and the error bar is the standard deviation.
As discussed in the manuscript, the XAS measured at low temperatures presents an additional contribution of the oxygen K-edge. The base pressure of the chamber during the measurements is below 5e-11mbar, so a residual gas containing lots of oxygen is unlikely to be the origin. One of the main differences between the 90 K or 300 K compared to the 2 K measurement is that the TEY drain current is strongly reduced at 2 K, since VI3 is a semiconductor. The TEY at the V-L3 absorption is reduced by a factor of almost 40 between 90 K and 2 K. This means that if there is a 2% illumination of the sample plate (which contains oxides) that will not be seen at 90 K since the Vanadium signal is much stronger, but, at 2 K, the Vanadium TEY has approximately the same magnitude as this small contribution from the sample holder. The xray beam used is out-of-focus to minimize radiation damage. In this condition, the 4 vertical footprint of the x-ray beam, which should contain about 95% of the x-ray intensity is roughly matching the crystal vertical size. Therefore a few % of x-rays illuminating the sample holder is quite plausible.
The O K-edge presence makes a direct calculation of the XAS integral difficult. There are different ways of estimating the XAS integral for the sum rule, which we discuss in sequence. Figure S4 shows the XAS spectra and corresponding integral for the 2 K and 300 K data. One way to calculate the XAS integral is to adjust the step function to the 2 K XAS spectrum around 527 eV, subtracting the step function from the XAS and taking the integral at this energy. The  A second method would be to take the integral of the data measured at 300 K, where the oxygen spectrum does not interfere, and the step function can be normalized at high energy. In this case, the step function is normalized around 570eV and the corresponding integral after the step function removal is 4.81. The resulting OM is 0.53 B. This method is better than the previous one because we can calculate the full XAS integral more precisely. The disadvantage is that the XAS at 300 K is not identical to the one at 2 K. As seen in figure S4, the relative ratio of L3 and L2 peaks is different.
For this reason, we suggest a third method, which is to estimate the difference in the XAS integral, if the step function is normalized at 527 eV or at 570 eV and use that correction factor to estimate the correct XAS integral for the 2 K spectrum if the step would be normalized at 570 eV. The correction factor was calculated using the 300 K data. As shown in figure S5 the XAS integral of the 300 K data taken at 527 eV, with the step function normalized to the data at this energy, is 3.53. Therefore, the integral taken just above the L2 edge corresponds to 3.53/4.81*100=73% of the integral taken at 570 eV. Using this relative difference between the 300 K integrals at two different energy points, we can attempt to correct the 2 K XAS integral, which can only be taken in a restricted energy range. We obtain: 3.76/0.73 = 5.12 as an estimate for the 2 K XAS integral value if there would be no Oxygen contribution to the spectrum. If we use this value for the XAS integral, we obtain an OM of 0.49 B.
In summary, we present three different methods to calculate the XAS integral and the OM.
Each method has its advantages and disadvantages. The main message of these calculations is to show that the OM for vanadium in VI3 is significantly large no matter which method we use.
As a representative estimate of OM, we took the average value of the results of the three aforementioned methods and used the standard deviation among them as the error bar. This comes down to OM of 0.6(1) B.

MULTIPLET SIMULATIONS
The X-ray absorption spectra and x-ray magnetic circular dichroism spectra were simulated for V 3+ d 2 using ligand field multiplet simulations based on atomic Hartree-Fock calculations using Cowan's code 6 followed by chain symmetry proposed by Butler 7 . This approach has been used for the first time by Thole and van der Laan for the description of x-ray absorption spectra 8 .
The local V symmetry used in the simulations was D3d which is reduced to C3 when the magnetic moment is added. Slater Integrals were reduced by 50% from Hartree-Fock values.
All parameters used are summarized in Table S1 -parameter used in multiplet simulations shown in Figure 4 of main text.
The energy levels in trigonal symmetry are related to the CFS parameters as shown in equations S(3-5).

XAS, XMCD MEASUREMENTS
X-ray magnetic circular dichroism measurements were carried out using the EPFL-PSI endstation in the X-Treme beamline 10 , Swiss Light Source. The crystals were mounted inside a glove box and cleaved for a fresh surface. Then the crystals were transferred in an inert atmosphere from the glove box to the load lock, where they were cleaved for a second time and then pumped down. To improve conductivity for the total electron yield measurements, the crystals were capped with a few nm carbon layer evaporated in situ in the ultra-high vacuum sample preparation chamber. The measurements shown were done with the x-rays perpendicular to the exfoliation plane (parallel to the c-axis in rhombohedral structure, c R ). Figure S6 shows two measurements done during the same beamtime with a time difference of 6.5 hours. These measurements show that no sign of change in the XAS nor XMCD is observed due to the long incidence time of the x-rays.

DFT CALCULATIONS
Density functional theory (DFT) calculations employed the full-potential linear augmented plane wave (FP-LAPW) method, as implemented in the band structure program ELK 11 . The generalized gradient approximation (GGA) parametrized by Perdew-Burke-Ernzerhof 12 has been used to determine the exchange-correlation potential. SOC is known to play a key role in the formation of OM and has been included in the calculation. DFT simulations utilizing GGA-PBE have already successfully described quasi-2D compound VI3 13 . The full Brillouin zone has been sampled by 10 × 10 × 5 k-points and the convergence w.r.t. k-mesh density has been verified. An increased accuracy of expansion into spherical harmonics has been used with lmax=14. The simulation was performed for a bulk system consisting of weakly coupled layers.
Since the material is known to be a The states | + , ⟩, | − , ⟩ originate from the following linear combinations of spherical harmonics, the former 2 states of the cubic system 17 : These are eigenstates of Lz operator with values ±1, and therefore remain degenerate w.r.t.
trigonal distortion, but are split by SOI.

NOTE S1
For temperature-induced effects energy 5.6 meV (~ 65 K) is high enough to be not excited at considered temperatures (2 K). However, even a small distortion of the lattice introduces significantly larger changes in total energy. Therefore, in the presence of various point defects or stacking faults, the state previously located 5.6 meV higher could become favorable. Another aspect is that in this situation we have a local minimum separated by 5.6 meV from the global minimum. These appear to be relatively far from each other in phase space, and there is also the metallic solution as another local minimum. This increases the risk that self-consistent calculation finds local minimum that is not a true global minimum, which could explain why different results were got by different groups.