The role of electron correlations in the electronic structure of putative Chern magnet TbMn$_6$Sn$_6$

A member of the RMn$_6$Sn$_6$ rare-earth family materials, TbMn$_6$Sn$_6$, recently showed experimental signatures of the realization of a quantum-limit Chern magnet. In this work, we use quantum Monte Carlo (QMC) and density functional theory with Hubbard $U$ (DFT$+U$) calculations to examine the electronic structure of TbMn$_6$Sn$_6$. To do so, we optimize accurate, correlation-consistent pseudopotentials for Tb and Sn using coupled-cluster and configuration-interaction (CI) methods. We find that DFT$+U$ and single-reference QMC calculations suffer from the same overestimation of the magnetic moments as meta-GGA and hybrid density functional approximations. Our findings point to the need for improved orbitals/wavefunctions for this class of materials, such as natural orbitals from CI, or for the inclusion of multi-reference effects that capture the static correlations for an accurate prediction of magnetic properties. DFT$+U$ with Mn magnetic moments adjusted to experiment predict the Dirac crossing in bulk to be close to the Fermi level, within $\sim 120$ meV, in agreement with the experiments. Our non-stoichiometric slab calculations show that the Dirac crossing approaches even closer to the Fermi level, suggesting the possible realization of Chern magnetism in this limit.

Experimental studies show that stable compounds form for R = {Sc, Y, Gd, Tb, Dy, Ho, Er, Tm, Yb, Ly} [18][19][20][21][22][23][24][25]. In this study, R = Tb is of particular interest since it is the only element of the previous set that forms an out-of-plane spin system. In contrast, other R-element materials display either an in-plane or canted spin direction behavior [1,2,8,18,23,24]. In particular, this case is interesting because the Chern magnetic phase predicted by Haldane's model requires an out-of-plane spin alignment [26]. Initially, Haldane considered a honeycomb lattice [26]; however, the model was later extended to kagome systems [27]. Therefore, the kagome lattice with out-of-plane spin alignment requirement seems to be met by Mn atoms of TbMn 6 Sn 6 and therefore opens * annaberdiyea@ornl.gov † subhasish.mandal@mail.wvu.edu ‡ ganeshp@ornl.gov the possibility of realizing a Chern magnet in this material.
Indeed, a few years ago, Yin et al. [1] carried out experimental transport measurements using scanning tunneling microscopy (STM) in high-magnetic fields and observed a series of Landau quantization states in TbMn 6 Sn 6 . They conducted AHE measurements and observed a large intrinsic contribution due to the large Berry curvature of Chern-gapped Dirac fermions. Further edge-state conductivity measurements showed significant signals within the Chern gap energy scale. Therefore, they concluded that TbMn 6 Sn 6 is close to realizing a quantum-limit Chern magnet, as predicted by the Haldane model. Indeed, kagome lattice geometry with an out-of-plane magnetization formed by the Mn atoms and the presence of strong spin-orbit coupling (SOC) stemming from the Tb and Sn atoms to open a Chern gap seems to suggest that TbMn 6 Sn 6 could realize such an exotic phase.
However, the possibility of realizing the Chern phase hinges on the ability to align the Fermi level (E F ) in the Dirac crossing (DC), or the Chern gap, seen in the kagome systems [1,2,27]. Yin and coworkers' experiments show that the Chern gap is located at 130(4) meV above E F [1]. However, Lee and coworkers reported density functional theory with Hubbard U (DFT+U ) calculations [2] predicting the quasi-2D DC to be 700 meV above E F , or at least 300 meV above E F , suggesting that the Dirac crossing cannot be related to the observed AHE signal. Therefore, there seems to be a disagreement in the literature about the exact position of DC with respect to E F and whether TbMn 6 Sn 6 could realize a Chern phase. In this work, we focus on whether the DC is indeed close to the Fermi level using quantum Monte Carlo (QMC), DFT+U , and other correlated methods.
Clearly, accurately estimating the energetic level of DC with respect to E F , which is on the order of a few hundred meV, would require a very accurate electronic structure treatment. Such a small energetic level could be affected by the accuracy of pseudopotentials [28], the level of electron correlation treatment, the accuracy of SOC effects, and the degree of multi-reference character intrinsic to the true electronic wave function of the system. In the following sections, we thoroughly investigate these aspects and try to enumerate their contributions and biases.
In Section II A, we optimize the required Tb and Sn element pseudopotentials using correlated methods and evaluate their errors via transferability tests. Then in Section II C, we carry out nearly-exact energy calculations for related molecular systems to estimate the biases present in single-reference QMC calculations. We then apply the developed pseudopotentials to TbMn 6 Sn 6 using DFT+U in Section II D and QMC in Section II E and demonstrate the issues with single-reference calculations. Section II F includes additional DFT+U results, which elucidate the position of DC. The findings are further analyzed and discussed in Section III. in-plane views are shown. The large red spheres are Tb atoms, medium blue spheres are Mn atoms, and small gray spheres are Sn atoms.

A. Pseudopotentials
One of the major approximations used in correlated methods is the pseudopotential or effective core potential (ECP) approximation [29]. High-accuracy correlated methods such as configuration interaction (CI) or coupled cluster (CC) are limited by the accuracy of the ECP as it modifies the Hamiltonian. Therefore, it's crucial that ECPs be optimized in a correlation-consistent way, especially for correlated methods and systems. Here, we develop correlation-consistent ECPs (ccECP) for Sn and Tb elements using CC and CI methods and use a ccECP for Mn from a previously published work [30]. In the next subsection, we show the results for the Tb element and compare its accuracy against fully correlated, relativistic, all-electron (AE) calculations. The full results for Sn, optimization methods, and further details are provided in Supplementary Note 2.

B. Tb ccECP
The optimization of the Tb ECP is more involved than Sn (Supplementary Note 2) or Mn [30]. For instance, in Sn, the electrons can be easily partitioned by quantum principle number n into core and valence spaces by treating n ≥ 5 in the valence space (5s, 5p, 5d, ...). However, in Tb, core↔valence partitioning is not straightforward since 4f orbital eigenvalues are on the same order as n = 5 principal number orbitals (5s, 5p, 5d), suggesting it should be included in the valence space, see Supplementary Tables 7 and 8. To overcome this, one could include the whole n = 4 electrons in the valence and treat only n ≤ 3 electrons as the core; however, this results in a prohibitively high computational cost in many-body QMC calculations. Conversely, 4f electrons are sometimes included in the core, but this can result in poor results [28]. In this work, we decided to include the 4f electrons in the valence (with N core = 46 e − ) and treat them fully self-consistently along with 5s, 5p, 5d, and 6s orbitals. However, the many-body accuracy of such partitioning is not well-known and must therefore be carefully tested. At least in 3d transition metals, similar large-core ECPs with valence spaces of 3d k 4s 2 (where the lowest n orbital is not an s) are well-known to produce inaccurate results [31][32][33].
Our extensive many-body calculations show that the chosen core choice for Tb ccECP should result in small errors once properly optimized. In fact, the accuracy of Tb ccECP is on par with the uncorrelated core allelectron (UC) where N core = 46 electrons are not active. This is illustrated by the atomic gap errors in Figure 2(a) and TbH 3 , TbO molecular binding energy errors in Figure 2(b-c). We observed that all atomic gap errors are below ∼0.07 eV up to the third ionization potential level, see Supplementary Table 9. Higher-order ionizations result in significantly larger errors on the order of an eV, showing the limits of such an ECP. On the other hand, very high ionization potential errors are probably not relevant for condensed matter applications since rare-earth elements form either R 2+ or R 3+ compounds [34]. Nevertheless, Tb ccECP is much more accurate than the existing ECPs with the same core choices, such as SBKJC [35] and an Opium-generated DFT-PBE [36] ECP using the Troullier-Martins scheme [37]. In fact, the ccECP errors are below chemical accuracy for the low-lying states (Figure 2(a)) and most of the bond lengths in hydride and oxide molecules (Figure 2(b-c)). Tb and Sn ccECPs are available online in Ref [38] for broader use.
The presented low biases of Tb and Sn ccECPs here and of Mn ccECP elsewhere [30] provide the opportunity to focus on the other systematic biases, such as the accuracy of the wavefunction and correlation treatment. Namely, we can assume that any deviations from the experiments are not due to the technical limitations of the representative Hamiltonians; rather, it has to be due to more fundamental reasons, such as inaccurate trial wave function in QMC or improper density functional approximation (DFA) and Hubbard U in DFT+U . These are the subjects of Sections II C and II D, respectively.

C. QMC Biases of Elemental Sub-Species
It is beneficial to know the inherent systematic biases of single-reference QMC calculations for the system of interest by comparing the QMC energies to CI or CCSD(T) energies at the complete basis set limit (CBS). Although CI and CCSD(T) calculations have been recently applied with great success to periodic systems [39][40][41][42][43][44][45], they are still out of reach for large systems such as TbMn 6 Sn 6 [40,45]. Therefore, we focus on constituent atoms and relevant molecular systems to estimate the QMC biases. Once the ECPs are made accurate, and other systematic errors such as timestep bias and walker population biases are under control, the main remaining errors are fixed-node and localization bias (if a nonlocal ECP is used as here). Usually, the fixed-node bias is the dominant of these two [46][47][48]; however, both of them go to zero as the trial wave function approaches the exact form [47,49,50].
We enumerate the biases of fixed-node diffusion Monte Carlo (DMC) [49] η as a percentage of the correlation energies in Figure 3(a) using nearly-exact energies from CCSD(T) at CBS, and single-reference DMC calculations with PBE [51], PBE0 [52], and Hartree-Fock (HF) references. This plot reveals a few insights about TbMn 6 Sn 6 solid. First, the Sn atom has the smallest fixed-node bias, ∼ 2% of the correlation energy. This error is close to the 1 − 2% errors of the iso-valent Si atom, Si x H y molecules, and diamond-structure bulk Si observed before [53,54]. Second, Mn shows a sizable error, ∼ 9% of the correlation energy missing in single-reference calculations, which also agrees with previous work [55]. Finally, the largest errors occur for systems with Tb, where ∼ 12% of the correlation energy is missing. However, we found that these large biases tend to cancel out for Tb when considering energy differences, such as binding energies of TbH 3 and TbO molecules (see Supplementary Table 17).
In TbMn 6 Sn 6 , the leading Tb atomic density of states (DOS) is well below the Fermi level (E F ), and Mn atomic DOS is dominant near E F as shown in Figure 3(b). Given that the Sn atom DMC errors are small and Tb states are well below E F , we argue that the largest errors in QMC will stem from an improper characterization of Mn atoms in our trial wave functions. In other words, low-lying conduction states near E F might give rise to wave functions of multi-reference character due to sizable single-reference DMC error in the Mn atom. This is supported by calculated atomic moments from DFT+U , which show severe sensitivity when U is applied on Mn-3d orbitals (shown later) and insensitivity on Tb-4f , Tb-5d, and Sn-5p orbitals, see Supplementary Table 15. In addition, we observe no changes in the band structure when +U is applied on Tb-4f orbitals (Supplementary Figure  3). Figure 3(a) is that for a given magnetization, Tb systems with HF reference result in the lowest energies signifying the importance of localization in 4f orbitals. On the other hand, Mn and Sn seem to favor more delocalized orbitals when compared to HF for constrained magnetization. These aspects will be discussed later in detail by means of applying Hubbard U on Mn-3d orbitals.

D. DFA Sensitivity Problem
This work uses the DFT+U approach with an effective Hubbard repulsion (U eff = U − J) [56,57] with the fully localized limit (FLL) double counting scheme [57][58][59][60] as implemented in quantum espresso code [61][62][63]. We explored two types of DFT+U calculations. In one case, the total magnetization of the cell was constrained to agree with the experimentally observed value of ∼ M s = 7 µ B in the SCF procedure (indicated with C-DFT, such as C-LDA or C-PBE). In the other case, the total magnetization was unconstrained and determined during the SCF cycles (simply indicated with the DFA name, such as LDA or PBE). Extended details of the employed DFT+U methodology are given in Supplementary Note 3.
TbMn 6 Sn 6 ground state is a ferrimagnetic phase (FiM), where the Mn atoms align ferromagnetically, and the Tb spins point in the opposite direction [24]. Figure 4 shows the bulk band structures using various DFAs for the FiM phase. For a pure 2D kagome lattice, tight binding models show that the DC occurs at E F , k = K point while the flat bands are above the E F [1,27]. Considering the LDA [64] band structure ( Fig. 4(a)), we see the expected 'flat' band in the minority spin channel at ∼0.5 eV. In addition, a few DCs occur at k = K in the minority spin channel. These DCs were labeled for clarity in Figure 4 and will be referred to as DCn where n is the label. When LDA is switched to PBE DFA ( Fig. 4(b)), the flat band and DC2, DC3 shift up by around 200 meV. Introducing a small effective Hubbard U 3d Mn = 0.5 eV to LDA (Fig. 4(c)) shifts the flat bands and DCs up even further. The upward shifts in the minority spin channel are accompanied by downward shifts in the majority spin channel. Obviously, this energetic sensitivity is a problem where one is interested in resolving an energy scale of ∼ 150 meV.
Another related issue is demonstrated in Figure 5 by considering the Mn magnetic moment and cell magnetization as U 3d Mn is varied. Strikingly, the magnetizations are severely overestimated compared to the experimental green-shaded region once a small positive value of U 3d Mn is introduced. In fact, the overestimation of magnetic moments with advanced DFAs or DFT+U is a well-known issue for metals displaying itinerant magnetism [2,[65][66][67]. Previous studies [65,66] showed that increasingly advanced DFAs in Jacob's ladder of DFA ranks [68] such as GGA PBE [51], meta-GGA SCAN [69], and hybrid PBE0 [52] tend to increasingly overestimate the magnetic moments in metallic ferromagnetic systems such as bcc Fe, hcp Co, and fcc Ni. Interestingly, even LDA could slightly overestimate the magnetic moments, such as in fcc Ni [66].
In Figure 5, it is evident that LDA+U requires a small negative U 3d Mn value to agree with both experimental Mn atomic moment and overall cell magnetization. We note that utilization of a negative U was suggested [70,71] and applied [72][73][74][75][76] before, especially on superconducting systems. Physically, this means that Hund's exchange J dominates the on-site Coulomb repulsion U , resulting in a minor negative U eff = U −J value. The result is more delocalization of orbitals and a tendency to pair electrons, which produce lower magnetic moments. The opposite extreme can be seen when using HF and hybrid DFA for metals, which over-localize the electrons and tend to unpair the electrons, resulting in overestimated magnetic moments. Even though DFT+(U > 0) physics is different than hybrid DFAs, the effect is a similar localization of orbitals [70]. In this case, the proper delocalization could be achieved by a minor negative U eff or a Hund's J with a magnitude larger than U . In addition, previous studies showed that the employed double-counting scheme in DFT+U could significantly change the predicted magnetic moments [2,77]. Ultimately, the key goal is to achieve a proper (de)localization of electrons and hence the correct magnetization driven by an ap- propriate tendency of electrons to pair. The difficulties of predicting magnetic moments using advanced DFAs raise the issue of whether single-reference QMC calculations can predict the correct magnetization in a complex correlated metal such as TbMn 6 Sn 6 .

E. Single-Reference QMC
As kinetic energy sampling via k-mesh integration is crucial in metals [78][79][80][81], we use canonical twistaveraging (CTA) and a supercell with two formula units ( Figure 1) in all TbMn 6 Sn 6 QMC calculations. To assess the validity of single-reference QMC, we consider results from two different QMC methods. These methods are fixed-node/fixed-phase DMC (captures most correlations), and σ 2 → 0 extrapolated VMC (VMC σ 2 →0 extrap , estimate of converged dynamic correlations) using singlereference Slater-Jastrow type trial wave functions. The VMC σ 2 →0 extrap method uses two data points with the same determinantal form to extrapolate to the σ 2 = 0 limit as shown in Figure 6 for TbO molecule where nearly-exact energies can be obtained via UCCSD(T)/CBS. Specifically, energies from the Slater trial wave function and Slater-Jastrow trial wave function with one-, two-, and three-body terms (SJ 123 ) were used in the extrapolation (black squares). In principle, as the variance goes to zero, the VMC energy should recover the exact correlation energy due to zero-variance principal [49]. However, since the Jastrow factor only captures the dynamic correlations [82], we do not expect this estimator to recover the exact correlation, as the trial wave function is improved only in the symmetric term (due to Jastrow), while the antisymmetric term (due to Slater) remains as an approximation. Although one can not strictly distinguish between static and dynamic correlations, the above separation is clearthe Jastrow term is the source of 'dynamic correlations', while the fermionic term is the source of 'static correlations'. In this sense, the VMC σ 2 →0 extrap method is expected to give the energy for the 'perfect Jastrow'. We would like to note that this is only an approximation, as σ 2 = 0 should only occur for the exact ground state. Provided that the trial wave function Ψ T has a significant overlap with the exact ground state only, one can show that the VMC energy must linearly depend on the variance for σ 2 → 0 [83]: In practice, it can be difficult to achieve near-ideal linear extrapolations, for instance, as seen in multi-determinant extrapolations of high-pressure hydrogen [50]. We investigated the convergence behavior of one-and two-body Jastrow (J 12 (n)) and also three-body Jastrow (J 3 (n)) as the number of optimizable parameters n is increased. A nonlinear convergence was observed for J 12 (n) in all cases of Figure 6. For large n, J 12 (n) reaches a plateau where the energy is converged with respect to n. A linear behavior was observed when using these converged J 12 functions and including J 3 (n) terms with various n values. The simple two-point extrapolations using the Slater trial wave function (free of Jastrow optimization imperfections) and the SJ 123 with the highest feasible optimizable parameters (the best Jastrow we could optimize) proved to be the most robust and consistent way to compare energies. The obtained VMC σ 2 →0 extrap value would therefore represent the VMC energy of a SJ (123...∞) , i.e, Jastrow with up to M -body interactions where M → ∞.
In Figure 6, even though VMC σ 2 →0 extrap estimator energy is below the single-reference DMC, it is still significantly above the estimated exact energy since the estimator is mainly probing for the effect of dynamic correlations due to single-reference form. We observed similar plots with the same energy ordering for the other molecules considered in this work; see Supplementary Figure 14.
The only exception to this was the Sn atom, where the VMC σ 2 →0 extrap estimated energy goes slightly below the UCCSD(T)/CBS value (by about 10 −4 Ha). This is not surprising since the fixed-node bias in the Sn atom is very small (∼ 2% of the correlation energy), and thus we might expect that the static correlation is captured well by single-reference and VMC σ 2 →0 extrap should give close to exact energies. More information about VMC energy extrapolations and extended details of QMC methods are given in Supplementary Note 4.
In Figure 7(a-b), we plot the bulk energies with the above-described methods where the orbitals are generated from DFA = {LDA, PBE}. Specifically, LDA DFT calculations obtain magnetizations close to the experiments, so it is interesting to see whether QMC predicts the correct magnetization when the LDA orbitals are reused for various QMC magnetizations. Namely, in Figure 7(a-b), the same orbitals are used, but a different overall magnetization M s is constructed in QMC calculations. Since QMC methods are variational, the lowest energy states correspond to the physical predictions from single reference QMC [84,85]. As shown in these plots, the lowest energy for both estimators corresponds to ∼ 15 µ B , approximately twice the correct value and well outside of the experimental region. This overestimation seems independent of the employed orbitals, as LDA and PBE results follow similar trends. In addition, extrapolating the dynamic correlations to the zero-variance limit does not change the energy curve, signifying a necessity for either significantly improved orbitals or including static correlations via multi-reference wavefunctions.
Next, we explore fixing the QMC cell magnetization to agree with experiments (at M s = 7 µ B ) and varying the orbitals by the introduction of Hubbard +U (Figure 7(c-f)). In Figures 7(c-d), this is done for unconstrained DFT (M s from SCF) and magnetizationconstrained DFT (M s = 7 µ B ) using an effective Hubbard U as described in Section II D. In Figure 7(c), the DMC method obtains the lowest energy for unconstrained LDA at U 3d Mn = 1 eV. The magnetization in this reference is severely overestimated in DFT with cell magnetization of ≈ 12 µ B and Mn moment of 3.38 µ B , see Figure 5. On the other hand, the lowest DMC energy for the constrained reference occurs at a much larger value of U 3d Mn = 2 eV. This is because in the constrained case, the Mn moments increase much more slowly as U 3d Mn is increased (Figure 5(a)) requiring a larger value of U 3d Mn . This suggests that the magnetization of Mn atoms dictates the QMC energy.
In Figure 7(d), similar plots are shown using VMC σ 2 →0 extrap method. The results are qualitatively similar to DMC in Figure 7(c), but the minima are shifted by about ∆U 3d Mn = 0.5 eV closer to zero. However, the corresponding DFT magnetizations are still considerably overestimated and outside the experimental measurements (see Figure 5 for corresponding U 3d Mn values). This shows that including the full dynamic correlations can improve the prediction of magnetic properties, but it is not enough for an accurate estimation in this case.
Finally, we explored the effect of J in the DFT+U + J formalism [60] in Figures 7(e-f). We consider LDA with U 3d Mn = 1.0 eV, which results in the lowest energy in DMC for this DFA (Figure 7(c)). In DMC/SJ 123 method, Figure 7(e), the energies increase as J is increased. This is because the effect of J is to discourage high-spin states, and we see a similar effect as in Figure 7(c).
The VMC σ→0 extrap in Figure 7(e) method predicts a value of J 3d Mn = 0.6 eV; however, the corresponding Mn moment of 2.61 µ B is still slightly overestimated relative to the experimental value of 2.39(8) µ B [24] (see Supplementary  Figure 8).
It is evident from Figure 7(c-f) that constraining the QMC cell magnetization seems to alleviate the overestimation problem. This is especially pronounced in the VMC σ→0 extrap case where full dynamic correlations are ex-pected to be captured. However, the lowest energies using constrained QMC cell magnetizations in Figure 7(c-f) are still much higher than the unrestricted, M s = 15 µ B case in Figure 7(a). This is true within each method, DMC/SJ 123 and VMC σ 2 →0 extrap . These calculations show that single-reference QMC overestimates the magnetic moments similar to meta-GGA SCAN and hybrid PBE0 DFAs. This overestimation trend persists even when the dynamic correlations are extrapolated to the zerovariance values, Figure 7(d). This points to a significant deficiency in static correlations, which possibly stem from a few sources: 1. Imperfect one-particle orbitals. Within the singlereference framework, the nodal surface is fully determined by the single-particle inputs. In case of the presence of static/multi-reference correlations, this can bias the expectations. Natural orbitals (NOs) from CI calculations can provide a better orbital set for improving the description within the single-reference model [86].
2. Lack of multi-reference trial wave function. Even with high-quality orbitals such as NOs, it is possible that multi-reference wave functions are needed to predict the magnetic moments correctly. This would not be surprising because systems with Mn elements are known to display multi-reference characters [87]. For instance, a previous study of bulk MnO using CCSD found improvement in band gaps and magnetic moments on the overestimated values of UHF [45].
An additional example where a similar problem occurs is the W atom. Its ground state is [Xe]4f 14 5d 4 6s 2 ( 5 D) configuration, while the first excited state configuration is a higher spin [Xe]4f 14 5d 5 6s 1 ( 7 S) [88]. Fixed-phase spin-orbit DMC, with minimal expansion of determinants, incorrectly predicts the high-spin state 7 S as the ground state [89]. However, as the trial wave function is expanded with more determinants, the correct 5 D ground state is recovered in agreement with the experiments [89].
3. Finite-size effects. In this work, we used a supercell with two formula units and canonical twistaveraging where each twist is occupied with the same M s . Although CTA should reach the same TDL energy as grand-canonical twist averaging (GCTA) for the same magnetization, it could have a considerable impact on finite supercell sizes. This overlaps with the first point above as it is a onebody effect. However, it has a different origin since the change is across-twist rather than within-twist.
Although the above effects can be present, it is out of the scope of this work to study TDL convergence and to apply multi-reference approaches here since we are interested in the single-reference model and its accuracy limit for these types of materials. Therefore, we leave these aspects for future study. Intuitively, the QMC overestimation issue could be explained as follows. Mn element in the atomic limit favors the 6 S 5/2 high-spin state [88], and this localized atomic energy scale could dominate the QMC energy in the metallic solid. The energy due to bonding is partially captured in the single-reference case; however, it is not properly balanced with the localized energy scale, which skews the optimal energy to the Mn atomic high-spin state. Therefore, a balanced trial wave function would also need to capture the multi-reference effects stemming from crystal field and orbital hybridization so as to counterweight the pronounced atomic high-spin limit. Indeed, in the HF method, where there are no correlations by definition, the Mn moment is close to the fully unpaired limit of 5 µ B , see Table 1. Compared to insulators, the interactions must be significantly screened in metals [90], and HF or hybrid DFAs are well-known to result in poor screening [91,92], and other unphysical features in partially-filled band metals [93,94]. We note that in the thermodynamic limit (TDL) of full-CI (FCI), the proper screening must be achieved in metallic systems due to the cancelation of static correlation effects and long-range interactions [92][93][94]. Remarkably, the approximate correlations in LDA seem to effectively capture the screening and multi-reference effects quite well, judging by the obtained magnetizations in Table 1. Indeed, a few previous studies found that LDA could mimic the long-range correlations due to multireference effects [95][96][97], which is facilitated by the larger self-interaction errors in LDA compared to hybrids.

F. DFT+U Results
In order to shed additional light on the intricate relationship of single-reference and magnetic order, we empirically adjust the magnetic moment of the Mn atom to the neutron diffraction experimental result of 2.39(8) µ B [24]. A recent muon spin rotation measurement, which is a highly powerful probe of local magnetism, also obtained an Mn moment of ∼ 2.4 µ B at low temperatures [6]. We note that using the Mn moments as a reference is more robust than comparing against the overall cell magnetization or Tb atomic moments due to the orbital magnetic moment contribution of the localized and well-screened 4f orbitals [34]. On the other hand, the 3d orbital magnetic moments are negligible due to orbital quenching [2,34], and the spin moments can be directly compared with the experiments. In addition, the Mn moments were shown to display much smaller spin fluctuations [7], making Mn more suitable to compare against experiments. In fact, using muon spin rotation experiments, Ref. [6] showed a critical slowing down of spin fluctuations below T * C1 ≃ 120 K. Further colling the systems below T C1 ≃ 20 K resulted in freezing the spin fluctuations into static patches of ideal out-of-plane FiM ordering, which persisted down to the lowest measured 1.7 K temperature. Therefore, although our calculations are T = 0 K and collinear FiM, they can be directly compared with the Mn moments obtained from T < T C1 ≃ 20 K without considering the effects of spin fluctuations which are seen in some other frustrated kagome materials [10,98,99]. We find that employing LDA DFA and Löwdin population analysis, a small effective value U 3d Mn ≈ −0.5 eV is required to approximately match with the neutron diffraction experimental value of 2.39 (8) Mn (−0.5 eV) results in Mn moment of 2.366 µ B ). This value of U 3d Mn is qualitatively corroborated by Figure 3(a), where Mn seems to favor more delocalized orbitals for the given magnetization, and previous studies showing that even LDA, which presumably produces the smallest orbital localization, shows a tendency to overestimate the magnetic moments [65,66].
Having established an appropriate DFA for the study of TbMn 6 Sn 6 , we now plot the bulk band structure for the FiM phase using LDA+U 3d Mn with U 3d Mn = −0.5 eV. In Figure 8(a), the bands are plotted with scalar relativistic pseudopotentials, namely, with averaged spin-orbit interactions. We find that the DC2 energetic level (E DC2 ) with respect to E F shifts down to about E DC2 = 120 meV when compared with LDA (Figure 4(a)). This is in excellent agreement with the experimental result of E exp DC = 130(4) meV from tunneling and quasiparticle scattering along the bulk crystal edge direction [1]. In addition, the inclusion of explicit SOC results in the opening of the Chern gap in DC2 and DC3 seen at 120 meV above E F , Figure 8(b). The Chern gap for DC2 is ∆ = 25 meV, which is also in reasonable agreement with the experimental value of ∆ exp = 34(2) meV.
Our calculations show that DC2 in bulk is close to E F , which could affect the experimental AHE measurements. An important distinction between the bulk theory and experimental transport measurements is that the experiments were carried out in a cleaved bulk terminated with an Mn surface. Specifically, the Mn-terminated surface showed a significant modulation in the dI/dV line map at B = 9 T , interpreted as a Landau quantization signature. On the other hand, the TbSn-terminated surface showed an almost homogenous dI/dV line map [1]. This was attributed to a possible geometry reconstruction due to dangling Tb 3+ bonds [1]. The Mn termination could also modify the band structure near the surface. To explore this, we carried out DFT+U calculations of a TbMn 12 Sn 10 non-stoichiometric semi-infinite slab terminated with Mn layers. In other words, the slab is periodic along in-plane directions, while a large vacuum is inserted in both out-of-plane ends, which are terminated with Mn layers, see Figure 9. To demonstrate that the Mn surface effects can be appropriately modeled by this slab, we calculated the atomic moments and charges in these two settings. Table 2 provides these values for Tb and Mn in various environments. We observed that the Tb spin moments and charges are largely unchanged for bulk vs. slab, indicating that the relevant changes occur only in the surface and sub-surface Mn layers as well as Sn layers in between. The sub-surface Mn shows a moderate change in the moments. However, a significant change occurs near the Mn surface, where we see much larger Mn moments of 3.33 µ B , in quantitative agreement with previous findings [2]. This shows that the surface effects can be reliably modeled by the constructed slab.  Figure 8(c) and 8(d) show the slab band structures. Compared to bulk bands, the notable differences are DC3 shifting up from ∼ 120 meV to ∼ 275 meV, while DC2 shifting down very close to E F . The shift in DC2 is likely due to charge transfer from the surface Mn layer to the subsurface Mn layer (see Table 2 for atomic charges). In Figure 8(c), where the SOC is not explicitly included, the DC2 is ∼ 15 meV below the E F . However, once the SOC is included, the Chern gap of ∆ = 26 meV opens up, and the Fermi level lies within this gap (Figure 8(d)). In reality, the Fermi level in the cleaved bulk might or might not exactly lie in the Chern gap of 34(2) meV, as the ab-initio methods used here combined with a level of experimental input could result in larger systematic bias than 34 meV. However, the change in E DC2 going from the bulk to the surface, ∆E DC2 = 135 meV must be a more robust quantity, as this is an energy difference, and we expect it to be much less sensitive to the quality of DFA or trial wave functions.
To obtain further insights into the DCs' origin, we investigated their orbital characters using projections onto the atomic Mn orbitals. In the bulk, there are two DCs near 120 meV (Figure 8(a), DC2 and DC3). DC3 has a heavy d zx and d zy character, while DC2 has a d xy and d x 2 −y 2 character (see Supplementary Figure  FIG in-plane views are shown. The large red spheres are Tb atoms, medium blue spheres are Mn atoms, and small gray spheres are Sn atoms. The slab is Mn-terminated due to the observed intense modulation in the experimental dI/dV line map [1]. The vacuum length between the periodic slabs is ≈ 13.5 Ang. 9). Since DC3 is significantly above the E F in the slab (∼ 275 meV), we focus on DC2 with d xy , d x 2 −y 2 character and plot the band structures of the surface, subsurface, and bulk Mn layers projected onto atomic d xy , d x 2 −y 2 orbitals in Figure 10. The subsurface Mn bands look similar to the bulk Mn bands, but the DC2 appears closer to E F . On the other hand, the surface bands look quite different, and the d xy , d x 2 −y 2 DC2 signal near the E F is very weak. We note that a DC2 near E F can still be seen on the surface Mn layer, albeit with strong d z 2 , d zx , and d zy characters (Supplementary Figure 11). Namely, the DC2 near E F in the slab (Figure 8(c)) is split between the surface and subsurface DCs in the d-orbital manifold with different characters. Perhaps it is not surprising that the surface DC2 takes on a substantial z character since it now faces the vacuum. Assuming that STM transport measurements involved conduction mainly in x, y directions (in-plane), the above DC2 split suggests that the experimental signal mainly consists of the bulk and subsurface layers with a stronger d xy , d x 2 −y 2 character, and a weaker signal coming from the surface Mn d zx , d zy orbitals. This is corroborated by the fact that the STM Landau fan diagram obtained from an Mn-terminated bulk shows the Chern gap at 130(4) meV above E F , not near E F as we see for the surface.
Another aspect of this material studied by Ref. [2] was the k z dispersion of the bulk TbMn 6 Sn 6 . Specifically, the authors mentioned that the k z -dispersion was too high in the bulk TbMn 6 Sn 6 , which is not desired for a 2D model. Indeed, we also find this to be true in bulk, and the band structure shows drastic changes as the k z is slowly varied (Supplementary Figure 12). In fact, the DC2 and DC3 occur only in k z = 0 when looked at increments of ∆k z = 0.1 in reciprocal lattice units and disappears for higher values. The 'flat' band in k z = 0 also gradually disperses for high k z values. The non-ideal 2D nature of the bulk can be seen in the Fermi surface as well (see Supplementary Figure 5). This demonstrates the additional challenge in realizing the Chern magnetism in the bulk Mn layers and motivates the experimental synthesis of TbMn 6 Sn 6 thin films to eliminate the heavy k zdependance and to move the DC2 closer to E F . Finally, besides E DC2 , we have also calculated the Dirac velocity (v D ) to compare with experiments. To obtain v D , we use the Dirac dispersion in the presence of a Chern gap due to SOC [1]: The fit of Equation 2 to the upper band of DC2 in Figure  8(b) provides v D = 2.25 × 10 5 m/s (see Supplementary Figure 4). This is only in qualitative agreement with the experimental value of v exp D = 4.2(3) × 10 5 m/s which was obtained by a fit to the Landau fan diagram [1]. In addition, the Dirac velocity extracted from tunneling data provided results close to v exp D = 4.2(3) × 10 5 m/s [1]. Therefore, LDA+U 3d Mn seems to genuinely underestimate v D (and relatedly, Fermi velocities). This v D underestimation is only mildly sensitive to the value of +U 3d Mn , as the increase of v D going from LDA+U 3d Mn (−0.5 eV) to LDA+U 3d Mn (+0.5 eV) is only ∆v D = 0.27 × 10 5 m/s. Namely, the main effect of +U 3d Mn in this material is to shift the band energies up/down while only moderately changing the dispersions (namely slopes or v D ). Therefore, we don't expect LDA+U 3d Mn (−0.5 eV) to obtain both E DC2 and v D correctly, nor do we expect it to produce accurate results for all other properties such as optimized geometries. We note that the v D underestimation is not only specific to TbMn 6 Sn 6 , and it has been observed in other materials such as graphene and Dirac semimetal Na 3 Bi, where LDA underestimates the Fermi velocity while hybrid DFT or GW methods show considerable improvements [100,101]. Ultimately, the inability to predict all properties using the same DFA shows the limitations of the single-particle picture, which supports our findings throughout this work.
To further investigate the nature of electron correlations and spin fluctuations in TbMn 6 Sn 6 , we carried out DFT+Dynamical Mean Field Theory (DMFT) calculations (see Supplementary Note 5 for DMFT methods). The preliminary results on the high-temperature paramagnetic (PM) phase of TbMn 6 Sn 6 show strong orbitaldependent electron correlations. We find that in the PM phase, the electron occupation of the Mn 3d-subspace to be 5.3 e − and the mean fluctuating local moments (< m z >) of Mn 3d electrons to be ∼ 3.7 µ B , which is in good agreement with a recent DFT+DMFT computation performed on YMn 6 Sn 6 [102], a sister compound in this family. To compute the fluctuating moment, we use where P i is the probability of the i th multiplet in the continuous time Monte Carlo impurity solver and |S z | i is the absolute value of the corresponding moment. The histogram plot (Supplementary Figure 16) shows that the most probable spin states are high-spin states with S z = 2.0, 2.5, and 1.5, indicating the correlation in the PM phase is likely due to Hund's rule coupling [103]. This prediction also agrees with other kagome metals, including YMn 6 Sn 6 [102,104]. Hund's metals were shown to display strong electron correlations [105], and this agrees with our QMC results, which show significant correlations due to multi-reference character. The underlying mechanism of the origin of ferromagnetism with the ordering of fluctuating moments is an important area for further studies and is beyond the scope of the present work.

III. DISCUSSION
This study has made advancements in three distinct but related directions. First, an accurate, correlationconsistent ECP was generated for the Tb element with a valence space of 4f 9 5s 2 5p 6 6s 2 . To the best of our knowledge, the accuracy of such a core-valence partitioning was not known for many-body methods such as CCSD(T) or DMC. In CCSD(T), we found that the low-lying atomic gaps and the binding energies of TbH 3 and TbO can be made chemically accurate (1 kcal/mol) with careful optimizations. In light of the obtained results, this seemingly too technical aspect of the effective core model is actually important for a clear delineation of subtle physical effects in the studied system. Single-reference DMC calculations showed that ≈ 12% of the correlation energy is missing in Tb, TbH 3 , and TbO. We believe the promising results obtained for the Tb element indicate that similarly accurate ccECPs could be generated for other rare-earth elements. Such a set of rare-earth ccECPs would open further possibilities to study the broader family of RMn 6 Sn 6 with significantly improved accuracy using many-body methods such as CCSD(T), CI, QMC, and also DFT with appropriate DFAs. In fact, the current progress in applying real-space QMC to f -element systems seems to be hindered by the unavailability of accurate enough rare-earth ECPs [28,106], as we were able to find only a few QMC studies in the literature [107][108][109][110]. Therefore, we hope this work will motivate further studies in this avenue.
Second, recent studies [65,66] showed that meta-GGA such as SCAN, hybrid DFT such as PBE0, and DFT+U methods overestimate the magnetic moments in simple elemental metals such as bcc Fe, hcp Co, fcc Ni, and other transition metals. Here, we observe a similar overestimation in Mn magnetic moments when using DFT+U . Importantly, we show that single-reference QMC calculations also severely overestimate the transition metal moments similar to the previously shown [65,66] advanced DFAs such as SCAN and PBE0. The overestimation persists even when the dynamic correlations are extrapolated to the zero-variance limit, suggesting that a multireference treatment is needed to predict the magnetic moments correctly. These results have broader implications for the ability of commonly used DFAs, such as LDA, to effectively capture the multi-reference effects and for the origins of effective weak interactions in metallic systems. A study of static correlation effects is possible for such systems with small primitive cells as multi-reference methods such as CCSD(T) and CI mature for use in periodic boundary conditions. We leave this aspect of the study to future work and hope that these findings will stimulate further research in this direction.
Finally, using a combination of ab-initio and neutron diffraction experimental results, we reveal key insights about TbMn 6 Sn 6 bulk and non-stoichiometric thin film limit. We show that DC2 is only ∼ 120 meV above E F in bulk, in agreement with experiments [1]. However, the realization of Chern magnetism in bulk is compli-cated by the heavy k z -dispersion as suggested previously [2,7]. We show that in the slab, the DC2 shifts even closer to E F . This motivates further experimental studies of TbMn 6 Sn 6 in the thin film limit to eliminate the k z -dispersion and to probe for the possible realization of Chern magnetism. Viable pathways to experimentally synthesize defect-free [111] thin films in a controlled manner to realize the Chern phase remains to be shown. Our work thus illuminates the aspects of TbMn 6 Sn 6 both from theoretical and experimental viewpoints and opens the door for future studies of this exciting material.

A. DFT+U Methods
DFT+U calculations were carried using quantum espresso package [61][62][63]. DFT+U is well-known to converge to meta-stable states depending on the given initial guess for the Hubbard occupation matrix. To overcome this issue, we used a method similar to the ramping method described in Ref. [112]. Specifically, first, we converged a DFT+U calculation with a very small U value, such as 10 −16 eV. Then, the converged charge density, orbitals, and Hubbard occupation matrices were directly provided as an initial guess for the desired value of U , skipping the adiabatic ramping of U . On a few occasions, we tested that this method indeed results in lower energies than the default values for Hubbard occupations.

B. QMC Methods
QMC calculations were carried out using the qmcpack package [113,114]. Most calculations were driven by the nexus automation tool [115]. All bulk QMC calculations use canonical twist-averaging (CTA) with [8×4×8] twists to sample the kinetic energy. In this approach, every twist is charge-neutral. In addition, every twist (θ) is enforced to have the same cell magnetization: namely, every twist has the same magnetization as the twist-averaged magnetization of the cell. At each twist, the orbitals are occupied using the lowest Kohn-Sham eigenvalues in each spin channel. All molecular and bulk Jastrow functions were optimized using energy minimization in VMC. For bulk calculations, the Jastrow was optimized at k = Γ twist and reused for other twists. This optimization was carried out for LDA orbitals at M s = 7 µ B and reused for all other DFAs and M s values. This was done to obtain improved energy differences originating from the determinantal part of the trial wave function.
The input files, output files, and supporting data generated in this work are published in Materials Data Facility [134,135] and can be found in Ref. [136]. The data is also available from the corresponding author upon reasonable request. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.
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 to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

VI. COMPETING INTERESTS
The authors declare no competing financial or nonfinancial interests.
VII. AUTHOR CONTRIBUTIONS P.G. conceived and supervised the project. A.A. carried out the majority of calculations and wrote the manuscript. S.M. carried out the DFT+eDMFT calculations and wrote the related text. L.M. and J.T.K. aided with methodology approaches and writing the manuscript.