On the Ferrimagnetic State of CrCl₂(pyz)₂

Abstract

Van der Waals layered ferromagnetic compounds with high two-dimensional electronic conductivity holds strong potential for quantum computing, future unconventional superconductors, catalysts, batteries, and fuel cells. We suggest a minimal theoretical model to understand the magnetic properties of the metal-organic framework CrCl₂(pyz)₂ (pyz=pyrazine). Using a Hubbard model we show that the groundstate is dominated by a specific configuration of delocalized electrons on the pyz sites with a ferrimagnetic coupling to the localized spins on the Cr sites. This model suggests a magnetic moment of $2\mu_{B}$ which is remarkably close to the experimental value of $1.8\mu_{B}$ [K. S. Pedersen et al., Nat. Chem. 10, 1056-1061 (2018)]. From Weiss mean-field theory we predict a weak ferromagnetic Cr-Cr coupling of $\approx 0.9$ meV. This is consolidated by second order perturbation theory of the RKKY interaction yielding a Cr-Cr coupling of $\approx 5$ meV. Understanding the interactions in these types of compounds can facilitate designs of metal-organic compounds with tailored magnetic properties.

I Introduction

Metal-organic frameworks (MOFs) have been predicted to serve as a platform for the next generation of quantum technology [1, 2]. The diversity of organic ligands compared to metal ligands allow for enhanced tunability of electronic and magnetic properties. These compounds may be key to harvest renewed insights into and exploitation of quantum phenomena such as the Quantum Hall Effect [3, 4, 5], topological insulators [6] and superconductivity [7]. The putatively first organic-inorganic hybrid 2D MOF, CrCl₂(pyz)₂ (pyz = pyrazine), was successfully synthesized in 2018 [8]. CrCl₂(pyz)₂ exhibits both long-range magnetic order and high 2D electronic conductivity. The compound has a Curie temperature of 55 K, and its magnetization at low temperatures saturates to 1.8 $\mu_{B}$ in strong external fields. This indicates an intralayer ferrimagnetic coupling in agreement with predictions by density functional theory (DFT) [9, 10]. First principle calculations indicate that CrCl₂(pyz)₂, like other hybrid materials [11, 12, 13], can be exfoliated making it relevant for 2D applications [9]. The majority of 2D materials are non-magnetic in their pristine form [14]. A magnetic 2D material, like CrCl₂(pyz)₂, could have applications within spintronics [15], multiferroics [16, 17, 18] and magnetoelectronics [19]. This material could possibly even have applications within dark matter detection [20, 21, 22]. The sibling coordination solids VCl₂(pyz)₂ and TiCl₂(pyz)₂ display contrasting physical properties. In the former, divalent V ions facilitate an antiferromagnetically ordered ground state, while in the latter trivalent Ti ions favor paramagnetism and a correlated Fermi liquid state [23]. In principle, this opens the possibility for tuning through exotic phase transitions via doping. This may be compared to the magnetic phase transitions proposed for Mn_2-xCr₂Sb [24, 25] and Mn₂Sb_1-xBi_x [26]. Finally, the perovskite S₂FeMoO₆ may be highlighted as another material in which the coexistence of itinerant and localized electrons cause an ordered, in this case putatively ferromagnetic, ground state [27, 28].

A crucial step towards synthesizing materials with desired magnetic properties is to understand the mechanism yielding these properties. Here we provide a simple explanation of the ferrimagnetic coupling and explain why the magnetic moment of CrCl₂(pyz)₂ is approximately 2$\mu_{B}$. As such, we propose a minimal model for monolayer CrCl₂(pyz)₂ that captures the essential features of chromium couplings mediated by intermediate and itinerant pyrazine electrons.

First, we compute a tight-binding model using a Slater–Koster procedure to obtain hopping parameters. From this we propose an effective model relevant for the valence electrons. The ground state configuration of the model is investigated before a magnetic field is introduced and the magnetic moment is computed. Using mean-field theory an estimate of the Cr-Cr coupling is found. Finally, the exchange interaction between electrons on neighboring chromium and pyrazine sites is calculated and an expression of the indirect exchange interaction coupling between neighboring chromium sites is obtained using second order perturbation of the RKKY interaction. From this we predict an estimate of the coupling between the localized Cr spins and the itinerant pyz electrons.

II Tight Binding Model of Monolayer $\textbf{CrCl}_{\textbf{2}}\textbf{(pyz)}_{2}$

Monolayer CrCl₂(pyz)₂ has an approximately four-fold rotational in-plane symmetry. Each layer is coupled through Van der Waals interactions. The pyrazine rings are tilted with two alternating orientations [8], see FIG. S1 in the Supplementary Material. One unit cell contains two Cl atoms, one Cr atom, and two pyrazine rings, see a fragment of a unit cell in Fig. 1 a).

We compute an effective tight-binding model for monolayer CrCl₂(pyz)₂, as summarized in FIG. 1 b). The lattice is approximated as a square lattice in the Cr-pyz plane, and each pyrazine ring is modeled by a single site with two orbitals tilted by $45^{\circ}$ compared to the in-plane Cr lattice, see FIG. S1. Due to the tilting of the pyrazine sites, the Cr sites and the pyz sites form a Lieb lattice [29]. A study of the symmetries of the model complex CrCl₂(pyz)₄ using a group theoretical approach can be found in the Supplementary Material. The basis includes five d-orbitals on the Cr site, three p-orbitals on both Cl sites, and two p${}_{z^{\prime}}$ orbitals tilted in two different orientations on the two pyz sites, $\mathbf{\Phi}_{\textbf{k}}=\mathbf{\Phi}_{\textbf{k}}^{\text{Cr}}\otimes\mathbf{\Phi}_{\textbf{k}}^{\text{Cl}}\otimes\Phi_{\textbf{k}}^{\text{pyz}}$. Here $z^{\prime}$ refers to the direction perpendicular to the plane spanned by the pyz rings. The Hamiltonian

$$H=\sum_{\textbf{k}}\mathbf{\Phi}_{\mathbf{k}}^{\dagger}\hat{\mathcal{H}}(\textbf{k})\mathbf{\Phi}_{\textbf{k}}^{\phantom{\dagger}}$$ (1)

describes inter- and intraorbital hopping between the aforementioned orbitals, as encoded by the $15\crossproduct 15$ matrix $\hat{\mathcal{H}}$. The form of $\hat{\mathcal{H}}$ is listed in Eq. (VI) in the Supplementary Material. The basis is not spin-resolved. The elements in the Hamiltonian has been obtained using the Slater–Koster decomposition of the tight binding model [30, 31].

The band structure along the path $\mathbf{M}\rightarrow\mathbf{\Gamma}\rightarrow\mathbf{X}$ is displayed in FIG. 1 c). The Fermi Level is located at $0$ eV. This captures the qualitative, semi-metallic properties in agreement with first principle calculations [8]. The two empty d-bands of Cr are located above 1.5 eV, whereas the filled d-bands are located below the p-bands (not visible in the plot). This is also seen from first principle calculations of the PDOS [8] and could indicate a strong $\pi$-d hybridization. The flatness of the d-bands indicates a localization of the d-electrons in accordance with the behavior of a Mott insulator. The p-levels from the pyz sites are dispersive which indicates a delocalization of the p-electrons. The Cl levels are located far below the Fermi level and will be safely neglected in the minimal model. From first principle calculations the summed spin population on the pyrazine rings is calculated to $-0.68$ delocalized over four pyz sites of the model complex CrCl₂(pyz)₄ [8].

FIG. 1: a) Illustration of a fragment of one unit cell of CrCl₂(pyz)₂. b) Monolayer CrCl₂(pyz)₂ viewed along the Cl-Cr-Cl axis. The chromium and pyrazine sites are illustrated with black and white circles, respectively. Each pyz ring has been approximated by a single sites. Further, the unit cell is marked with a dashed square and the relevant lattice vectors, $\mathbf{a_{1}}$ and $\mathbf{a_{2}}$, are indicated. c) Band structure of monolayer CrCl₂(pyz)₂ obtained from tight binding along the path $\mathbf{M}\rightarrow\mathbf{\Gamma}\rightarrow\mathbf{X}$. Inset show the Brillouin zone and the relevant symmetry points.

In this paper the two pyrazine rings in a unit cell will be approximated as two single sites over which one electron delocalize. Further, Ref. [8] concludes that the total spin on the Cr site is $S=3/2$. For these reasons, the Cr sites will be assumed to have a localized spin $S=3/2$, and we assume one spin $S=1/2$ to be delocalized over the two pyz sites.

III Minimal Model

To shed light on the microscopic mechanism of the magnetic ground state of CrCl₂(pyz)₂ we investigate the model shown in FIG. 2. This is not identical to two unit cells of a single layer of CrCl₂(pyz)₂ since the pyz sites are moved with respect to the two chromium atoms. The model does, however, capture the competition between kinetic energy and exchange energy, which comprise the essential low-energy physics relevant for explaining the magnetic properties of CrCl₂(pyz)₂. In this model the two chromium atoms are assumed to be localized $S=3/2$ spins. There are two electrons that can delocalize over the four pyrazine sites. The Cl sites are omitted in this description, yet effectively taken into account by reducing the number of electrons on the Cr sites.

FIG. 2: Model containing two chromium sites (red circles labeled A and B) and four pyrazine sites (blue circles labeled 1-4). Each chromium spin couples only to three out of the four pyrazine sites (A to 1-3 and B to 2-4).

We include the exchange coupling between the chromium spins and the spins of the two pyrazine electrons. Moreover, the pyrazine electrons are allowed to hop between the pyrazine sites (indicated with dashed lines in FIG. 2). Lastly, we include an on-site repulsive Hubbard interaction, $U$, penalizing electrons residing on the same site. The effective Hamiltonian can thus be written as

$$\hat{\mathcal{H}}= -2J\hat{\vec{S}}_{A}\cdot\left(\hat{\vec{S}}_{1}+\hat{\vec{S}}_{2}+\hat{\vec{S}}_{3}\right)-2J\hat{\vec{S}}_{B}\cdot\left(\hat{\vec{S}}_{2}+\hat{\vec{S}}_{3}+\hat{\vec{S}}_{4}\right)$$

$$-t\sum\limits_{\sigma=\uparrow,\downarrow}\Big{(}\hat{c}_{1\sigma}^{\dagger}\hat{c}_{2\sigma}^{\phantom{\dagger}}+\hat{c}_{2\sigma}^{\dagger}\hat{c}_{3\sigma}^{\phantom{\dagger}}+\hat{c}_{1\sigma}^{\dagger}\hat{c}_{3\sigma}^{\phantom{\dagger}}$$

$$\quad\quad\quad\quad\,\,+\hat{c}_{2\sigma}^{\dagger}\hat{c}_{4\sigma}^{\phantom{\dagger}}+\hat{c}_{3\sigma}^{\dagger}\hat{c}_{4\sigma}^{\phantom{\dagger}}+\mathrm{h.c.}\Big{)}$$ (2)

$$+2U\sum\limits_{i=1}^{4}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow},$$

where $\hat{\vec{S}}_{A}$ and $\hat{\vec{S}}_{B}$ represent the spin of the two chromium atoms and $\hat{\vec{S}}_{i}$ with $i=1,2,3,4$ labeling the spin of the four pyrazine sites. Moreover, $\hat{c}_{i\sigma}^{\dagger}$ and $\hat{c}_{i\sigma}$ are the creation and annihilation operators of an electron on site $i$ and spin $\sigma$, respectively. $\hat{n}_{i\sigma}=\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}\phantom{\dagger}$ is the spin density on site $i$ with spin $\sigma$. $J$ is the exchange coupling and $t$ the hopping constant.

The spins of each chromium atom can take four $S_{z}$ values, i.e., $-3/2$, $-1/2$, $1/2$ and $3/2$. Details of the diagonalization are provided in the Supplementary Material. The basis is given by $|ms_{A},ms_{B}\rangle\otimes|\psi\rangle$, where $|\psi\rangle$ and $|ms_{A},ms_{B}\rangle$ are the given states for the pyrazine and chromium electrons, respectively. To invoke an antiferromagnetic coupling between the spin of the pyrazine electrons and the chromium spins, $J$ must be negative. The eigenstates are generally superpositions of 448 states but when all three variables, $J$, $t$ and $U$, are comparable in size the ground state is two-fold degenerate. Surprisingly, however, two basis states in these superpositions turns out to dominate, namely

$$|-3/2,-3/2\rangle\otimes|0\uparrow\uparrow 0\rangle\quad\text{and}$$ (3)

$$|3/2,3/2\rangle\otimes|0\downarrow\downarrow 0\rangle.$$

This is a robust feature of the model checked for the physically relevant parameter ranges of $t,U$ and $J$. For $t=U=-J=1$ the amplitudes of these states are $-0.747$, and is $\geq 0.5$ for all parameters considered. A two level model is exploited in the Supplementary Material which indicates that one unit cell acts as a spin-2 particle. For $t=U=-J=1$ the occupation probability of the basis states for one of these degenerate states is illustrated, see FIG. 3, clearly showing the dominance of one basis state.

FIG. 3: A sorted list of the occupation probability of each basis state for one of two degenerate ground states for $t=U=-J=1$. For the illustrated ground state the dominant basis state is $|-3/2,-3/2\rangle\otimes|0\uparrow\uparrow 0\rangle$. It is evident that the probability of all other basis states is small compared to the most dominant basis state. For symmetry reasons the plot of occupation probabilities for the other ground state (when sorted) looks identical, though the dominant basis state here is $|3/2,3/2\rangle\otimes|0\downarrow\downarrow 0\rangle$.

This result has a simple, intuitive explanation. The exchange interaction is maximal when both the pyrazine electrons and the chromium spins are antiparallel and maximal. Electrons located on sites 2 and 3 couple to both chromium spins, unlike electrons located on site 1 and 4 which only couple to one of the Cr spins. Thus, the exchange interaction is further maximized by the electrons being located at site 2 and 3. For $t=0$ this effect is further enhanced and the amplitudes of the states shown in Eq. (3) increases to $-0.845$. When $t\neq 0$ the electrons can gain energy by delocalizing and the amplitudes of the triplet states with $S_{z}=0$ and the singlet states increases. The two pyz electrons can delocalize and still avoid each other in two different configurations for a given spin configuration. One electron can delocalize over sites 1 and 2 while the other delocalizes over 3 and 4, see FIG. 4 a), or one electron delocalizes over sites 1 and 3 while the other one delocalizes over sites 2 and 4, see FIG. 4 b).

FIG. 4: The two different configurations of the spins of the pyrazine sites in the minimal model. One electron delocalizes over two sites with different exchange energy and in this way gain kinetic energy while still minimizing the on-site repulsion.

A magnetic field can be introduced by adding the term

$$\hat{\mathcal{H}}_{B}=g\mu_{B}\mathbf{B}\cdot\hat{\mathbf{S}}_{\mathrm{tot}}=g\mu_{B}B^{z}\hat{S}_{\mathrm{tot}}^{z},$$ (4)

to the Hamiltonian. We choose the magnetic field to point out of the plane. The expectation value of total spin in the minimal model is found to be $\langle S_{\mathrm{tot}}^{z}\rangle=\partial E/\partial B=2$. This is derived in the Supplementary Material. Since the model contains the same number of pyrazine rings and chromium atoms as two unit cells of CrCl₂(pyz)₂ this would correspond to one unit cell having $S=1$ which corresponds to a magnetic moment of

$$\mu_{\text{model}}=2\mu_{B}.$$ (5)

This is true both for $t=0$ and $t\neq 0$. This is strikingly similar to both experimental results [8] and DFT results [8, 9]. Experimentally, it is found that the saturation magnetization at high pressure and low temperature is $1.8\mu_{B}$ [8]. The slightly lower value can be ascribed to the canted radical spins in the real compound and the higher spin occupation on the pyrazines in DFT calculations (0.68 compared 0.5 in this model). Thus, we have here made a minimal model focusing on the relevant physics in the system that is quantitatively precise in its predictions on physical properties.

If we generalize this model to a layer of CrCl₂(pyz)₂ we would expect the ground state to look something like what is depicted in FIG. 5. Again, all Cr spins are localized and couple antiferromagnetically to the spins on the pyrazine electrons which have delocalized over two neighboring sites. Obviously, this can be configured in several different ways, but it will focus on one of these, since it will be very costly to change configuration.

FIG. 5: A sketch of one of the configurations the delocalized spins on the pyrazine ligands could choose.

Weiss mean-field theory for a Cr³⁺ square lattice that ignores the pyrazine rings yields a fairly weak ferromagnetic coupling, $J_{\text{Cr-Cr}}=\frac{3k_{B}T_{C}}{\mu_{0}zS(S+1)}\approx 0.9$ meV ($J_{\text{Cr-Cr}}/hc\approx 7$ cm^-1), see the Supplementary Material for further details. A more sophisticated estimate of the Cr-Cr coupling can be obtained by means of second order perturbation theory of the RKKY interaction mediated by intermediate pyrazine rings. This calculation yields $J_{\text{Cr-Cr}}\approx 0.15J^{2}/2t$, where $t$ is an effective square lattice hopping strength. Using a strong antiferromagnetic coupling of $J_{\text{Cr-pyz}}\equiv J\approx-0.25$ eV as suggested by ab initio calculations [8] and $t\approx 1$ eV (see FIG. 1), we predict that $J_{\text{Cr-Cr}}\approx 5$ meV $\left(J_{\text{Cr-Cr}}/hc\approx 40\leavevmode\nobreak\ \text{cm}^{-1}\right)$. Inserting an effective Cr spin of $S=1$, accounting for the presence of the pyz spins, brings the mean-field result closer to the perturbative estimate of the Cr-Cr coupling. Here, two independent calculations unify in a meV range prediction of the ferromagnetic Cr-Cr coupling. This coupling may be probed by neutron scattering techniques and spin wave experiments. Uniaxial strain experiments on CrCl₂(pyz)₂ may be employed to strain-tune $t$ and reveal potentially interesting insight into the importance of the itinerant electrons and the broken four-fold symmetry. Derivations for both expressions of $J_{\text{Cr-Cr}}$ is included in the Supplementary Material.

IV Conclusions

We have proposed both realistic and minimal models capable of explaining the ferrimagnetic state of monolayer CrCl₂(pyz)₂. From the Slater–Koster-derived band structure monolayer CrCl₂(pyz)₂ is found to exhibit an insulating ground state, in agreement with first principle calculations and the experimental results [8]. Furthermore, the d-bands are observed to be flat indicative of these electrons being localized whereas the pyrazine electrons are found to be very delocalized.

To investigate the magnetic properties of monolayer CrCl₂(pyz)₂ a simple model was proposed. Here, the electrons on the chromium sites are viewed as localized spins, whereas the electrons on the pyrazine sites can delocalize over a few sites. This suggests a magnetic moment of one unit cell to be $2\mu_{B}$ which corresponds to an antiferromagnetic interaction between the chromium spins and the pyrazine spins. Due to the uneven magnitudes of the spin on the Cr and pyz sites, this results in a ferrimagnetic ordering. The calculated magnetic moment is strikingly close to the experimental value of $1.8\mu_{B}$ [8]. Weiss mean-field theory estimates a direct ferromagnetic coupling between neighboring Cr sites of $J_{\text{Cr-Cr}}\approx 0.9$ meV $\left(J_{\text{Cr-Cr}}/hc\approx 7\leavevmode\nobreak\ \text{cm}^{-1}\right)$. Second order perturbation theory of the RKKY interaction between Cr and pyz yields a Cr-Cr coupling of the same order of magnitude, namely $J_{\text{Cr-Cr}}\approx 5$ meV. This Cr-Cr coupling may be probed by neutron scattering and spin wave experiments.

— Supplementary Material —

This is based on results that are previously reported in Ref. [32].

V Group Theory for Model Complex $\textbf{CrCl}_{\textbf{2}}\textbf{(pyz)}_{\textbf{4}}$

First, we investigate the energy levels and symmetries of the model complex CrCl₂(pyz)₄, see Fig. S1. For this to apply to monolayer CrCl₂(pyz)₂, the stoichiometry used in this section will be that of CrCl₂(pyz)₂. One unit cell of CrCl₂(pyz)₂ consists of one Cr ion, two Cl ions, and two pyrazine rings. In the main text and when investigating monolayer CrCl₂(pyz)₂ in the following this corresponds to only including pyz${}^{\text{1,r}}$ and pyz${}^{\text{2,t}}$ of the pyz sites in Fig. S1. Each chlorine atom reduce the number of electrons on the chromium atom by one. The two pyrazine rings reduce the Cr atom of an extra electron, such that three electrons remain in an unfilled d-shell on the Cr ion. If these spins occupy three distinct energy levels and align parallel to each other, to obey Hund’s rules, the spin on the Cr atoms is $S=3/2$.

FIG. S1: Illustration of the model complex CrCl₂(pyz)₄. The nomenclature of the six relevant atoms, the coordinate system in use, and the orientation of the pyrazine rings are illustrated.

As for Benzene, the 6 atoms in the pyrazine rings are sp² hybridized so we are left with one unhybridized p-orbital per site perpendicular to the plane of the ring. This will be referred to as the $z^{\prime}$ direction. Due to the presence of the N atoms in pyrazine, compared to Benzene, there are no degenerate states within the ring. The 6 energy levels for a pyrazine ring, where the H atoms are neglected and only the unhybridized p-orbitals are included, is shown in Fig. S2 a). In Fig. S2 b) the eigenstates for the $n=1,...,4$ levels are illustrated. Each eigenstate is visualized as circles plotted on top of the pyrazine ring with a radius proportional to the amplitude of the p-orbital on the given atom. The red and blue color corresponds to a positive and negative sign on the basis state, respectively. For the HOMO level ($n=3$) and the LUMO+1 level ($n=5$) the amplitude on the N atoms are equal to zero. These orbitals do not couple to the Cr atoms and will be discarded in the following calculations. The coupling between the carbon atoms in these states and the chromium atom will be negligibly small due to distance. The $n=1$ ($n=6$) level is so high (low) in energy that these can be neglected. Therefore, we only include the two levels $n=2$ and $n=4$, and each pyrazine ring is safely approximated as a single atom with p${}_{z^{\prime}}$-orbitals of two different amplitudes, $\alpha_{2}$ and $\alpha_{4}$, where $\alpha_{2}$ and $\alpha_{4}$ are the amplitudes on the N atoms of the $n=2$ and $n=4$ levels, respectively. The amplitude on the C atoms are neglected. Note, that the former is an odd linear combination whereas the latter is an even linear combination of the two N atoms.

FIG. S2: a) Energy levels of a pyrazine ring calculated from tight binding. We have excluded the 4 hydrogen atoms and only included the unhybridized p-orbitals on the 4 carbon atoms and 2 nitrogen atoms. b) The eigenstates of the 4 intermediate energy levels $n=2$ to $n=4$. The eigenstates are illustrated by circles with a radius proportional to the amplitude of the orbital occupying the given atom with a positive (red) or negative (blue) sign, respectively.

The pyrazine rings are canted relative to the Cl-Cr-Cl axis. In Ref. 8 the rotation angles are measured to $42.5^{\degree}$ and $43.9^{\degree}$ with respect to the crystallographic c direction for the pyrazine ring 1 and 2, respectively. This will be approximated as $45^{\degree}$ in both cases. This means that the $p_{z^{\prime}}$-orbitals of the pyrazine “atoms” are instead linear combinations of both $p_{x}$, $p_{y}$ and $p_{z}$ orbitals with respect to the global coordinate system, see FIG. S1.

The energy spectrum of model complex CrCl₂(pyz)₄ can be obtained using the tight binding model, which in second quantization can be written as

$$\hat{\mathcal{H}}_{t}=\xi\sum\limits_{i,\sigma}\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i,\sigma}-t\sum\limits_{\langle i,j\rangle}\sum\limits_{\sigma}\left(\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{j,\sigma}+h.c.\right).$$ (S1)

Here $\xi$ refer to the onsite energies, $t$ is the hopping integral and h.c. includes the Hermitian conjugate terms. $\hat{c}_{i,\sigma}^{\dagger}$ and $\hat{c}_{i,\sigma}$ are the creation and annihilation operators, respectively, for an electron on site $i$ with spin $\sigma$. Only nearest neighbor interactions are included. For simplicity the spin-indices will be omitted in the following.

The full Hamiltonian for a single molecule can be written in terms of the Slater–Koster elements in the basis

$$\Psi=\{ d_{xy},d_{yz},d_{zx},d_{x^{2}-y^{2}},d_{z^{2}},p_{x}^{t},p_{y}^{t},p_{z}^{t},p_{x}^{b},p_{y}^{b},p_{z}^{b},$$ (S2)

$$p^{1,r}_{z^{\prime}}(\alpha_{2}),p^{1,r}_{z^{\prime}}(\alpha_{4}),p^{2,t}_{z^{\prime}}(\alpha_{2}),p^{2,t}_{z^{\prime}}(\alpha_{4}),p^{1,l}_{z^{\prime}}(\alpha_{2}),p^{1,l}_{z^{\prime}}(\alpha_{4}),p^{2,b}_{z^{\prime}}(\alpha_{2}),p^{2,b}_{z^{\prime}}(\alpha_{4})\},$$

where the first five elements refer to the d-orbitals of the chromium atom followed by the three $p$-orbitals of both chlorine atoms (top and bottom), and lastly the p${}_{z^{\prime}}$-orbitals for the four pyrazine “atoms”. The full Hamiltonian can be written in matrix form as

$$\mathcal{H}=\begin{bmatrix}\mathcal{H}_{d}&V_{dp}\\ V_{pd}&\mathcal{H}_{p}\end{bmatrix}.$$ (S3)

$\mathcal{H}_{d}$ is a diagonal 5$\crossproduct$5 matrix with $\xi_{d}$ on the diagonal whereas $\mathcal{H}_{p}$ is a 14$\crossproduct$14 matrix with $\xi_{p}$ on the diagonal, where $\xi_{d}$ and $\xi_{p}$ are the onsite energies of the d- and p-orbitals, respectively. Thus, the onsite energies of all p-orbitals are assumed to be degenerate prior the interaction with the d-orbitals and similar for the d-orbitals. The two off-diagonal blocks of the Hamiltonian $V_{dp}$ and $V_{pd}$ couple the d-orbitals to the p-orbitals and vice versa.

V.1 Subspace of the d-orbitals

The Hamiltonian describing the system can be split up into a part that describes the d-orbitals, a part that describes the p-orbitals and then a term that describes the interaction between these orbitals,

$$\hat{\mathcal{H}}=\hat{\mathcal{H}}_{p}+\hat{\mathcal{H}}_{d}+\hat{V}.$$ (S4)

An eigenstate can be expressed as a sum of a p- and d-part of the state as

$$|\psi\rangle=|\psi_{p}\rangle+|\psi_{d}\rangle.$$ (S5)

The projection operators $\hat{\mathscr{P}}_{d}$ and $\hat{\mathscr{P}}_{p}$ for the d-subspace and the p-subspace, respectively, can be applied to rewrite the Hamiltonian. These projection operators must sum to identity, $\hat{\mathscr{P}}_{p}+\hat{\mathscr{P}}_{d}=\mathds{1}$, and are defined to project the eigenstate of the full Hamiltonian into the subspace of either the d- or the p-orbitals, as

$$\hat{\mathscr{P}}_{d}|\psi\rangle=E|\psi_{d}\rangle\quad,\quad\hat{\mathscr{P}}_{p}|\psi\rangle=E|\psi_{p}\rangle.$$ (S6)

Using these two operators yields

$$\hat{\mathcal{H}}_{d}|\psi_{d}\rangle+\hat{V}_{dp}|\psi_{p}\rangle =E|\psi_{d}\rangle,$$ (S7)

$$\hat{\mathcal{H}}_{p}|\psi_{p}\rangle+\hat{V}_{pd}|\psi_{d}\rangle =E|\psi_{p}\rangle.$$ (S8)

Here $V_{dp}$ and $V_{pd}$ refer to the 5$\crossproduct$14 and 14$\crossproduct$5 off-diagonal blocks of Eq. (S3), respectively. Isolating $|\psi_{p}\rangle$ in Eq. (S8) and inserting this in Eq. (S7) leads to the expression

$$\hat{\mathcal{H}}_{d}|\psi_{d}\rangle+\hat{V}_{dp}\frac{1}{E-\hat{\mathcal{H}}_{p}}\hat{V}_{pd}|\psi_{d}\rangle=E|\psi_{d}\rangle.$$ (S9)

Using that $|\psi_{d}\rangle$ and $|\psi_{p}\rangle$ are eigenstates of $\hat{\mathcal{H}}_{d}$ and $\hat{\mathcal{H}}_{p}$, respectively, and assuming all p-levels to be degenerate prior to the interaction, this can be written as

$$\xi_{d}|\psi_{d}\rangle+\frac{\xi^{2}}{E-\xi_{p}}|\psi_{d}\rangle=E|\psi_{d}\rangle.$$ (S10)

Here, we have introduced the eigenenergy of $V_{dp}V_{pd}$ as $\xi^{2}$. Setting $\xi_{d}=0$, $E$ is isolated as

$$E_{\pm}=\frac{\xi_{p}}{2}\left(1\pm\sqrt{1+\left(\frac{2\xi}{\xi_{p}}\right)^{2}}\right).$$ (S11)

Since we are seeking the bonding orbitals and $\xi_{p}$ is chosen to be positive, the $E_{+}$ solution is discarded. Knowing the eigenenergies of the 5$\crossproduct$5 Hamiltonian $V_{pd}^{\dagger}V_{pd}$, namely $\xi^{2}$, in the d-subspace allow us to solve Eq. (S11). Investigating this matrix analytically we find that the matrix can be written in a block diagonal in the basis $\{d_{yz},d_{zx},d_{xy},d_{x^{2}-y^{2}},d_{z^{2}}\}$,

$$V_{dp}V_{pd}=\begin{bmatrix}M_{11}&M_{12}&0&0&0\\ M_{12}&M_{22}&0&0&0\\ 0&0&M_{11}&M_{34}&M_{35}\\ 0&0&M_{34}&M_{44}&M_{45}\\ 0&0&M_{35}&M_{45}&M_{55}\\ \end{bmatrix},$$ (S12)

with

$$M_{11} =\frac{1}{2}|V_{pd\pi}|^{2}(4+\alpha),$$ (S13)

$$M_{12} =-\frac{\sqrt{3}}{2}V_{pd\sigma}V_{pd\pi}^{\dagger}\alpha,$$ (S14)

$$M_{22} =\frac{3}{2}|V_{pd\sigma}|^{2}\alpha,$$ (S15)

$$M_{34} =\frac{1}{4}V_{pd\pi}^{\dagger}\alpha\left(\sqrt{2}V_{pd\pi}-\sqrt{3}V_{pd\sigma}\right),$$ (S16)

$$M_{35} =-\frac{1}{4}V_{pd\pi}^{\dagger}\alpha\left(\sqrt{6}V_{pd\pi}+V_{pd\sigma}\right),$$ (S17)

$$M_{45} =\frac{1}{8}\left[4\sqrt{3}|V_{pd\sigma}|^{2}-\alpha\left(2\sqrt{3}V_{pd\pi}+\sqrt{2}V_{pd\sigma}\right)\left(V_{pd\pi}^{\dagger}-\sqrt{\frac{3}{2}}V_{pd\sigma}^{\dagger}\right)\right],$$ (S18)

$$M_{44} =\frac{1}{8}\left[12\sqrt{3}|V_{pd\sigma}|^{2}+\alpha\left(2V_{pd\pi}-\sqrt{6}V_{pd\sigma}\right)\left(V_{pd\pi}^{\dagger}-\sqrt{\frac{3}{2}}V_{pd\sigma}^{\dagger}\right)\right],$$ (S19)

$$M_{55} =\frac{1}{8}\left[4|V_{pd\sigma}|^{2}+\alpha\left(V_{pd\pi}^{\dagger}\left(6V_{pd\pi}+\sqrt{6}V_{pd\sigma}\right)+V_{pd\sigma}^{\dagger}(V_{pd\sigma}+\sqrt{6}V_{pd\pi})\right)\right],$$ (S20)

where $\alpha\equiv|\alpha_{2}|^{2}+|\alpha_{4}|^{2}$ has been defined. Here, $V_{ijb}$ is the $b=\sigma$ or $b=\pi$ bond integral between two orbitals $i=\text{(s, p, d, f, ...)}$ and $j=\text{(s, p, d, f, ...)}$. Note that $M_{11}=M_{33}$ by symmetry. The two blocks in Eq. (S12) are trivial to diagonalize. It is evident from Eq. (S12) that the $d_{yz}$ orbital only couples to the $d_{zx}$ orbital, and vice versa, and further that $d_{xy}$, $d_{x^{2}-y^{2}}$, and $d_{z^{2}}$ only couples to each other. To understand this, we elaborate on the symmetries of the molecule in the following.

The fact that the pyrazine rings are rotated by two different angles decreases the symmetry of the molecule $C_{\text{4h}}\rightarrow C_{\text{2h}}$. As always, the group describing the spatial symmetries of CrCl₂(pyz)₄ includes identity. Only three other operations leave the molecule invariant, namely inversion (i), a $\pi$ rotation of about the $z$-axis ($C_{2}(z)$) and a reflection in the $xy$ plane ($\sigma_{\mathrm{h}}$), see FIG. S1. These four symmetry-operations form the point group $C_{2\mathrm{h}}$. This group is Abelian and thus both the order of the group and the number of irreducible representations are equal to 4. The character table can now be constructed, see Table 1.

Table 1: Character table of the $C_{2\mathrm{h}}$ group.

$C_{2\mathrm{h}}$	$E$	$C_{2}(z)$	i	$\sigma_{\mathrm{h}}$	Quadratic Functions
$\mathrm{A}_{\mathrm{g}}$	1	1	1	1	$x^{2}$, $y^{2}$, $z^{2}$, $xy$
$\mathrm{B}_{\mathrm{g}}$	1	-1	1	-1	$yz$, $zx$
$\mathrm{A}_{\mathrm{u}}$	1	1	-1	-1
$\mathrm{B}_{\mathrm{u}}$	1	-1	-1	1

From the character table it is observed that $yz$ and $zx$ belong to the same representation and similarly for the last three $d$-orbitals. This explains the block diagonal form of Eq. (S12). Further, $\matrixelement{xy}{\mathcal{H}}{xy}=\matrixelement{yz}{\mathcal{H}}{yz}$ which is due to the chosen coordinate system conforming with the symmetry of the model complex such that these matrix elements are identical.

VI Tight Binding Model of Monolayer $\textbf{CrCl}_{\textbf{2}}\textbf{(pyz)}_{\textbf{2}}$

Monolayer CrCl₂(pyz)₂, see FIG. 1 b), is now modeled. Using the Rietveld refinement technique [33] the three lattice constants are in Ref. [8] found to be $a=6.90351(4)$ Å, $b=6.97713(4)$ Å and $c=10.82548(6)$ Å. The close affinity of $a$ and $b$ motivates the lattice to be approximated as a square lattice in the Cr-pyz plane. The relevant symmetry-points are thus approximated to be

$$\mathbf{\Gamma}=(0,0,0),\quad\mathbf{M}=\left(0,0,\frac{1}{2}\right),\quad\mathbf{X}=\left(-\frac{1}{2},0,\frac{1}{2}\right).$$ (S21)

These symmetry points are indicated on the Brillouin zone illustrated in FIG. 1 c).

A unit cell includes one chromium atom, two chlorine atoms and two pyrazine rings, which will be labeled 1 and 2. Assuming that the two pyrazine rings have the same onsite energy, and similarly for the chlorine atoms, a tight binding Hamiltonian for this system in real-space, including only nearest neighbor hopping between p- and d-orbitals, can be written as

$$\hat{\mathcal{H}}= \xi^{\mathrm{Cr}}\sum\limits_{\mathbf{r}_{i}}\sum\limits_{\mu=1}^{5}\hat{d}_{i,\mu}^{\dagger}\hat{d}_{i,\mu}^{\phantom{\dagger}}+\xi^{\mathrm{Cl}}\sum\limits_{\mathbf{r}_{i}}\sum\limits_{\nu=1}^{3}\left[(\hat{p}_{i,\nu}^{t})^{\dagger}\hat{p}_{i,\nu}^{t}+(\hat{p}_{i,\nu}^{b})^{\dagger}\hat{p}_{i,\nu}^{b}\right]$$

$$+\xi^{\mathrm{pyz}}\sum\limits_{\mathbf{r}_{i}}\left[(\hat{p}_{i,z^{\prime}}^{1})^{\dagger}\hat{p}_{i,z^{\prime}}^{1}+(\hat{p}_{i,z^{\prime}}^{2})^{\dagger}\hat{p}_{i,z^{\prime}}^{2}\right]$$

$$+\frac{1}{2}\sum\limits_{\mathbf{r}_{i}}\sum\limits_{\mathbf{r}_{j_{\pm}}=\mathbf{r}_{i}\pm\mathbf{a}_{1}}\sum\limits_{\mu=1}^{5}\left[t_{\mu,\alpha}^{1,r}\hat{d}_{i,\mu}^{\dagger}\hat{p}_{j_{+}}^{1}+t_{\mu,\alpha}^{1,l}\hat{d}_{i,\mu}^{\dagger}\hat{p}_{j_{-}}^{1}+\mathrm{h.c.}\right]$$ (S22)

$$+\frac{1}{2}\sum\limits_{\mathbf{r}_{i}}\sum\limits_{\mathbf{r}_{j_{\pm}}=\mathbf{r}_{i}\pm\mathbf{a}_{2}}\sum\limits_{\mu=1}^{5}\left[t_{\mu,\alpha}^{2,t}\hat{d}_{i,\mu}^{\dagger}\hat{p}_{j_{+}}^{2}+t_{\mu,\alpha}^{2,b}\hat{d}_{i,\mu}^{\dagger}\hat{p}_{j_{-}}^{2}+\mathrm{h.c.}\right]$$

$$+\sum\limits_{\mathbf{r}_{i}}\sum\limits_{\mathbf{r}_{j_{\pm}}=\mathbf{r}_{i}\pm\mathbf{a}_{3}}\sum\limits_{\mu=1}^{5}\left[t_{\mu,\nu}^{t}\hat{d}_{i,\mu}^{\dagger}\hat{p}_{j_{+}}^{t}+t_{\mu,\nu}^{b}\hat{d}_{i,\mu}^{\dagger}\hat{p}_{j_{-}}^{b}+\mathrm{h.c.}\right].$$

Here $\hat{d}_{i,\mu}^{\dagger}$ ($\hat{d}_{i,\mu}^{\phantom{\dagger}}$) creates (annihilates) an electron on the chromium site at $\mathbf{r}_{i}$ in unit cell $i$ and orbital $\mu$. Similarly, the operator $(\hat{p}_{j,\nu}^{t,b})^{\dagger}$ and $\hat{p}_{j,\nu}^{t,b}$ creates and annihilates an electron in the orbital $\nu$ at the chlorine atom at site $\mathbf{r}_{i}\pm\mathbf{a}_{3}$, respectively. Lastly, the operators $(\hat{p}_{j,z^{\prime}}^{1(2)})^{\dagger}$ and $\hat{p}_{j,z^{\prime}}^{1(2)}$ creates and annihilates an electron in a $p_{z^{\prime}}$ orbital on the pyrazine site 1 (2) at $\mathbf{r}_{i}+\mathbf{a}_{1}$ and $\mathbf{r}_{i}+\mathbf{a}_{2}$, respectively.

The three first terms of the Hamiltonian states the onsite energies of the d- and p-orbitals of the five atoms in the unit cell. The lower three lines of Eq. (S22) describe electrons hopping from a p-orbital to any nearest neighboring d-orbital. The factor $1/2$ compensate for double counting. For simplicity the spin indices have been omitted. The following on-site energies and hopping integrals are defined as

$$\xi^{\mathrm{Cr}}_{e_{g}} =\matrixelement{\mathbf{r}_{i};d_{\mu}}{\hat{\mathcal{H}}_{t}}{\mathbf{r}_{i};d_{\mu}},$$ (S23)

$$\xi^{\mathrm{Cl}} =\matrixelement{\mathbf{r}_{i}\pm\mathbf{a}_{3};p_{\nu}}{\hat{\mathcal{H}}_{t}}{\mathbf{r}_{i}\pm\mathbf{a}_{3};p_{\nu}},$$ (S24)

$$\xi^{\mathrm{pyz}} =\matrixelement{\mathbf{r}_{i}+\mathbf{a}_{1};p_{z^{\prime}}(\alpha^{\prime})}{\hat{\mathcal{H}}_{t}}{\mathbf{r}_{i}+\mathbf{a}_{1};p_{z^{\prime}}(\alpha^{\prime})},$$ (S25)

$$t_{\mu,\nu}^{t} =\matrixelement{\mathbf{r}_{i};d_{\mu}}{\hat{\mathcal{H}}_{t}}{\mathbf{r}_{i}+\mathbf{a}_{3};p_{\nu}},$$ (S26)

$$t_{\mu,\nu}^{b} =\matrixelement{\mathbf{r}_{i};d_{\mu}}{\hat{\mathcal{H}}_{t}}{\mathbf{r}_{i}-\mathbf{a}_{3};p_{\nu}},$$ (S27)

$$t_{\mu,z^{\prime}}^{1,r} =\matrixelement{\mathbf{r}_{i};d_{\mu}}{\hat{\mathcal{H}}_{t}}{\mathbf{r}_{i}+\mathbf{a}_{1};p_{z^{\prime}}(\alpha^{\prime})},$$ (S28)

$$t_{\mu,z^{\prime}}^{1,l} =\matrixelement{\mathbf{r}_{i};d_{\mu}}{\hat{\mathcal{H}}_{t}}{\mathbf{r}_{i}-\mathbf{a}_{1};p_{z^{\prime}}(\alpha^{\prime})},$$ (S29)

$$t_{\mu,z^{\prime}}^{2,t} =\matrixelement{\mathbf{r}_{i};d_{\mu}}{\hat{\mathcal{H}}_{t}}{\mathbf{r}_{i}+\mathbf{a}_{2};p_{z^{\prime}}(\alpha^{\prime})},$$ (S30)

$$t_{\mu,z^{\prime}}^{2,b} =\matrixelement{\mathbf{r}_{i};d_{\mu}}{\hat{\mathcal{H}}_{t}}{\mathbf{r}_{i}-\mathbf{a}_{2};p_{z^{\prime}}(\alpha^{\prime})}.$$ (S31)

$\hat{\mathcal{H}}_{t}$ includes hopping between the two relevant orbitals. Here, we have used that $\alpha^{\prime}=\alpha_{2},\alpha_{4}$. The numerical values for these can be obtained using the Slater–Koster decomposition [30, 31]. $t^{t}$ and $t^{b}$ are 3$\crossproduct$5 matrices that includes hopping between the 3 p-orbitals on the Cl atoms and the 5 d-orbitals on the Cr atoms. The last four lines in Eq. (S31) includes the 2$\crossproduct$5 hopping matrices between pyrazine orbitals and the chromium atoms that take the form

$$t_{\mu,z^{\prime}}^{1,r} =\begin{bmatrix}\alpha_{2}E_{z^{\prime},xy}&\alpha_{2}E_{z^{\prime},yz}&\alpha_{2}E_{z^{\prime},zx}&\alpha_{2}E_{z^{\prime},x^{2}-y^{2}}&\alpha_{2}E_{z^{\prime},z^{2}}\\ \alpha_{4}E_{z^{\prime},xy}&\alpha_{4}E_{z^{\prime},yz}&\alpha_{4}E_{z^{\prime},zx}&\alpha_{4}E_{z^{\prime},x^{2}-y^{2}}&\alpha_{4}E_{z^{\prime},z^{2}}\end{bmatrix}_{\mu,z^{\prime}},$$ (S32)

$$t_{\mu,z^{\prime}}^{1,l} =\begin{bmatrix}-\alpha_{2}E_{z^{\prime},xy}&-\alpha_{2}E_{z^{\prime},yz}&-\alpha_{2}E_{z^{\prime},zx}&-\alpha_{2}E_{z^{\prime},x^{2}-y^{2}}&-\alpha_{2}E_{z^{\prime},z^{2}}\\ \alpha_{4}E_{z^{\prime},xy}&\alpha_{4}E_{z^{\prime},yz}&\alpha_{4}E_{z^{\prime},zx}&\alpha_{4}E_{z^{\prime},x^{2}-y^{2}}&\alpha_{4}E_{z^{\prime},z^{2}}\end{bmatrix}_{\mu,z^{\prime}}.$$ (S33)

The above Slater–Koster elements in the global coordinate system for d-orbital $\mu$ on Cr and the p${}_{z^{\prime}}$-orbitals on pyz is given by

$$E_{z^{\prime},\mu}=\left[\frac{1}{4}\left(\sqrt{2}-2\right)E_{x,\mu}+2E_{y,\mu}+\left(2+\sqrt{2}\right)E_{z,\mu}\right].$$ (S34)

The bonding matrices $t^{2,t}$ and $t^{2,b}$ take on a similar form. The Hamiltonian can be formulated in reciprocal space as

$$\hat{\mathcal{H}}= \sum\limits_{\mathbf{k}}\bigg{\{}\xi^{\mathrm{Cr}}\sum\limits_{\mu=1}^{5}\hat{d}_{\mathbf{k},\mu}^{\dagger}\hat{d}_{\mathbf{k},\mu}^{\phantom{\dagger}}+\xi^{\mathrm{Cl}}\sum\limits_{\nu=1}^{3}\left[(\hat{p}_{\mathbf{k},\nu}^{t})^{\dagger}\hat{p}_{\mathbf{k},\nu}^{t}+(\hat{p}_{\mathbf{k},\nu}^{b})^{\dagger}\hat{p}_{\mathbf{k},\nu}^{b}\right]$$

$$+\xi^{\mathrm{pyz}}\left[(\hat{p}_{\mathbf{k},z^{\prime}}^{1})^{\dagger}\hat{p}_{\mathbf{k},z^{\prime}}^{1}+(\hat{p}_{\mathbf{k},z^{\prime}}^{2})^{\dagger}\hat{p}_{\mathbf{k},z^{\prime}}^{2}\right]$$

$$+\sum\limits_{\mu=1}^{5}\left[\hat{d}_{\mathbf{k},\mu}^{\dagger}\hat{p}_{\mathbf{k},z^{\prime}}^{1\phantom{\dagger}}\left(t_{\mu,z^{\prime}}^{1,r}e^{i\mathbf{k}\cdot\mathbf{a}_{1}}+t_{\mu,z^{\prime}}^{1,l}e^{-i\mathbf{k}\cdot\mathbf{a}_{1}}\right)+\mathrm{h.c.}\right]$$ (S35)

$$+\sum\limits_{\mu=1}^{5}\left[\hat{d}_{\mathbf{k},\mu}^{\dagger}\hat{p}_{\mathbf{k},z^{\prime}}^{2\phantom{\dagger}}\left(t_{\mu,z^{\prime}}^{2,t}e^{i\mathbf{k}\cdot\mathbf{a}_{2}}+t_{\mu,z^{\prime}}^{2,b}e^{-i\mathbf{k}\cdot\mathbf{a}_{2}}\right)+\mathrm{h.c.}\right]$$

$$+\sum\limits_{\mu=1}^{5}\sum\limits_{\nu=1}^{3}\left[t_{\mu,\nu}^{t}\hat{d}_{\mathbf{k},\mu}^{\dagger}\hat{p}_{\mathbf{k},\nu}^{t\phantom{\dagger}}e^{i\mathbf{k}\cdot\mathbf{a}_{3}}+t_{\mu,\nu}^{b}\hat{d}_{\mathbf{k},\mu}^{\dagger}\hat{p}_{\mathbf{k},\nu}^{b\phantom{\dagger}}e^{-i\mathbf{k}\cdot\mathbf{a}_{3}}+\mathrm{h.c.}\right]\bigg{\}}.$$

The eigenenergies of the $15\crossproduct 15$ matrix $\hat{\mathcal{H}}$ are obtained numerically for the path $\mathbf{M}\rightarrow\mathbf{\Gamma}\rightarrow\mathbf{X}$ in reciprocal space, yielding the band structure in FIG. 1 c). The following values, obtained from a manual fit to DFT results [8], were used to calculate the energy bands, $\xi^{\mathrm{Cr}}_{t_{2g}}=-5.7$ eV, $\xi^{\mathrm{Cr}}_{e_{g}}=1.56$ eV, $\xi^{\mathrm{Cl}}=-6.24$ eV, $\xi^{\mathrm{pyz}}_{\alpha_{2}}=-1.02$ eV, $\xi^{\mathrm{pyz}}_{\alpha_{4}}=-0.08$ eV, $V_{pd\sigma}=1.7$ eV, $V_{pd\pi}=-4.75$ eV, $\alpha_{2}=0.67$, $\alpha_{4}=-0.54$. Knowing that 3 of the d-levels are half-filled and the last two are empty, we employ two different on-site energies for the chromium orbitals; one value for the three $t_{2g}$ orbitals and one for the two $e_{g}$ orbitals. Further, we know from FIG. S2 that the two “orbitals” on the pyrazine sites have two distinct values. All Cl orbitals are assumed to have the same on-site energy.

The empty d-bands are located above 1.5 eV (the $e_{g}$ orbitals) whereas the three filled d-bands ($t_{2g}$) are located below the p-levels (not visible in the plot). This is evident when studying the eigenstates of each band but even more clearly seen from first principle calculations of the PDOS [8]. These d-bands are very flat and thus, the d-electrons are localized on the Cr sites. This could indicate a strong $\pi$-d hybridization. In contrary the p-bands mainly occupied by pyz electrons, located above and below the Fermi energy, are dispersive indicating that the pyz electrons are delocalized.

VII The Minimal Model and its ground states

Here, we present an effective model of monolayer CrCl₂(pyz)₂, see FIG. 2. Note that this is not identical to two unit cells of a single layer of CrCl₂(pyz)₂ since it contain fewer pyrazine sites and they are shifted with respect to the two chromium atoms. However, it essentially describes the competition between the kinetic energy and the exchange energy. In this model the two chromium atoms are viewed as localized $S=3/2$ spins. There are two electrons that can delocalize over the four pyrazine sites. The chlorine sites are neglected in this model, except for the fact that they have reduced the number of electrons on the chromium sites.

We include the exchange coupling between the chromium spins and the spins of the two pyrazine electrons. Furthermore, the pyrazine electrons can hop between the pyrazine sites (indicated with dashed lines in FIG. 2). Lastly, we include onsite Hubbard repulsion, $U>0$, for electrons located on the same site. The full Hamiltonian can thus be written as

$$\hat{\mathcal{H}}= -2J\hat{\vec{S}}_{A}\cdot\left(\hat{\vec{S}}_{1}+\hat{\vec{S}}_{2}+\hat{\vec{S}}_{3}\right)-2J\hat{\vec{S}}_{B}\cdot\left(\hat{\vec{S}}_{2}+\hat{\vec{S}}_{3}+\hat{\vec{S}}_{4}\right)$$ (S36)

$$-t\sum\limits_{\sigma=\uparrow,\downarrow}\left(\hat{c}_{1\sigma}^{\dagger}\hat{c}_{2\sigma}^{\phantom{\dagger}}+\hat{c}_{2\sigma}^{\dagger}\hat{c}_{3\sigma}^{\phantom{\dagger}}+\hat{c}_{1\sigma}^{\dagger}\hat{c}_{3\sigma}^{\phantom{\dagger}}+\hat{c}_{2\sigma}^{\dagger}\hat{c}_{4\sigma}^{\phantom{\dagger}}+\hat{c}_{3\sigma}^{\dagger}\hat{c}_{4\sigma}^{\phantom{\dagger}}+\mathrm{h.c.}\right)+U\sum\limits_{i=1}^{4}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}^{\phantom{\dagger}},$$

where $\hat{\vec{S}}_{A}$ and $\hat{\vec{S}}_{B}$ represent the spin of the two chromium atoms and $\hat{\vec{S}}_{i}$ with $i=1,2,3,4$ represent the spin of the four pyrazine electrons. The operators $\hat{c}_{i\sigma}^{\dagger}$ and $\hat{c}_{i\sigma}$ create and annihilate an electron on site $i$ with spin $\sigma$, respectively. $J$ is the exchange integral, $t$ the hopping constant and $U$ is the strength of the Hubbard repulsion.

VII.1 Basis

The spins of each chromium atom can take four $S_{z}$ values, i.e., $-3/2$, $-1/2$, $1/2$ and $3/2$. The basis for the two electrons on the pyrazine sites will be written in terms of the singlet and triplet states. For this purpose I will introduce the following notation. A singlet state will be represented by a line , e.g.,

$$|-00\rangle=\frac{1}{\sqrt{2}}\left(|\uparrow\downarrow 00\rangle-|\downarrow\uparrow 00\rangle\right),$$ (S37)

$$|\,\leavevmode\hbox to17.18pt{\vbox to9.43pt{\pgfpicture\makeatletter\hbox{\hskip 8.59045pt\lower-3.22221pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{ {{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$0$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{}{ {}{}{}}{}{{}{}}{{}}{}{}{}{}{}{}{} {}{{}{}{}}{}{}{}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-8.39046pt}{-2.26825pt}\pgfsys@curveto{-8.39046pt}{-2.26825pt}{-6.05342pt}{6.00554pt}{0.0pt}{6.00554pt}\pgfsys@curveto{6.05342pt}{6.00554pt}{8.39046pt}{-2.26825pt}{8.39046pt}{-2.26825pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{ {}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,0\rangle=\frac{1}{\sqrt{2}}\left(|\uparrow 0\downarrow 0\rangle-|\downarrow 0\uparrow 0\rangle\right), 0$$ (S38)

$$|-000\rangle=\frac{1}{\sqrt{2}}\left(|\uparrow\downarrow 000\rangle-|\downarrow\uparrow 000\rangle\right).$$ (S39)

The first state is a singlet state between an electron on site 1 and an electron on site 2, the second between two electrons on site 1 and 3, and the last is a singlet state of two electrons both occupying site 1. Similarly, we can represent a triplet state with $S_{z}=0$ as an arrow,

$$|\rightarrow 00\rangle=\frac{1}{\sqrt{2}}\left(|\uparrow\downarrow 00\rangle+|\downarrow\uparrow 00\rangle\right),$$ (S40)

$$|\,\leavevmode\hbox to18.92pt{\vbox to9.43pt{\pgfpicture\makeatletter\hbox{\hskip 8.59045pt\lower-3.22221pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{ {{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$0$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{}{ {}{}{}}{}{{}{}}{{}}{}{}{}{}{}{}{} {}{{}{}{}}{}{}{}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }{}{{ {\pgfsys@beginscope \pgfsys@setdash{}{0.0pt}\pgfsys@roundcap\pgfsys@roundjoin{} {}{}{} {}{}{} \pgfsys@moveto{-2.07999pt}{2.39998pt}\pgfsys@curveto{-1.69998pt}{0.95998pt}{-0.85318pt}{0.28pt}{0.0pt}{0.0pt}\pgfsys@curveto{-0.85318pt}{-0.28pt}{-1.69998pt}{-0.95998pt}{-2.07999pt}{-2.39998pt}\pgfsys@stroke\pgfsys@endscope}} }{}{}{{}}\pgfsys@moveto{-8.39046pt}{-2.26825pt}\pgfsys@curveto{-8.39046pt}{-2.26825pt}{-6.05342pt}{6.00554pt}{0.0pt}{6.00554pt}\pgfsys@curveto{6.05342pt}{6.00554pt}{8.39046pt}{-2.26825pt}{8.28174pt}{-1.88332pt}\pgfsys@stroke\pgfsys@invoke{ }{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{0.27184}{-0.96236}{0.96236}{0.27184}{8.33609pt}{-2.07579pt}\pgfsys@invoke{ }\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}{{}}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{ {}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,0\rangle=\frac{1}{\sqrt{2}}\left(|\uparrow 0\downarrow 0\rangle+|\downarrow 0\uparrow 0\rangle\right). 0$$ (S41)

The reason for representing the triplet state as an arrow can be seen by anticommuting operators, $\hat{c}_{1\uparrow}^{\dagger}\hat{c}_{2\downarrow}^{\dagger}+\hat{c}_{1\downarrow}^{\dagger}\hat{c}_{2\uparrow}^{\dagger}=-\left(\hat{c}_{2\uparrow}^{\dagger}\hat{c}_{1\downarrow}^{\dagger}+\hat{c}_{2\downarrow}^{\dagger}\hat{c}_{1\uparrow}^{\dagger}\right)$, such that the overall sign is meaningfully giving the direction of the arrow. For singlets the two states would be identical. All triplet states are defined such that the operator to the left creates an electron on a site labeled by a lower numerical value than the operator to the right (also applicable to $S_{z}=1$ and $S_{z}=-1$ triplets). The basis for the two electrons are thus

$$|\psi\rangle=\Big{\{} |-00\rangle,|0-0\rangle,|00-\rangle,|\,\leavevmode\hbox to17.18pt{\vbox to9.43pt{\pgfpicture\makeatletter\hbox{\hskip 8.59045pt\lower-3.22221pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{ {{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$0$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{}{ {}{}{}}{}{{}{}}{{}}{}{}{}{}{}{}{} {}{{}{}{}}{}{}{}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-8.39046pt}{-2.26825pt}\pgfsys@curveto{-8.39046pt}{-2.26825pt}{-6.05342pt}{6.00554pt}{0.0pt}{6.00554pt}\pgfsys@curveto{6.05342pt}{6.00554pt}{8.39046pt}{-2.26825pt}{8.39046pt}{-2.26825pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{ {}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,0\rangle,|0\,\leavevmode\hbox to17.18pt{\vbox to9.43pt{\pgfpicture\makeatletter\hbox{\hskip 8.59045pt\lower-3.22221pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{ {{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$0$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{}{ {}{}{}}{}{{}{}}{{}}{}{}{}{}{}{}{} {}{{}{}{}}{}{}{}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-8.39046pt}{-2.26825pt}\pgfsys@curveto{-8.39046pt}{-2.26825pt}{-6.05342pt}{6.00554pt}{0.0pt}{6.00554pt}\pgfsys@curveto{6.05342pt}{6.00554pt}{8.39046pt}{-2.26825pt}{8.39046pt}{-2.26825pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{ {}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,\rangle,|\,\leavevmode\hbox to22.18pt{\vbox to9.43pt{\pgfpicture\makeatletter\hbox{\hskip 11.09047pt\lower-3.22221pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{ {{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-5.00002pt}{-3.22221pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$00$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{}{ {}{}{}}{}{{}{}}{{}}{}{}{}{}{}{}{} {}{{}{}{}}{}{}{}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-10.89047pt}{-2.26825pt}\pgfsys@curveto{-10.89047pt}{-2.26825pt}{-7.8571pt}{6.00554pt}{0.0pt}{6.00554pt}\pgfsys@curveto{7.8571pt}{6.00554pt}{10.89047pt}{-2.26825pt}{10.89047pt}{-2.26825pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{ {}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,\rangle, 0 0 00 0 0 00$$

$$|-000\rangle,|0-00\rangle,|00-0\rangle,|000-\rangle,$$

$$|\downarrow\downarrow 00\rangle,|0\downarrow\downarrow 0\rangle,|00\downarrow\downarrow\rangle,|\downarrow 0\downarrow 0\rangle,|0\downarrow 0\downarrow\rangle,|\downarrow 00\downarrow\rangle,$$ (S42)

$$|\rightarrow 00\rangle,|0\rightarrow 0\rangle,|00\rightarrow\rangle,|\,\leavevmode\hbox to18.92pt{\vbox to9.43pt{\pgfpicture\makeatletter\hbox{\hskip 8.59045pt\lower-3.22221pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{ {{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$0$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{}{ {}{}{}}{}{{}{}}{{}}{}{}{}{}{}{}{} {}{{}{}{}}{}{}{}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }{}{}{}{}{{}}\pgfsys@moveto{-8.39046pt}{-2.26825pt}\pgfsys@curveto{-8.39046pt}{-2.26825pt}{-6.05342pt}{6.00554pt}{0.0pt}{6.00554pt}\pgfsys@curveto{6.05342pt}{6.00554pt}{8.39046pt}{-2.26825pt}{8.28174pt}{-1.88332pt}\pgfsys@stroke\pgfsys@invoke{ }{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{0.27184}{-0.96236}{0.96236}{0.27184}{8.33609pt}{-2.07579pt}\pgfsys@invoke{ }\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}{{}}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{ {}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,0\rangle,|0\,\leavevmode\hbox to18.92pt{\vbox to9.43pt{\pgfpicture\makeatletter\hbox{\hskip 8.59045pt\lower-3.22221pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{ {{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$0$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{}{ {}{}{}}{}{{}{}}{{}}{}{}{}{}{}{}{} {}{{}{}{}}{}{}{}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }{}{}{}{}{{}}\pgfsys@moveto{-8.39046pt}{-2.26825pt}\pgfsys@curveto{-8.39046pt}{-2.26825pt}{-6.05342pt}{6.00554pt}{0.0pt}{6.00554pt}\pgfsys@curveto{6.05342pt}{6.00554pt}{8.39046pt}{-2.26825pt}{8.28174pt}{-1.88332pt}\pgfsys@stroke\pgfsys@invoke{ }{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{0.27184}{-0.96236}{0.96236}{0.27184}{8.33609pt}{-2.07579pt}\pgfsys@invoke{ }\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}{{}}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{ {}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,\rangle,|\,\leavevmode\hbox to23.71pt{\vbox to9.43pt{\pgfpicture\makeatletter\hbox{\hskip 11.09047pt\lower-3.22221pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{ {{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-5.00002pt}{-3.22221pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$00$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{}{ {}{}{}}{}{{}{}}{{}}{}{}{}{}{}{}{} {}{{}{}{}}{}{}{}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }{}{}{}{}{{}}\pgfsys@moveto{-10.89047pt}{-2.26825pt}\pgfsys@curveto{-10.89047pt}{-2.26825pt}{-7.8571pt}{6.00554pt}{0.0pt}{6.00554pt}\pgfsys@curveto{7.8571pt}{6.00554pt}{10.89047pt}{-2.26825pt}{10.75279pt}{-1.8927pt}\pgfsys@stroke\pgfsys@invoke{ }{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{0.34422}{-0.93889}{0.93889}{0.34422}{10.82164pt}{-2.08047pt}\pgfsys@invoke{ }\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}{{}}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{ {}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,\rangle, 0 0 00$$

$$|\uparrow\uparrow 00\rangle,|0\uparrow\uparrow 0\rangle,|00\uparrow\uparrow\rangle,|\uparrow 0\uparrow 0\rangle,|0\uparrow 0\uparrow\rangle,|\uparrow 00\uparrow\rangle\Big{\}}.$$

The full state is then given by $\ket{ms_{A}}\otimes\ket{ms_{B}}\otimes\ket{\psi}$ and the Hilbert space is 448-dimensional.

VII.2 Hamiltonian representation

The Hamiltonian is diagonalized numerically. First, we rewrite the Hamiltonian, using $\hat{\vec{S}}_{1}\cdot\hat{\vec{S}}_{2}=\hat{S}_{1}^{z}\hat{S}_{2}^{z}+\frac{1}{2}\left(\hat{S}_{1}^{+}\hat{S}_{2}^{-}+\hat{S}_{1}^{-}\hat{S}_{2}^{+}\right)$, such that

$$\hat{\mathcal{H}}= -2J\hat{S}_{A}^{z}\left(\hat{S}_{1}^{z}+\hat{S}_{2}^{z}+\hat{S}_{3}^{z}\right)-J\hat{S}_{A}^{+}\left(\hat{S}_{1}^{-}+\hat{S}_{2}^{-}+\hat{S}_{3}^{-}\right)-J\hat{S}_{A}^{-}\left(\hat{S}_{1}^{+}+\hat{S}_{2}^{+}+\hat{S}_{3}^{+}\right)$$ (S43)

$$-2J\hat{S}_{B}^{z}\left(\hat{S}_{2}^{z}+\hat{S}_{3}^{z}+\hat{S}_{4}^{z}\right)-J\hat{S}_{B}^{+}\left(\hat{S}_{2}^{-}+\hat{S}_{3}^{-}+\hat{S}_{4}^{-}\right)-J\hat{S}_{B}^{-}\left(\hat{S}_{2}^{+}+\hat{S}_{3}^{+}+\hat{S}_{4}^{+}\right)$$

$$-t\sum\limits_{\sigma=\uparrow,\downarrow}\left(\hat{c}_{1\sigma}^{\dagger}\hat{c}_{2\sigma}^{\phantom{\dagger}}+\hat{c}_{2\sigma}^{\dagger}\hat{c}_{3\sigma}^{\phantom{\dagger}}+\hat{c}_{1\sigma}^{\dagger}\hat{c}_{3\sigma}^{\phantom{\dagger}}+\hat{c}_{2\sigma}^{\dagger}\hat{c}_{4\sigma}^{\phantom{\dagger}}+\hat{c}_{3\sigma}^{\dagger}\hat{c}_{4\sigma}^{\phantom{\dagger}}+\mathrm{h.c.}\right)$$

$$+U\sum\limits_{i=1}^{4}\sum\limits_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}\hat{c}_{i\sigma^{\prime}}^{\dagger}\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}^{\phantom{\dagger}}\hat{c}_{i\sigma^{\prime}}^{\phantom{\dagger}}.$$

In the following we will explore how the full Hamiltonian can be written in the basis described above by investigating some specific examples. The less enthusiastic reader can potentially skip this part and jump to section VII.3.

We focus on the first three exchange terms including $\hat{S}_{A}^{z}$. The chromium part of this is trivial in the basis $\{\ket{-3/2},\ket{-1/2},\ket{1/2},\ket{3/2}\}$:

$$M_{\mathrm{Cr}}^{z}=\begin{bmatrix}-\frac{3}{2}&0&0&0\\ 0&-\frac{1}{2}&0&0\\ 0&0&\frac{1}{2}&0\\ 0&0&0&\frac{3}{2}\end{bmatrix}.$$ (S44)

Now, let us see how this affects the basis states for the electrons on the pyrazine sites. Applying this operator on basis states where neither of the two electrons are located on site 4 results in the expected $S_{z}=0$ for a singlet state, e.g.,

$$\left(\hat{S}_{1}^{z}+\hat{S}_{2}^{z}+\hat{S}_{3}^{z}\right)\ket{-00}=\frac{1}{\sqrt{2}}\left(\frac{1}{2}\ket{\uparrow\downarrow 00}-\frac{1}{2}\ket{\uparrow\downarrow 00}-\frac{1}{2}\ket{\downarrow\uparrow 00}+\frac{1}{2}\ket{\downarrow\uparrow 00}\right)=0.$$ (S45)

The same applies to triplet states with $S_{z}=0$ with no electrons occupying site 4 and for singlet states where both electrons occupy the same site. On the other hand, when one of the electrons is occupying site 4, singlet and triplet states couple. For instance:

$$\left(\hat{S}_{1}^{z}+\hat{S}_{2}^{z}+\hat{S}_{3}^{z}\right)\ket{00-}=\frac{1}{\sqrt{2}}\left[\frac{1}{2}\ket{00\uparrow\downarrow}-\left(-\frac{1}{2}\ket{00\downarrow\uparrow}\right)\right]=\frac{1}{2}\ket{00\rightarrow}.$$ (S46)

The triplet states with $S_{z}=1$ or $S_{z}=-1$ are eigenstates of this operator; the eigenvalue depends on whether one of the electrons is occupying site 4 or not,

$$\left(\hat{S}_{1}^{z}+\hat{S}_{2}^{z}+\hat{S}_{3}^{z}\right)\ket{0\downarrow\downarrow 0}=-\ket{0\downarrow\downarrow 0},$$ (S47)

$$\left(\hat{S}_{1}^{z}+\hat{S}_{2}^{z}+\hat{S}_{3}^{z}\right)\ket{00\downarrow\downarrow}=-\frac{1}{2}\ket{00\downarrow\downarrow},$$ (S48)

$$\left(\hat{S}_{1}^{z}+\hat{S}_{2}^{z}+\hat{S}_{3}^{z}\right)\ket{0\uparrow\uparrow 0}=\ket{0\uparrow\uparrow 0},$$ (S49)

$$\left(\hat{S}_{1}^{z}+\hat{S}_{2}^{z}+\hat{S}_{3}^{z}\right)\ket{00\uparrow\uparrow}=\frac{1}{2}\ket{00\uparrow\uparrow}.$$ (S50)

This procedure is similar for the $S_{B}^{z}$ part of the Hamiltonian but in this case singlet and triplet states are coupled when one of the electrons is occupying site 1.

Now, let us look at the other part of the exchange interaction, i.e., the terms
$-J\hat{S}_{A}^{+}\left(\hat{S}_{1}^{-}+\hat{S}_{2}^{-}+\hat{S}_{3}^{-}\right)-J\hat{S}_{A}^{-}\left(\hat{S}_{1}^{+}+\hat{S}_{2}^{+}+\hat{S}_{3}^{+}\right)$. The raising and lowering operators for the chromium spins can be written as

$$M_{\mathrm{Cr}}^{+}=\begin{bmatrix}0&0&0&0\\ \sqrt{3}&0&0&0\\ 0&2&0&0\\ 0&0&\sqrt{3}&0\end{bmatrix}\quad\text{and}\quad M_{\mathrm{Cr}}^{-}=\begin{bmatrix}0&\sqrt{3}&0&0\\ 0&0&2&0\\ 0&0&0&\sqrt{3}\\ 0&0&0&0\end{bmatrix}.$$ (S51)

For the electrons on the pyrazine site it is not as simple. As before, the outcome depends on whether one of the electrons is occupying site 4 or not. Focusing on the electron part including raising operators for the pyrazine electrons, we look at a few cases to see how the basis states are affected by this operator,

$$\left(\hat{S}_{1}^{+}+\hat{S}_{2}^{+}+\hat{S}_{3}^{+}\right)\ket{-00}=0,$$ (S52)

$$\left(\hat{S}_{1}^{+}+\hat{S}_{2}^{+}+\hat{S}_{3}^{+}\right)\ket{00-}=-\frac{1}{\sqrt{2}}\ket{00\uparrow\uparrow},$$ (S53)

$$\left(\hat{S}_{1}^{+}+\hat{S}_{2}^{+}+\hat{S}_{3}^{+}\right)\ket{\rightarrow 00}=\sqrt{2}\ket{\uparrow\uparrow 00},$$ (S54)

$$\left(\hat{S}_{1}^{+}+\hat{S}_{2}^{+}+\hat{S}_{3}^{+}\right)\ket{00\rightarrow}=\frac{1}{\sqrt{2}}\ket{00\uparrow\uparrow},$$ (S55)

$$\left(\hat{S}_{1}^{+}+\hat{S}_{2}^{+}+\hat{S}_{3}^{+}\right)\ket{00\downarrow\downarrow}=\ket{00\uparrow\downarrow},$$

$$\quad\quad=\frac{1}{2}\left(\ket{00\uparrow\downarrow}+\ket{00\downarrow\uparrow}\right)+\frac{1}{2}\left(\ket{00\uparrow\downarrow}-\ket{00\downarrow\uparrow}\right)$$

$$\quad\quad=\frac{1}{\sqrt{2}}\ket{00\rightarrow}+\frac{1}{\sqrt{2}}\ket{00-}.$$ (S56)

The final example is quite interesting since this results in a superposition between a triplet state and a singlet state. The same method can be used to obtain the resulting states when applying the operator $\left(\hat{S}_{1}^{-}+\hat{S}_{2}^{-}+\hat{S}_{3}^{-}\right)$ and similarly for the B-part of the exchange coupling.
Now, we look at the tight binding part of the Hamiltonian. Since the spin is unchanged by this part, the singlet and the triplet states will not couple, and thus we can look at these parts of the basis states separately. First, let us focus on the 10 singlet states. Here, one example is sufficient to illustrate what happens to all of them:

$$\hat{\mathcal{H}}_{t}\ket{-00}=-t\left(\sqrt{2}\ket{-000}+\sqrt{2}\ket{0-00}+|\,\leavevmode\hbox to17.18pt{\vbox to9.43pt{\pgfpicture\makeatletter\hbox{\hskip 8.59045pt\lower-3.22221pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{ {{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-3.22221pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$0$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{}{ {}{}{}}{}{{}{}}{{}}{}{}{}{}{}{}{} {}{{}{}{}}{}{}{}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-8.39046pt}{-2.26825pt}\pgfsys@curveto{-8.39046pt}{-2.26825pt}{-6.05342pt}{6.00554pt}{0.0pt}{6.00554pt}\pgfsys@curveto{6.05342pt}{6.00554pt}{8.39046pt}{-2.26825pt}{8.39046pt}{-2.26825pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{ {}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,0\rangle+\ket{0-0}+|\,\leavevmode\hbox to22.18pt{\vbox to9.43pt{\pgfpicture\makeatletter\hbox{\hskip 11.09047pt\lower-3.22221pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{ {{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-5.00002pt}{-3.22221pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$00$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{}{ {}{}{}}{}{{}{}}{{}}{}{}{}{}{}{}{} {}{{}{}{}}{}{}{}{}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }{}\pgfsys@moveto{-10.89047pt}{-2.26825pt}\pgfsys@curveto{-10.89047pt}{-2.26825pt}{-7.8571pt}{6.00554pt}{0.0pt}{6.00554pt}\pgfsys@curveto{7.8571pt}{6.00554pt}{10.89047pt}{-2.26825pt}{10.89047pt}{-2.26825pt}\pgfsys@stroke\pgfsys@invoke{ } \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{ {}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,\rangle\right). 0 00 0 00$$ (S57)

What is important to note here is that a factor of $\sqrt{2}$ arises when a singly occupied state become a doubly occupied site, since

$$\sum\limits_{\sigma}\hat{c}_{1\sigma}^{\dagger}\hat{c}_{2\sigma}\ket{-00} =\sum\limits_{\sigma}\hat{c}_{1\sigma}^{\dagger}\hat{c}_{2\sigma}\frac{1}{\sqrt{2}}\left(\hat{c}_{1\uparrow}^{\dagger}\hat{c}_{2\downarrow}^{\dagger}-\hat{c}_{1\downarrow}^{\dagger}\hat{c}_{2\uparrow}^{\dagger}\right)\ket{\mathrm{vac}}$$

$$=\frac{1}{\sqrt{2}}\left(\hat{c}_{1\uparrow}^{\dagger}\hat{c}_{1\downarrow}^{\dagger}-\hat{c}_{1\downarrow}^{\dagger}\hat{c}_{1\uparrow}^{\dagger}\right)\ket{\mathrm{vac}}=\sqrt{2}\ket{-000}.$$ (S58)

Here $\ket{\mathrm{vac}}$ represents the vacuum state. Due to the antisymmetry of a singlet state we do not need to think about signs. Again, one example is sufficient:

$$\hat{\mathcal{H}}_{t}\ket{\uparrow\uparrow 00} =\sum\limits_{\sigma}\left(\hat{c}_{3\sigma}^{\dagger}\hat{c}_{1\sigma}+\hat{c}_{3\sigma}^{\dagger}\hat{c}_{2\sigma}+\hat{c}_{4\sigma}^{\dagger}\hat{c}_{2\sigma}\right)\ket{\uparrow\uparrow 00}$$

$$=\left(\hat{c}_{3\uparrow}^{\dagger}\hat{c}_{1\uparrow}\hat{c}_{1\uparrow}^{\dagger}\hat{c}_{2\uparrow}^{\dagger}+\hat{c}_{3\uparrow}^{\dagger}\hat{c}_{2\uparrow}\hat{c}_{1\uparrow}^{\dagger}\hat{c}_{2\uparrow}^{\dagger}+\hat{c}_{4\uparrow}^{\dagger}\hat{c}_{2\uparrow}\hat{c}_{1\uparrow}^{\dagger}\hat{c}_{2\uparrow}^{\dagger}\right)\ket{\mathrm{vac}}$$

$$=\left(-\hat{c}_{2\uparrow}^{\dagger}\hat{c}_{3\uparrow}^{\dagger}+\hat{c}_{1\uparrow}^{\dagger}\hat{c}_{3\uparrow}^{\dagger}+\hat{c}_{1\uparrow}^{\dagger}\hat{c}_{4\uparrow}^{\dagger}\right)\ket{\mathrm{vac}}$$

$$=-\ket{0\uparrow\uparrow 0}+\ket{\uparrow 0\uparrow 0}+\ket{\uparrow 00\uparrow}.$$ (S59)

The last term in the Hamiltonian (onsite repulsion) is trivial to include for the four singlet states. Using this basis the $448\crossproduct 448$ Hamiltonian can be diagonalized numerically yielding the results described in the main text. The main conclusion when diagonalizing this Hamiltonian is that the two degenerate ground state are dominated by the two basis states viewed in Eq. (3) and Fig. 4.

VII.3 Two Level Model

Focusing on the dominant basis state in the ground state obtained in the model described above, see Eq. (3) and Fig. 4, i.e., going to the limit $J\gg t$, we can make a simple model to understand how the ground state energy changes when increasing the hopping parameter $t$. For $t=0$ the two electrons will localize on sites 2 and 3 with spins antiparallel compared to the chromium spins. When $t$ increases the two electrons can then hop between two different sites and still avoid each other, which they can do in two different configurations. One electron can hop between site 1 and 2 while the other hop between 3 and 4, see FIG. 4 a), or one electron can hop between site 1 and 3 while the other can hop between 2 and 4, see FIG. 4 b).

For an electron hopping between two sites, e.g., 1 and 2, the two levels do not have the same energy, since the electron can couple to both Cr spins at site 2 but only couple to one Cr spin at site 1. The gap between these two levels is $\Delta=\frac{3J}{4}$ since one level has the energy $\left(\frac{3}{2}+\frac{3}{2}\right)J=3J$ and the other has energy $\frac{3}{2}J$. Thus, we employ the simple Hamiltonian

$$\hat{\mathcal{H}}=\begin{bmatrix}\frac{3J}{4}&t\\ t&\frac{-3J}{4}\end{bmatrix},$$ (S60)

with the eigenvalues

$$E=\pm\sqrt{\left(\frac{3J}{4}\right)^{2}+t^{2}}.$$ (S61)

The negative solution is the ground state. In FIG. S3 we show this function alongside the exact eigenvalues as a function of $t$ with $J=-0.2529$ eV and $U=0.2529$ eV. The value chosen for $J$ is motivated by Ref. [8], with $U$ chosen to be of the same magnitude as $J$. Note that Eq. (S61) has been shifted to fit with the exact eigenvalues in the figure. The exact ground state energies follow this simple model quite closely in the entire $t$ regime investigated. The exact solution is a bit lower in energy, which can be explained by the fact that in this model the system has two possible configurations with the same low energy, see FIG. 4. This only lowers the energy slightly since it has to be coordinated such that the electrons can still avoid each other. These features are robust towards changes in $U$ (even turning off repulsion completely) since electrons occupying the same site is not favorable in this model.

FIG. S3: The five lowest, degenerate eigenenergies plotted as a function of t with $J=-0.2529$ eV and $U=0.2529$ eV, compared with the function in Eq. (S61), shifted to match at t = 0.

When $t=0$ the ground state energy becomes five-fold degenerate. Setting $t\neq 0$ this energy splits in three of which two are two-fold degenerate. The five-fold degeneracy in the limit where $J\gg t$ indicates that the system behaves as a spin-2 particle. This makes sense since

$$\hat{S}_{\mathrm{tot}}^{z}\ket{3/2}\otimes\ket{3/2}\otimes\ket{0\downarrow\downarrow} =\left(\hat{S}_{A}^{z}+\hat{S}_{B}^{z}+\hat{S}_{2}^{z}+\hat{S}_{3}^{z}\right)\ket{3/2}\otimes\ket{3/2}\otimes\ket{0\downarrow\downarrow}$$ (S62)

$$=\left(\frac{3}{2}+\frac{3}{2}-\frac{1}{2}-\frac{1}{2}\right)\ket{3/2}\otimes\ket{3/2}\otimes\ket{0\downarrow\downarrow}$$

$$=2\ket{3/2}\otimes\ket{3/2}\otimes\ket{0\downarrow\downarrow}.$$

All these five lowest states are plotted in FIG. S3.

VII.4 Zeeman splitting

To see how these states split up in a magnetic field we add the term

$$\hat{\mathcal{H}}_{B}=g\mu_{B}\mathbf{B}\cdot\hat{\mathbf{S}}_{\mathrm{tot}}=g\mu_{B}B^{z}\hat{S}_{\mathrm{tot}}^{z},$$ (S63)

where we choose the magnetic field to point in the $z$-direction, i.e. out of the plane in FIG. 2. As seen in FIG. S4 a) the five-fold degeneracy of the ground state for $t=0$ is lifted when a magnetic field is applied. Furthermore, we see that for $t\neq 0$ the energy splits into three, two of which are two-fold degenerate. The two doubly degenerate states act as a doublet and the non-degenerate state behaves as a singlet.

FIG. S4: The five lowest eigenenergies of Eq. (S36) + Eq. (S63) depicted as a function of $B_{z}$ with $J=-0.2529$ eV$=-U$ for $t=0$ eV (a) and $t=0.2529$ eV (b). The degeneracy is lifted by Zeeman splitting for $B\neq 0$.

Adding a magnetic field also allows us to calculate the expectation value of the total spin using

$$\expectationvalue{S_{\mathrm{tot}}}=\frac{\partial E}{\partial B}.$$ (S64)

Both for $t=0$ and $t\neq 0$ this is found numerically to be 2.

VIII Estimate of the effective chromium coupling from mean-field theory

From mean-field theory we can derive an expression of the exchange coupling constant between neighboring chromium spins $\mathbf{S}_{i}$. We depart from

$$\hat{\mathcal{H}}=-J_{\text{Cr-Cr}}\sum\limits_{\expectationvalue{i,j}}\mathbf{S}_{i}\cdot\mathbf{S}_{j}.$$ (S65)

Using mean-field theory and neglecting the constant shift, the interaction of the $i$th spin is

$$\hat{\mathcal{H}}_{i}=-J_{\text{Cr-Cr}}z\expectationvalue{\mathbf{S}}\cdot\mathbf{S}_{i},$$ (S66)

where $z$ denotes the number of neighboring spins (the coordination number), i.e., $z=4$ for a square lattice. Defining

$$\mathbf{H}\equiv-\frac{zJ_{\text{Cr-Cr}}}{g\mu_{B}}\expectationvalue{\mathbf{S}},$$ (S67)

the Hamiltonian for the $i$th spin can be written as

$$\hat{\mathcal{H}}_{i}=g\mu_{B}\mathbf{S}_{i}\cdot\mathbf{H}.$$ (S68)

Thus, $\mathbf{H}$ can be viewed as an exchange field affecting the $i$th spin, which arises from the polarization of the spins on the neighboring chromium sites. The magnetization from this exchange field is then given by

$$\mathbf{M}=-g\mu_{B}N\expectationvalue{\mathbf{S}},$$ (S69)

where $N$ is the number of magnetic chromium atoms in the sample. From this an expression for $\expectationvalue{\mathbf{S}}$ can be obtained, such that

$$\mathbf{H}=A\mathbf{M},\quad\text{with}\quad A=\frac{zJ_{\text{Cr-Cr}}}{g^{2}\mu_{B}^{2}N}.$$ (S70)

Including an external magnetic field $\mathbf{B}$ the total magnetization magnitude can be written as

$$M=\frac{C}{T}\left(B+AM\right),$$ (S71)

where $C$ is the Curie constant, which for $L=0$, where $L$ is the total orbital angular momentum, is given by

$$C=\frac{\mu_{0}\mu_{B}^{2}}{3k_{B}}Ng^{2}S(S+1).$$ (S72)

The magnetization magnitude can now be isolated:

$$M=\frac{C}{T-AC}B.$$ (S73)

Using the Curie Weiss law, we can conclude that $AC=T_{C}$, where $T_{C}$ is the Curie temperature, from which we obtain the expression

$$J_{\text{Cr-Cr}}=\frac{3k_{B}T_{C}}{\mu_{0}zS(S+1)}.$$ (S74)

IX Exchange Interaction between Cr and pyz

Here we study the interaction between chromium spins and neighboring pyrazine spins to provide an explanation as to why CrCl₂(pyz)₂ is ferrimagnetic. An expression for the change in energy, that arise from the interaction, is obtained through second order perturbation theory [34],

$$\Delta E^{(2)}=-\frac{|\langle\psi|\hat{\mathcal{H}}_{t}|\phi\rangle|^{2}}{-\Delta+U}-\frac{|\langle\psi|\hat{\mathcal{H}}_{t}|\phi\rangle|^{2}}{\Delta}=-|\langle\psi|\hat{\mathcal{H}}_{t}|\phi\rangle|^{2}\left(\frac{1}{(1-\Delta/U)\Delta}\right)\equiv J.$$ (S75)

$|\psi\rangle$ represents the highest occupied energy level on the chromium site in the environment. This is expressed as a superposition of the five d-orbitals and is obtained from the band structure studied previously. $|\phi\rangle$ can be written as

$$|\phi\rangle=\frac{1}{\sqrt{2}}\left(|\phi_{1}\rangle+|\phi_{2}\rangle\right),$$ (S76)

where $|\phi_{1}\rangle$ and $|\phi_{2}\rangle$ are the eigenstates of the relevant energy level on the two pyrazine sites in the unit cell, also obtained from the band structure. We assume that $J/t\ll 1$. $\Delta$ denotes the difference in energy of the two relevant energy levels and $U$ is the strength of the Hubbard repulsion for two electrons occupying the same chromium site, see FIG. S5. The Hubbard repulsion for the pyrazine sites is neglected due to their much larger spatial extent compared to chromium.

FIG. S5: Diagram illustrating the two ways in which an electron on the chromium site of one unit cell (five left energy levels in each illustration) can exchange location with an electron on a pyrazine site (6 right energy levels).

In the exchange interaction two electrons exchange location with one intermediate step. This can be done in two ways, both are illustrated in FIG. S5. The exchange interaction constant $J$ is a sum over these two possibilities, which results in Eq. (S75). For the coupling between the Cr sites and the pyrazine sites to be antiferromagnetic, it is required that $\Delta>U$.

IX.1 RKKY-Interaction

In this section, we derive an expression for the effective chromium coupling, as determined by means of second order perturbation theory of the RKKY interaction mediated through intermediate pyrazine sites. The starting point of this is the exchange interaction between the spin of the Cr ions, $\mathbf{S}$, and the spin of the valence electron, $\mathbf{\hat{s}}$, on the pyrazine sites, described above, with the interaction Hamiltonian

$$\hat{\mathcal{H}}^{\prime}=-J\sum\limits_{i,\delta}\mathbf{S}_{i}\cdot\mathbf{\hat{s}}_{i+\delta},$$ (S77)

for unit cell $i$. $\delta$ sums over the vectors pointing to the four neighboring pyrazine sites. We treat this interaction as a perturbation to the ground state configuration. The second order correction to the energy is given by

$$E_{\mathbf{k}}^{(2)}=\sum\limits_{\mathbf{k}^{\prime}\neq\mathbf{k}}\frac{|\langle\psi_{\mathbf{k}^{\prime}}^{0}|\hat{\mathcal{H}}^{\prime}|\psi_{\mathbf{k}}^{0}\rangle|^{2}}{\xi_{\mathbf{k}}^{0}-\xi_{\mathbf{k}^{\prime}}^{0}}=\sum\limits_{\mathbf{k}^{\prime}\neq\mathbf{k}}\frac{\langle\psi_{\mathbf{k}}^{0}|\hat{\mathcal{H}}^{\prime}|\psi_{\mathbf{k}^{\prime}}^{0}\rangle\langle\psi_{\mathbf{k}^{\prime}}^{0}|\hat{\mathcal{H}}^{\prime}|\psi_{\mathbf{k}}^{0}\rangle}{\xi_{\mathbf{k}}^{0}-\xi_{\mathbf{k}^{\prime}}^{0}},$$ (S78)

where $\ket{\psi_{\mathbf{k}}^{0}}$ is the ground state of the delocalized electrons, i.e., Bloch waves, and where $\xi_{\boldsymbol{k}}^{0}$ is the normal-state, unperturbed dispersion relation. The chromium spins in $\hat{\mathcal{H}}^{\prime}$ can be taken outside of the averages since these are not considered as a part of the perturbation. Now, we are left with a product of two correlation functions. Recall that we can write the spin of the valence electrons as

$$\sum\limits_{i}\mathbf{\hat{s}}_{i+\delta}=\frac{1}{2}\hat{c}_{i+\delta,\sigma}^{\dagger}\mathbf{\tau}_{\sigma\sigma^{\prime}}\hat{c}_{i+\delta\sigma^{\prime}}=\frac{1}{2N}\sum\limits_{\mathbf{kk}^{\prime}}\hat{c}_{\mathbf{k}\sigma}\mathbf{\tau}_{\sigma\sigma^{\prime}}\hat{c}_{\mathbf{k}\sigma^{\prime}}.$$ (S79)

In the textbook version of the RKKY model the valence electrons are considered as a free electron gas. In this case we consider the electrons to be delocalized but anchored to the four pyz sites. For the spin operators in each correlation function we must sum over all bonds to the relevant chromium spin. The dominant terms here are thus the terms where the neighboring chromium atoms couple through an electron located at a pyz site in between. Thus, we find

$$E_{n}^{(2)}=J_{\text{Cr-Cr}}=\frac{J^{2}}{2N^{2}}\mathbf{S}_{1}\cdot\mathbf{S}_{2}\,\chi,$$ (S80)

where $\chi$ is the Lindhard function given by

$$\chi=\,\,\leavevmode\hbox to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } \feynman \vertex(a); \vertex[right=of a] (b); \vertex[right=of b] (c); \vertex[right=of c] (d); \diagram* { (a) -- [photon] (b) -- [fermion, half right, looseness=1.5] (c) -- [fermion, half right, looseness=1.5] (b), (c) -- [photon] (d), }; \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{ {}{}{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,\,=\sum\limits_{\mathbf{kk}^{\prime}}\frac{n_{F}(\xi_{\mathbf{k}})-n_{F}(\xi_{\mathbf{k}^{\prime}})}{\xi_{\mathbf{k}}-\xi_{\mathbf{k}^{\prime}}}.$$ (S81)

The factor of $1/N^{2}$ stems from the Fourier transformation and the factor of two stems from the fact that

$$\langle\psi_{\mathbf{k}\sigma}|\mathbf{\tau}^{\alpha}\mathbf{\tau}^{\beta}|\psi_{\mathbf{k}\sigma}\rangle=\mathrm{Tr}[\mathbf{\tau}^{\alpha}\mathbf{\tau}^{\beta}]=2\delta_{\alpha\beta}.$$ (S82)

At $T=0$ the Fermi functions in Eq. (S81) are step functions. The Lindhard function can then be written as

$$\chi=2\sum\limits_{\mathbf{k}>\mathbf{k}_{F}}\sum\limits_{\mathbf{k}^{\prime}<\mathbf{k}_{F}}\frac{1}{\xi_{\mathbf{k}}-\xi_{\mathbf{k}^{\prime}}}=2\int\limits_{\xi_{F}}d\xi\int\limits^{\xi_{F}}d\xi^{\prime}\frac{D(\xi)D(\xi^{\prime})}{\xi-\xi^{\prime}},$$ (S83)

where $D(\xi)$ is the density of states. The density of states for a square lattice with spectrum $\xi=2t(\cos k_{x}+\cos k_{y})$ is given by

$$D(\xi)=\frac{N^{2}}{2(\pi t)^{2}}K\left[1-\left(\frac{\xi}{4t}\right)^{2}\right].$$ (S84)

Here $K$ is the complete elliptic integral of the first kind. In the subspace of the pyrazine sites there is $1/2$ electron in each unit cell and thus the band has filling $1/4$ due to spin.

The density of states is normalized such that $\int\limits_{-4t}^{4t}d\xi D(\xi)=1.$ Thus, the Fermi energy can be found by requiring that $\int\limits_{-4t}^{\xi_{F}}d\xi D(\xi)=\frac{1}{4}.$ The Fermi energy is obtained numerically and the Lindhard function is found to have a numerical value of $0.15$ for $t=1$, such that

$$J_{\text{Cr-Cr}}=0.15\frac{J^{2}}{2t}.$$ (S85)

The value $\xi_{F}/t=0.15$ obtained within the above approximations demonstrates that the Fermi surface is well-approximated by a small circle.