Nb$_3$Cl$_8$: A Prototypical Layered Mott-Hubbard Insulator

The Hubbard model provides an idealized description of electronic correlations in solids. Despite its simplicity, the model features a competition between several different phases that have made it one of the most studied systems in theoretical physics. Real materials usually deviate from the ideal of the Hubbard model in several ways, but the monolayer of Nb$_3$Cl$_8$ has recently appeared as a potentially optimal candidate for the realization of such a single-orbital Hubbard model. Here we show how this single orbital Hubbard model can be indeed constructed within a"molecular"rather than atomic basis set using ab initio constrained random phase approximation calculations. This way, we provide the essential ingredients to connect experimental reality with ab initio material descriptions and correlated electron theory, which clarifies that monolayer Nb$_3$Cl$_8$ is a Mott insulator with a gap of about 1 to 1.2eV depending on its dielectric environment. By comparing with an atomistic three-orbital model, we show that the single molecular orbital description is indeed adequate. We furthermore comment on the expected electronic and magnetic structure of the compound and show that the Mott insulating state survives in the low-temperature and bulk phases of the material.


I. INTRODUCTION
The first theoretical studies of correlation effects in solids date back as early as 1934 with studies on the so-called polar model [1].In this model, the crystal was represented as a periodic array of single-orbital sites such that each site can be empty, double occupied, or single occupied with a spin up or down electron, and various types of hopping and interaction processes were taken into account [2,3].Further simplification resulted in the appearance of the Hubbard model [4][5][6][7][8] in the 1960s, in which only basic processes of single-electron hopping and on-site Coulomb repulsion are taken into account.Although the Hamiltonian of this model is simple, the physics that emerges from the competition between interaction and kinetic hopping, or between localization and delocalization, is extremely rich and complicated.As a function of temperature, lattice structure, and number of electrons, the phase diagram of the model is believed to contain metallic and insulating phases, magnetic ordering, and more exotic phases such as unconventional superconductivity and charge-density waves [9].In general, the model is too difficult to be solved exactly in the thermodynamic limit (except in one [10] and infinite dimensions [11,12]), but many efficient computational strategies have been devised, and some regions of the phase space are well-understood nowadays [9,13,14].
Given its theoretical importance and computational difficulty, experimental realizations of the Hubbard model have long been sought since then nature itself would solve the model for us.Some unconventional superconductors, including cuprates [15,16] and nickelates [17,18], as well as Na x CoO 2 [19] have been described using the single-orbital Hubbard model.However, in these cases, most simplified single-orbital models are controversial, and there are indications that multiorbitals models are needed [9].Single-orbital models seem to be, nevertheless, adequate to describe transi-tion metal dichalcogenides undergoing charge-densitywave transitions, with one low-energy orbital per supercell [20].Artificial lattices, in the form of ad-atoms on surfaces [21,22] or atoms trapped in optical lattices [23] are other prominent examples.However, an actual solid that features the ideal realization of Hubbard model physics in the unit cell is of course even more desired.There are several difficulties associated with finding such material.First of all, most materials have several bands relatively close to the Fermi level, which casts doubt on a single-orbital low-energy description.Secondly, background screening in the solid should be so efficient that the Coulomb interactions between electrons on different sites are effectively zero, which is unrealistic for many two-dimensional materials [24].Thirdly, the Coulomb interaction between electrons on the same site should be strong enough to find the desired strong correlation effects.
It has been proposed that monolayer Nb 3 Cl 8 [25] is a suitable candidate for a true Hubbard material [26,27].At the DFT level, it has a single well-isolated band crossing the Fermi level.At the same time, experimental angle-resolved photo emission observations [26,28] and transport measurements [29] show a gap in bulk structures, which could be an indication of strong correlation effects.However, until now, no first-principles calculations of the strength of the local and non-local Coulomb matrix elements have been performed.Here, we use the constrained Random Phase Approximation (cRPA) [30] to calculate the relevant partially-screened Coulomb matrix elements and show that the local Hubbard interaction is indeed sufficient to open a Mott gap in the monolayer.Furthermore, by comparing models with three and one orbitals per unit cell, we demonstrate that the singleorbital model is indeed applicable.We provide all parameters, including interactions and hopping matrix elements, as necessary for many-body calculations of the material's electronic properties and investigate the de- rived models on the level of the Hubbard-I approximations.Furthermore, we show that the bulk structure can also be described by a slightly more involved two-orbital Hubbard model and investigate it in various experimentally observed high-and low-temperature structures.

II. RESULTS
The high-temperature lattice structure of monolayer Nb 3 Cl 8 is shown in Fig. 1.The Nb atoms form a distorted kagome lattice, with small (red) and big (blue) triangles.The distortion is 16.6%, i.e., the sides of the small triangle are 16.6% shorter and those of the big triangle 16.6% longer than for an undistorted kagome lattice, in which both would be a/2 = 3.36 Å.Every unit cell contains three Nb atoms, such that the unit cell can be chosen to contain exactly one complete small (red) triangle, which we refer to hereafter as a trimer.At low-temperature, further distortions take places [31][32][33], reducing the symmetry, which we discuss together with the bulk structures in more detail below.
Fig. 2(a) shows the electronic band structure at the DFT level, i.e., without considering strong electronic correlations.There are three Nb d t 2g bands close to the Fermi level, of which one is crossing the Fermi level, while the other two are slightly above.At Γ these states can be described by their 2a 1 and 2e symmetries, respectively [26,31,32].Additional t 2g and e g bands are above and below (green and blue).Finally, a block of Cl p bands can be found below −2.7 eV.A Wannier construction for the 2a 1 and 2e t 2g bands provides three equivalent orbitals ψ 1,2,3 each localized at one of the three Nb atoms in the trimer, as shown in Fig. 2(c)-(e).In each case, the orbital lobes have a particular orientation with respect to the triangle of atoms, as indicated in Fig. 2(b).
An isolated trimer with only internal hopping t 0 > 0 between the three atomic Wannier orbitals has three single-particle states with energies −2t 0 (2a 1 ) and twice-TABLE I. Single molecular orbital Hubbard model parameters for monolayer Nb3Cl8.The upper panel depicts the triangular lattice with lattice site indices.The lower panel lists the distance r between the trimer lattice sites, their single-particle hopping t, local and non-local Coulomb interactions U and effective magnetic exchange interactions J (s) = −2(t (s) ) 2 /U * .All parameters are given in meV.degenerate t 0 (2e) as obtained from diagonalization of the isolated trimer Hamiltonian.Then, adding a hopping t t 0 between trimers leads to one low-energy band and a pair of higher energy bands, each with a bandwidth of order t and a gap of 3t 0 .The strong distortion of the kagome lattice thus allows us to describe the 2a 1 and 2e bands observed in Fig. 2(a) by a triangular lattice tightbinding model in a basis of molecular (trimer) orbitals.
For an isolated trimer, the lowest molecular orbital wave function is a symmetric combination ψ 0 = 1 √ 3 (ψ 1 + ψ 2 + ψ 3 ) of the atomic wave functions, as shown in Fig. 2(f).Correspondingly, a Wannier construction for the lowest 2a 1 band using only one orbital per trimer yields such a molecular orbital, as shown in Fig. 2(g).This confirms that the trimer molecular orbital picture is a simple explanation of the low-energy single-particle electronic structure.
In Tables I and II, we summarize the hopping parameters corresponding to the single molecular orbital (MO) and three-band atomic orbital (AO) Wannier constructions, respectively.The magnitudes of the hopping parameters are consistent with the picture of weakly hybridized trimers on a triangular lattice (see Methods for more details).

A. Coulomb interaction using constrained Random Phase Approximation
Having established two Wannier orbital basis sets, we come to our main result: the partially screened Coulomb interactions, which we calculate from first prin- ciples using the constrained Random Phase Approximation (cRPA) [30] (see Methods).To this end, the electronic structure is divided into a target space and a rest space.We calculate the effective Coulomb interaction in the target space by considering all screening processes involving rest space electrons according to the Random Phase Approximation.The cRPA is an approximation for the screening by the rest-space electrons.This approximation should be applicable when the rest space has a large gap [34], and the two spaces can be disentangled cleanly.Here, the symmetry of the distorted kagome lattice greatly helps the disentanglement.Furthermore, the rest space has a gap of 2.5 and 2 eV for the atomicAO and MO target spaces, respectively.This is sufficiently large compared to the bandwidths of the individual bands to have confidence in the cRPA accuracy.In this sense, Nb 3 Cl 8 is an extraordinarily promising candidate for realizing a most simple Hubbard model.
In the present case, there are two relevant choices for the target space: all of the three lowest 2a 1 and 2e states spanned by three atomic Nb d orbitals, or just the lowest metallic 2a 1 band, best described by the lowest MO.We have performed cRPA calculations for both cases, and the resulting Coulomb matrix elements U (R) are given in Tables I and II.For the onsite R = 0 Coulomb interactions we find U ≈ 1.9 eV in the MO model and U 0000 ≈ 2.8 eV (U 0011 ≈ 1.8 eV) in the atomic model.U and U 0000 are obviously not identical.There are two origins for these differences, one purely geometrical and one physical.The MO is constructed as a linear combination of atomic orbitals, meaning the Coulomb tensor should also be transformed to a new basis involving all intra-trimer Coulomb matrix elements yielding U ≈ 1/3 U 0000 +2/3 U 0011 (see Methods).Given the numbers in Table II, more than half of the molecular Coulomb interaction comes from the interatomic interaction U 0011 within the trimer.The transformation to the MO reduces the Coulomb matrix element, since the molecular Wannier orbital, as depicted in Fig. 2(g), is more spread out than the individual atomic Wannier orbitals in Fig. (c)-(e).In addition to this geometrical effect, the Coulomb interactions in the MO model are further reduced because of the additional screening provided by the two 2e bands that are integrated out.
Fig. 3 shows additional information about the Coulomb interaction in the MO.Formally, the cRPA produces a frequency-dependent interaction [30,35] U (ω).In Nb 3 Cl 8 , the partially-screened Coulomb interaction is, however, well approximated as frequency-independent in the range ω ∈ [0, 10] eV, and only reaches the asymptotic, bare value around ω = 30 eV.This shows that retardation effects are unimportant at the relevant electronic energy scales [35].Next to these weak retardation effects, we find significant non-local values U (R) for R > 0, as shown in Fig. 3(b), depicting the 1/R tail of the static Coulomb interaction.In this figure, we also show how the non-local interaction of the molecular orbital is affected by a dielectric encapsulation of the monolayer (see Methods for details on the calculations).From this, we see that in the extrapolated fully free-standing ε = 1 case U (R = 0) ≈ 2 eV and U (R = 1) ≈ 0.9 eV, which can be reduced by an environmental screening of ε = 10 to approx.1.4 and 0.4 eV, respectively.U (R) is thus rather susceptible to the environment.However, the effective local interaction U * = U (R = 0) − U (R = 1), which is the adequate quantity to use in purely local TABLE II.Atomic Hubbard model parameters for monolayer Nb3Cl8.The upper panel depicts the atomic positions within the triangular lattice.The lower panel lists the distance r between the atomic trimer lattice sites, their singleparticle hopping t, local and non-local density-density Uiijj and Hund's exchange Uijji Coulomb interactions.All parameters are given in meV.single-band Hubbard models [24], is barely affected by the substrate screening as a result of similar screening effects to U (R = 0) and U (R = 1).This is illustrated in Fig. 3(c), which shows that U * (ε) reduces from about 1.15 to only 1 eV upon increasing the environmental dielectric screening from ε = 1 to 10.

B. Electronic correlations
Solving the single-orbital triangular lattice Hubbard model is a formidable task [36] which requires advanced computational techniques and is beyond the scope of this paper.Our model derivations result in consistent parameter sets summarized in Table I, which should be used for such calculations.Nevertheless, inspecting these model parameters given by our Wannier construction and cRPA calculations, we observe that the Coulomb interaction is substantially larger than the bandwidth of the partially filled low-energy band, U/|9t (1) | ≈ 10 for the molecular model.Thus, strong coupling techniques can be used to get a first impression of the correlation effects.The so-called Hubbard-I approximation, a multi-band generalization [37] of the approach introduced by Hubbard [6], was suggested as an adequate approximation by default in the narrow-band limit [37,38] such that we will exploit it here.
In Fig. 4(a) and (c), we show the spectral functions of the atomic and MO models, i.e., with 3 and 1 orbitals per trimer, calculated using the Hubbard-I approximation, taking into account the intratrimer interactions only.For the MO model, we take into account U (R = 0) only, while for the three-orbital model, we include the full trimer-local Coulomb tensor U ijkl where i, j, k, l are restricted to nearest-neighbour Nb positions on a single trimer.We also show in Fig. 4(d) results for the MO model using U * to account for long-range intertrimer Coulomb interactions.In all cases, the interaction leads to a significant frequency dependence in the local Hubbard-I self-energy, which splits the partially filled 2a 1 band into an upper and lower Hubbard band, with a gap given by the effective interaction strength acting on the lowest half-filled molecular orbital.For the atomic model, we have projected onto the trimer eigen- state ψ 0 to show that this is indeed the low-energy state.This comparison clearly proves that the system is a Mott insulator in the Hubbard-I approximation, regardless of the model that is used.
Comparing the Mott gaps in the three models, we find noticeable differences.The U * model obviously has a smaller gap since it uses a reduced Coulomb interaction to account for inter-trimer screening.This effect is not included at all in the approximate Hubbard-I treatment of the other models.Comparing the gaps in the atomic three-orbital and the MO model, we find a small reduction of the gap in the MO model.Both approaches differ in their diagrammatic content with regard to the two higher-energy trimer states, and an exact match is, therefore, not expected.The three-orbital model in the Hubbard-I approximation contains purely intratrimer effects in an exact way.In contrast, the molecular model contains both intra-and intertrimer screening effects, but only on the (c)RPA level, while the correlation effects are subsequently taken into account using Hubbard-I, i.e., intratrimer only.The quantitative similarity is a sign that the appearance of the gap is largely caused by intratrimer correlation physics.
Since it starts with three bands, the atomic model obviously has additional spectral weight at higher energies as well.Fig. 4 shows that this spectral weight predominantly remains in the tight-binding bands, with a limited amount of spectral weight transferred to localized excitations at higher energy.
Finally, in all cases, we observe a significant bandwidth renormalization in the lower and upper Hubbard bands as compared to the non-interacting 2a 1 bandwidth.This is a result of the frequency dependence of the dynamical self-energy and is in agreement with recently obtained angle-resolved photo emission spectra [26].

C. Magnetic Properties
Given the strong Coulomb repulsion and the filling of one electron per trimer, there is a tendency to form single S = 1  2 magnetic moments on every trimer.In this strong coupling limit, a simple estimate for the magnetic exchange interactions is given by the Hubbard Stratonovich transformation J (s) = t (s) /U * , and the resulting values are given in Table II up to third nearest neighbors.Fig. 5 shows possible magnetic structures based on these interactions, including a ferromagnetic, striped anti-ferromagnetic, and 120 • anti-ferromagnetic ordering together with their energy densities of the exchange interaction given by Here, i = 1, • • • , N i stands for the magnetic moments in the magnetic unit cell, and neighbors within the same distance |r (s) ij | characterized by shell index s.From these structures, the 120 • antiferromagnetic structure has the lowest energy, suggesting it adequately represents the magnetic ground state.
Experimentally, a paramagnetic state was observed at high temperatures for bulk samples [39].In quasi-twodimensional magnets, magnetic order exists only at temperatures much lower than typical energy of exchange interactions due to strong thermal spin fluctuations.At intermediate temperatures, one can expect paramagnetism with strong short-range order, which can be described using, for example, self-consistent spin-wave theory [40].initio models, we showed that monolayer Nb 3 Cl 8 in the high-temperature (HT) phase is a Mott insulator deep in the Mott regime as a result of the integer one-electron occupation of each trimer and the strong trimer-local Coulomb interactions relative to the small inter-trimer hopping.Below we further show that this Mott benhaviour also survives in the distorted low-temperature (LT) phases as well as within bulk stacks.Previous studies suggested using U ≈ 1 eV as an ad hoc choice together with a Wannierized single-band to describe bulk structures [26,27].Our first-principles cRPA study shows that this value is indeed reasonable, once the system is mapped onto an effective Hubbard model using the renormalized U * = U (R = 0) − U (R = 1) [24], which effectively takes the non-local Coulomb interaction into account.We find slightly larger values of U * ≈ 1.15 eV for an isolated monolayer, with a reduction towards 1 eV for monolayers in a dielectric environment [41,42].
Screening by the environment is thus relatively inefficient due to the intra-and inter-trimer Coulomb interaction being both screened in approximately the same way.
Based on the effective single-band Hubbard model describing ML Nb 3 Cl 8 , we furthermore showed that a 120 •antiferromagnetic structure forms between the magnetic moments localized on each trimer at the ground state.

A. Spin-Orbit Coupling
In the distorted kagome crystal structure yielding the correlated trimer lattice, Nb is formally in an [Nb 3 ] 8+ state.As the atomic spin-orbit coupling of Nb 2+ and Nb 3+ is below 100 meV [43] and thus significantly smaller than the local Coulomb of U 0000 ≈ 2.8 eV, spin-orbit coupling can be treated as a perturbation to the Mott insulating state here and will not change qualitatively the drawn conclusions.The interplay of spin-orbit coupling with the crystal field could lead however to subtleties in the magnetic anisotropy [44].

B. Bulk Structures
In the bulk structure, the details of the stacking between the layers affect the electronic properties via interlayer hybridization.In the orthorhombic lattice structure of the HT phase, the correlated trimers in the different layers are laterally shifted as visualized in Fig. 6(a).The resulting DFT band structure is shown in Fig. 7(a).It is similar to the ML HT band structure with a doubled amount of bands, which are mildly affected by interlayer hybridization.We now find two half-filled bands around the Fermi level with finite dispersion in k z direction.The corresponding atomic and MO models now host 6 Nb d states and 2 MOs, respectively.The nearest-neighbour intra-layer hopping matrix element t ≈ 25 meV of the corresponding MO model is similar to the one in the ML.The nearest-neighbour inter-layer hopping t ⊥ 1,2 ≈ −15 to 17 meV is on a similar order, such that the inter-layer hybridization does not qualitatively change the small-band width character of the metallic bands around the Fermi level.The direction dependence of t ⊥ 1,2 is a result of the Cl positions, one of which can be either above or below the trimer.
Next to the kinetic parts of the model Hamiltonian, the Coulomb interaction matrix elements are also affected by the bulk structure due to enhanced screening.In the MO model, the local Coulomb interaction on each trimer is reduced from U ≈ 1.9 eV in the ML to U ≈ 1.5 eV in the bulk structure.Similar trends hold for the long-range interactions within the individual layers.Additionally, there are now inter-layer Coulomb matrix elements, of which only the density-density ones are of non-vanishing magnitude.In the MO model, these are ca.360 meV.Given the half-filled flat bands and still significant trimerlocal Coulomb interaction matrix elements, it is thus reasonable to expect that the Mott behaviour in the HT bulk phase survives.The corresponding results, taking only trimer-local interaction terms into account, are presented in Fig. 7 for both models.In the atomistic model, we find a Mott gap of about 1.6 eV, and in the MO model of about 1.5 eV.We note that these calculations do not take non-local Coulomb interactions into account, which can further influence the Mott gap as U * did in the ML model.For the bulk structures, it is at the moment, however not clear how to construct U * , since U (R = 1) is different along the in-plane and out-of-plane directions.

C. Low-Temperature Distorted Phases
For low temperatures, below 100 K, bulk Nb 3 Cl 8 is known to undergo a structural phase transition accompanied by a paramagnetic to non-magnetic transition [31][32][33]39].Two mildly different low-temperature structures with C2/m [31] and R3 [32,33] point group symmetry have been experimentally observed.The corresponding DFT band structures are shown in Fig. 7.The main difference between the C2/m and R3 structures is a de- formation within the trimers, yielding slightly elongated trimers in the C2/m structure, and modifications to the trimer sizes between the different layers in the R3 structure, as depicted in Fig. 6.As a result, the degeneracies between the two higher-lying trimer states are lifted in both LT structures.The main difference between these LT and HT bulk structures is, however, the relative shift between the layers.In the LT structures, trimers from different layers are partially significantly closer to each other than in the HT structure, c.f. Fig. 6.This enhances the inter-layer hybridization, which leads to a splitting of the two metallic bands around the Fermi level, which yields a full but small gap between these states.On the DFT level, it is important to note that none of the symmetry breakings in the C2/m and R3 structures yield significant charge re-distributions between the Nb atoms of the trimers.
The molecular orbital description is furthermore still valid and allows to quantify of the nearest-neighbor intraand inter-layer hoppings yielding t ≈ 15meV, t ⊥ 1 ≈ 127 meV, t ⊥ 2 ≈ 3 meV, respectively for the R3 structure.Notably, t ⊥ 2 is much smaller than t ⊥ 1 , which is in turn much larger than t .The LT structures is thus best described as weakly hybridized stacks of bilayer Nb 3 Cl 8 , whereby the two layers forming the bilayers are strongly hybridized.The Coulomb interaction matrix elements are only barely affected by the relative shift between the two layers of the LT bulk unit cell yielding U ≈ 1.5 eV in the R3 structures, respectively.A similar description holds for the C2/m structure, however, with reduced symmetry.All model parameters are given in Table IV.
We have also applied the Hubbard-I approximation to the bulk structures.As shown in Ref. [45], the Hubbard-I approach corresponds to the formal first-order-in hopping correction to the Green's function, assuming that the interaction is local, which leads to vanishing vertex corrections.This assumption is still clearly valid in the case of the LT structures as the interactions between the trimers, as well as the hopping between them, is much smaller than the intertrimer Coulomb interaction.The only other factor which could yield vertex corrections is correlated hopping [1,2,46].Explicit calculations show here, however, that these contributions to the hopping are small: The largest U iiij element with i and j being located on different trimers, which could be responsible for correlated hopping terms, is on the order of just 3 meV and can thus be safely neglected.In Fig. 7, we, therefore, present the corresponding interacting spectral functions for both atomic and MO models.Similar to the situations studied above, we find Mott gaps on the order of 1.6 eV driven by the dynamical properties of the Hubbard-I self-energy.Due to the small splitting of the two bands around the Fermi level, we see a similarly small splitting in the lower and upper Hubbard bands as well.We again note, that these calculations do not include the effects of non-local Coulomb interaction terms.
It is worthwhile to note that the LDA+U approximation, popular due to its simplicity, is much less accurate than the Hubbard-I approximation for narrow-band systems such as Nb 3 Cl 8 .As discussed in detail for elemental rare-earth and their compounds [38,47], Hartree-Fock calculations utilized within LDA+U [6] cannot account for the significant bandwidth renormalizations, which are captured by Hubbard-I calculations and which have been observed experimentally in Nb 3 Cl 8 [26,28].

IV. CONCLUSIONS
Based on detailed ab initio down folding calculations, we were able to derive various minimal (generalized) Hubbard models for Nb 3 Cl 8 in several structures.For all of them, the derived model parameters, together with the Hubbard-I approximation suggest that Nb 3 Cl 8 is a Mott insulator, independent of the details of the low-or high-temperature bulk or monolayer structures.For the monolayer structures, we were able to derive the most simple single-band Hubbard model.For bulk structures in all known lattice structures, we derived two-orbital Hubbard models with long-range Coulomb interactions.
The sister compound Nb 3 Br 8 has been suggested to be an "obstructed atomic insulator" [48] already at the DFT level, based on symmetry analysis.Our results for Nb 3 Cl 8 suggest that correlation effects are also important in Nb 3 Br 8 and that these would give a substantial contribution to any experimentally observed gap.To quantify this, cRPA calculations for Nb 3 Br 8 are needed.
Non-local Coulomb interactions are substantial in this layered compound.For the monolayer, it is possible to use an effective Hubbard model with a modified Hubbard interaction to account for this [24].We need to stress that the role of non-local Coulomb interactions in the bulk structure should be investigated in more detail, since the effective Hubbard model approach assumes that the nonlocal Coulomb interaction is the same in all directions, which is not the case in this layered compound.
The analysis of electronic correlations in this work (using Hubbard-I) is restricted to the undoped case, with TABLE III.Wyckoff positions for bulk Nb3Cl8 characterized by P 3m1, R3, and C2/m space groups.The data is taking from Refs.[31,32] and modified using allowed symmetry operations in such a way that the transition from HT phase (P 3m1) to LT phase (R3 or C2/m) is clear, as demonstrated in Fig. 6 one electron per trimer.Doping away from this integer filling will destroy the Mott insulating phase and thus requires a more advanced many-body treatment of the arising correlation phenomena.The model parameters derived here make this possible and thus form a promising basis to study the iconic Hubbard model and its physics hand-in-hand between theory and experiment.
V. METHODS A. Electronic structure calculations

Cyrstal structures
We study the electronic and magnetic properties of Nb 3 Cl 8 using the unit cells of P 3m1, R3, and C2/m space groups [31,32].These unit cells have the following bases (unit cell vectors): with a 1 = 6.74566Å and c 1 = 12.28056 Å, Crystal structures with primitive bases A R3(p) and A C2/m(p) , similarly to the structure with basis A P 3m1 , contain only 22 atoms.The atomic positions for the unit cells in P 3m1, R3, and C2/m space groups are taken from Refs.[31,32] and their values in a new basis are obtained as To get the electronic structures in terms of pure t 2g and e g states for a triangular lattice, the above basis vectors were rotated such that all lattice triangle sides in each layer form the cube face diagonals in a rotated frame, see Fig. 2(b).Monolayers are constructed from these bulk structures.

Wannier Basis Sets and Non-Interacting Hamiltonians
We utilize two different Wannier basis sets representing the lowest trimer molecular orbital alone and an atomistic one describing the three relevant Nb d orbitals on the trimer-corner positions, which we refer to as the molecular orbital (MO) and atomic orbial (AO) models, respectively.To this end, we start with conventional DFT calculations utilizing the Perdew-Burke-Ernzerhof GGA exchange-correlation functional [49] within a PAW basis [50,51] as implemented in the Vienna Ab initio Simulation Package (vasp) [52,53] using the structures given above.We use (20 × 20 × 1) and (8 × 8 × 4) k-meshes for monolayer and bulk calculations, respectively, as well as an energy cutoff of 350 eV and Methfessel-Paxton smearing of σ = 0.05 eV.
For the MO model, we project the flat Kohn-Sham bands around the Fermi level to an (initial) Wannier orbital with s symmetry located at the center of the trimer.For the atomistic model, we similarly project all Kohn-Sham states between −1.0 eV and +1.5 eV to three individual Nb-centered d orbitals.Afterward, we maximally localize the orbitals using the wannier90 package [54] and applying an inner (frozen) window including either only the flat bands around E F or all states between −1.0 eV and +1.5 eV.Using the resulting localized orbitals ψ i (r), we obtain the single-particle Wannier Hamiltonian from the hopping matrix elements The resulting Wannier models perfectly interpolate all Kohn-Sham states in their respective windows, as indicated in Fig. 2 and Fig. 7.

Coulomb Matrix Elements from Constrained Random Phase Approximation Calculations
The Coulomb interaction matrix elements U ijkl are evaluated within the Wannier orbital basis sets from above and using the constrained random phase approximation (cRPA) scheme [30] according to which utilizes the partially screened Coulomb interaction Here v is the bare Coulomb interaction and Π cRPA is the partial polarization as defined with Π target describing all polarization contributions from within the target sub-space as defined by the Wannier functions.This way, we avoid any double counting of screening channels in subsequent solutions of the derived generalized Hubbard models.Π cRPA correspondingly describes screening from all other states except the targetones including from a significant amount of empty states from the full Kohn-Sham basis.In detail, we use in total 160 bands and define Π cRPA using Kaltak's projector method as recently implemented in vasp [55].
From this we derive the full retarded and non-local rank-4 Coulomb tensor U ijkl (ω) for both models.Casula et al. showed that using U (ω = 0) instead of the fully retarded U (ω) is justified when renormalized hopping parameters are utilized [56].The corresponding renormalization factor was shown to scale with the characteristic cRPA plasmon frequency ω p .To this end, we show in Fig. 3 U (ω) of the MO orbital model in the HT ML structure.This indicates that ω p ≈ 25 eV is large such that we can safely neglect retardation effects here and use only U ijkl (ω = 0).
All monolayer calculations are performed within a supercell of 25 Å height.All tabulated and mentioned Coulomb matrix elements refer to these calculations.To extrapolate to the effectively free-standing monolayer and to take dielectric screening from the environment into account, we use our Wannier Function Continuum Electrostatics (WFCE) approach [41] as applied in Ref. [57].eigenvalues E 0 = 2t 0 , E ± = −t 0 (note that t 0 is negative, so E 0 < E ± ) and corresponding eigenvectors where ψ 1 are the atomic orbital wavefunctions and λ = exp(2πi/3).For the Nb 3 Cl 8 trimer, t 0 = −0.325eV, so the trimer gap 3|t 0 | ≈ 1 eV.
In the lattice of trimers, we have to take into account the hopping between trimers (i.e., along the blue bonds in Fig. 1).In the molecular orbital model, the hopping be-tween trimers is found to be t 1 ≈ 22 meV.In the atomic orbital model, the hopping from the corner of one trimer to the nearest atom on the neighboring trimer has an amplitude of 84 meV.Based only on this, one might have expected t 1 ≈ 84/3 = 28 meV, but sub-leading matrix elements (with alternating signs) are responsible for corrections.For a triangular lattice with nearest-neighbor hopping, the bandwidth is 9t 1 ≈ 200 meV.
For the Coulomb interaction, there are two effects that play a role when going from the atomic orbitals to the molecular orbitals.First, there is a purely geometric effect, since the molecular orbitals are linear combinations of the atomic orbitals, which leads to a transformation

FIG. 1 .
FIG. 1. Lattice structure of Nb3Cl8.(a) shows a side view, and (b) shows a top view.Red (blue) triangles indicate the first (second) nearest Nb neighbors, while brighter and dimmer spheres represent atoms (of Nb and Cl) in opposite layers.

FIG. 2 .
FIG. 2. Electronic structure and Wannier orbitals of the atomic and molecular orbital models.(a) Shows DFT and Wannierinterpolated electronic structures without spin polarization.Bands in red, green, and blue colors stand for Nb-t 1 2g , Nb-(t 2 2g +eg), and Cl p states, respectively.t 1 2g orbitals are schematically illustrated in (b) for each Nb atom (in rotated coordinates), while t 2 2g represents all other t2g orbitals.Opened and filled markers stand for the atomic (AO) and molecular (MO) orbital models, respectively.(c-e) Show Wannier atomic orbitals for each Nb atom, (f) their sum, and (g) the Wannier molecular orbital.

FIG. 3 .
FIG. 3. Coulomb interactions of the molecular single-orbital model calculated using the constrained Random Phase Approximation.(a) Frequency-dependence of the local interaction, (b) non-local static interactions for a monolayer in a dielectric environment.(c) Effective local Hubbard interaction U * = U (R = 0)−U (R = 1) as a function of the dielectric environmental screening.

FIG. 4 .
FIG. 4. Interacting spectral functions of Nb3Cl8 monolayer within the Hubbard-I approximation.(a) Atomic three-orbital model (AO), (b) its projection to the lowest molecular-orbital, (c) molecular single-orbital model (MO), (d) MO model based on U * , and (e) cross-section of those spectral functions at one selected momentum.The chosen chemical potential constraints the occupation to one electron per trimer and to align the upper Hubbard band.

E 1 FIG. 5 .
FIG. 5. Magnetic structures based derived from the singleorbital Hubbard model for monolayer Nb3Cl8.Corresponding energy densities of the exchange interaction are indicated in the Figure.

FIG. 6 .
FIG. 6. Crystal structures of bulk Nb3Cl8.(a) At high temperature (HT, space group P 3m1) and (b) and (c) at low temperatures (LT, space groups C2/m and R3, respectively).The on-top views depict layer 1, layer 2, layer 1/layer 2, and layer 2/layer 1 , respectively.Spheres in fade and bright colors (green and yellow) represent atoms (of Nb and Cl) in neighboring layers.Brown arrows indicate the interlayer couplings t ⊥ n (n = 1, 2) between the nearest trimers.Numbers around trimers are distances between the nearest Nb atoms in each layer (in Å).

FIG. 7 .
FIG. 7. (a)DFT electronic structure and (b) interacting spectral functions obtained using the Hubbard-I approximation for bulk Nb3Cl8.The top, middle, and bottom panels represent results in P 3m1 (HT phase), C2/m (LT phase), and R3 (LT phase) space groups, respectively.The electronic structures in (a) are weighted between Nb-t 1 2g (in red), Nb-t 2 2g (in green), and Nb-eg (in blue) states, respectively.The spectral functions in (b) stand for the atomic six-orbital model (first panels), molecular two-orbital model (second panels), and the cross-section of those spectral functions at one selected momentum (last panels).

B
. Trimer modelA single trimer can be described in tight-binding theory using the Hamiltonian Ĥ = t 0 .

TABLE IV .
Generalized Hubbard model parameters for bulk Nb3Cl8 in HT and LT phases for the two molecular orbital models.Tables list in-plane t parameters together with local and non-local Coulomb matrix elements.i and j are molecular orbital positions separated by distance |r s ij |/a = |(rj +Rij)−ri|/a (in units of the in-plane lattice parameter a), where Rij is a unit cell translation vector, and s is shell index (represents bonds with the same symmetries). hopping