Dispersionless orbital excitations in (Li,Fe)OHFeSe superconductors

The superconducting critical temperature Tc of intercalated iron-selenide superconductor (Li,Fe)OHFeSe (FeSe11111) can be increased to 42 from 8 K of bulk FeSe. It shows remarkably similar electronic properties as the high-Tc monolayer FeSe and provides a bulk counterpart to investigate the origin of enhanced superconductivity. Unraveling the nature of excitations is crucial for understanding the pairing mechanism in high-Tc iron selenides. Here we use resonant inelastic x-ray scattering (RIXS) to investigate the excitations in FeSe11111. Our high-quality data exhibit several Raman-like excitations, which are dispersionless and isotropic in momentum transfer in both superconducting 28 K and 42 K samples. Using atomic multiplet calculations, we assign the low-energy ~0.3 and 0.7 eV Raman peaks as local eg − eg and eg − t2g orbital excitations. The intensity of these two features decreases with increasing temperature, suggesting a dominating contribution of the orbital fluctuations. Our results highlight the importance of the orbital degree of freedom for high-Tc iron selenides.


INTRODUCTION
After more than a decade of the discovery of Fe-based superconductors 1 , the origin of superconductivity is still under debate 2 . Spin fluctuations between the hole pockets at the Brillouin zone center and the electron pockets at the zone edges have been considered as a candidate medium of the electron pairing, which typically leads to sign-reversed s ± pairing 3,4 . This mechanism has been supported by the spin resonance mode observed by neutron scattering measurements below the superconducting gap 5 . Apart from the important role of magnetism 6 , other factors, and mechanisms have also been proposed in Febased superconductors [7][8][9][10][11][12] . For example, orbital fluctuations have been proposed to enhance spin fluctuation-mediated pairing 8 or directly account for the high T c in iron pnictides 7,9 . Angle-resolved photoemission spectroscopy (ARPES) measurements on iron chalcogenides revealed the disappearance of d xy -orbital spectral weight with increasing temperature, implying universal orbitalselective correlation effects 12 .
Among the Fe-based superconductors, FeSe has the simplest crystal structure and provides an ideal platform for the study of the pairing mechanism 13 . While the T c of bulk FeSe is only~8 K, it can be enhanced by inter-layer interactions: both the intercalation between FeSe layers and SrTiO 3 (STO) substrates for the monolayer FeSe can enhance the T c by a factor of over four [14][15][16][17] . Two possible pairing mechanisms were proposed for the enhanced superconductivity in monolayer FeSe-interfacial charge-transfer and electron-phonon coupling (EPC) 14,[18][19][20][21][22] . However, these scenarios call for verification in an effective monolayer material without the substrate and related phonons.
LiOH-intercalated FeSe (FeSe11111) is a pure bulk singlecrystalline superconductor with a T c of over 40 K 16,17 , which satisfies these conditions. Due to the intercalation, the distance between FeSe layers becomes as large as 9 Å in FeSe11111, leading to a highly two-dimensional structure. ARPES experiments also reflect the similarity between FeSe11111 and monolayer FeSe/STO in terms of Fermi surface topology, band structure, and the gap symmetry 23 . As such, FeSe11111 resembles the monolayer FeSe, avoiding the interface to the STO substrate, and therefore provides a unique opportunity to elucidate the origin of the increased T c evolving from bulk FeSe to monolayer FeSe/STO. Studies of FeSe11111 recently lead to an even more exciting observation of a Majorana zero mode implying its nontrivial topology 24 . Therefore, it is appealing to investigate the excitations in FeSe11111 as a crucial step toward understanding the pairing mechanism in iron chalcogenides.
To identify the nature of excitations in FeSe11111 and their connection to superconductivity, we employ the resonant inelastic x-ray scattering (RIXS) technique, which is sensitive to multiple excitations 25 and has been widely used to study Fe-based superconductors [26][27][28][29][30][31] . We uncover the excitations in FeSe11111 films grown on a LaAlO 3 substrate with T c ≃ 42 K (SC42K) and T c ≃ 28 K (SC28K) ("METHODS"). The higher T c accompanies a larger c lattice parameter and a lower Fe vacancy concentration in the FeSe-layer [32][33][34] . We identified six excitations in both samples, including four Raman-like excitations and two fluorescence-like excitations. The four Raman features are located at~0.1,~0.3, 0.7, and~2.5 eV. Through their momentum and temperature dependence, we attribute the~0.3 and~0.7 eV features to orbital excitations and the~2.5 eV feature to multiple spin excitations. These orbital excitations are sensitive to the distorted tetrahedral crystal field and, accordingly, the lattice geometry. Therefore, our study suggests the lattice-induced orbital fluctuations as an ingredient to enhance the T c of iron chalcogenides. Figure 1a shows the crystal structure of the FeSe-based superconductor FeSe11111 in the one-Fe unit notation. It consists of anti-PbO-type FeSe layers sandwiched by anti-PbO-type (Li/Fe)OH layers along c-direction and belongs to the space group P4/nmm (No. 129) 16 . The schematics of the RIXS experimental geometry is shown in Fig. 1b. The accessible reciprocal space at Fe L 3 -edge is displayed in Fig. 1c. We state coordinates in units of 2π/a, i.e., 1 corresponds to one reciprocal lattice unit (r.l.u.). Figure 1d shows a representative RIXS spectrum of FeSe11111 collected at the Fe L 3 -edge resonant energy (709 eV) using πpolarized incident x-rays. Intriguingly, the RIXS spectrum of FeSe11111 highly resembles that of monolayer FeSe, including the weak fluorescence background and broad features at around 0.5 and 3 eV 31 . This indicates that these features arise from the FeSe-layer instead of the Li-Fe layer. Moreover, we obtained much better statistics of RIXS spectra in FeSe11111 than those in monolayer FeSe due to the larger sample volume. It is noteworthy that, due to the metallic nature of Fe-based superconductors, the fluorescence background is substantially strong in the RIXS spectra of bulk FeSe 28,29,31 and iron pnictides 26,27 . Here, the weak fluorescence background is probably due to the absence of hole Fermi surface pockets at the Γ point, which closes the particle-hole scattering channel and allows us to clearly identify the excitations in FeSe11111. To quantitatively analyze the data, we fitted the RIXS spectra with Gaussian profiles and identified six excitations, as shown in Fig. 1d Since the incident-energy detuning RIXS measurements can distinguish between Raman-like excitations and fluorescence-like excitations 35,36 , we measured a series of RIXS spectra across Fe L 3edge in steps of 0.25 eV for SC42K (Fig. 1e). For Raman-like features, the emitted photon energy shifts by the same energy as the incident photon energy. For fluorescence-like features, the energy loss changes by the same energy as the incident photon energy 25 . The features at around 1.1 and 4.8 eV at the Fe L 3 -edge resonant energy track the incident energy [see Fig. 1e, f], which are typical fluorescence behaviors. They correspond to the Fe 3d to 2p fluorescence decay reported also in iron pnictides and Fe 3d to Se 4p fluorescence decay, respectively 26,37,38 . The Raman-like features include an evident shoulder at~0.1 eV and three excitations at around 0.3, 0.7, and 2.5 eV, which are captured by our atomic multiplet calculations as discussed below. We obtained similar detuning results in SC28K (see Supplementary Fig. 4).

Momentum and T c dependence
To understand the Raman-like excitations, we measured the detailed momentum-resolved RIXS maps along two highsymmetry directions,  Fig. 7), which is much broader than our energy resolution (~0.08 eV). They resemble those around 0.5 eV in monolayer FeSe, which were attributed to spin excitations 31 . It should be noticed that this comparison is not limited by energy resolution since the energies of these excitations (200-700 meV) are more than 2. times the resolution of our RIXS experiment (~80 meV). These excitations have similar energies in both SC28K and SC42K, irrespective of varying T c (Fig. 2g). We also measured the momentum-dependent RIXS spectra at 95 K and 200 K ( Fig. 2c-f), which showed the persistence of these Raman-like excitations to high temperatures. While their energies do not change with temperature ( Fig. 2g-i), the intensities descend with increasing temperature ( Fig. 2a-f). We will discuss the temperature dependence in more detail below. The dispersionless behavior of the Raman-like features suggests that the spectral peaks should be of local nature. Different from d 9 copper oxide superconductors where the low-energy states lie in a single 3d orbital, the hybridization and interactions among five partially filled 3d orbitals in Fe allow rich atomic excitations, which can be reflected in the RIXS spectra. Using atomic multiplet theory with a distorted tetrahedral crystal field environment (CFE), we were able to reproduce theoretical spectral peaks at the experimentally observed excitation energies [see comparisons in Fig. 3a, b]. By projecting the excited-state wavefunctions to the maximal-overlapped Fock states, the nature of each excitation can be addressed approximately in a single-particle picture [see Fig. 3c]. The~0.3 and 0.7 eV features, which are more evident in experiments, correspond to d − d excitations between the e g − e g and e g − t 2g orbitals, respectively. These excitations contain both spin-flip and non-spin-flip contributions, which can be further distinguished through polarization control (See Fig. 3d and SI). Due to the locality of these low-energy d − d excitations, their energy scales are primarily determined by crystal field and lattice geometry instead of in-plane hoppings or Coulomb interactions. Therefore, the energy scales of these excitations are robust against varying T c . We notice that RIXS measurements on FeS complexes with a distorted tetrahedral CFE for Fe 2+ also observed two peaks at 0.32 and 0.58 eV attributed to d − d excitations 39 . RIXS experiments on PrFeAsO 0.7 (T c = 42 K) revealed a~0.5 eV feature assigned to d − d orbital excitation as well 40 .
In addition to these two low-energy peaks, our simulation exhibits substantial intensity in a 2 eV window around the energy loss ω~2.5 eV, which consists of multiple spectral peaks when the spectral broadening is reduced. These excited states mainly exhibit a change of total spin quantum number ΔS = 1. The energy scale for these high-energy excitations is mainly determined by the Racah parameters B and C between d − d orbitals. In our calculation, we used B = 0.12 eV and C = 0.4 eV, which lead to a typical Hund's coupling~0.8−0.9 eV for Fe 2+ . When only a single Fe site is present, the Racah parameter A provides only an overall energy shift of the atomic multiplet levels, and it is simply set to zero here. Since B and C arise from multipole interactions, they are more difficult to be screened (unlike A due to monopole interaction). Therefore, it is expected that the high-energy 2.5 eV feature will not be strongly affected by charge screening and the change of T c , which is also consistent with the experiment.

Temperature dependence
To further explore the nature of the Raman excitations, we performed temperature-dependent measurements from 25 to 275 K on both SC28K and SC42K (Fig. 4). It is clear that the intensity at the energy-loss region of 0.  Fig. 1c at 20, 95, and 200 K, respectively. Each spectrum is shifted vertically for clarity. The spectra were collected at the Fe L 3 resonant peak energy of 709 eV. Spectra were normalized to the high energy-loss range of [4,12] eV. The elastic peak was subtracted from the data to visualize the low-energy features. b, d, f Similar measurements on SC42K. Examples of the six peak decompositions described in the text are shown in the bottom spectra. g, h, i Summary of fitted energies of four Raman features along two high-symmetry directions for both SC28K and SC42K at 25, 95, and 200 K, respectively. To avoid specular reflection, we did not measure RIXS spectra at zero momentum transfer. The size of the markers corresponds to the spectral weight of the features. Error bars were the standard deviation of the fits. temperature (Fig. 4a, c). Moreover, the ratio of the intensity between two excitations at~0.3 and~0.7 eV varies with temperature (Fig. 4b,  d), which supports our fitting with two peaks (details of the fitting shown in Supplementary Fig. 8). The temperature dependence of the integrated intensity of the four Raman features is shown in Fig. 4e-h. Intriguingly, the excitations at~0.3 and~0.7 eV decrease by more than 50% from 25 to 275 K, while the features at~0.1 and 2.5 eV change only slightly with temperature in both samples. It is also known that the change of T c is connected to different Fe vacancy concentrations [32][33][34] . This can affect the oxidation state so that Fe valency may be reduced from 2+, which can slightly affect the ratio between the 0.3 and 0.7 eV peaks. The temperature dependence of the 0.3 and 0.7 eV spectral features suggests the primary orbital nature of these excitations. Previous magnetic susceptibility measurements on SC28K and SC42K displayed a drop or deviation from the Curie-Weiss behavior around~120 K, implying a magnetic transition or fluctuation 17,33 . In Fe-based superconductors, magnetic transitions or fluctuations are usually accompanied by an orbital weight redistribution. Previous studies have demonstrated that the spectral intensity can exhibit a more apparent temperature dependence when orbital fluctuations are present 41 . To demonstrate the relevance of orbital fluctuations, we simulate the dynamical spin and (non-spin-flip) orbital structure factors S(q = 0, ω) and O(q = 0, ω) separately (see Fig. 3d). Orbital excitations with spin-flip character is included in the S(q, ω). As expected from the single-particle picture in Fig. 3c, the two low-energy Raman peaks at~0.3 and~0.7 eV are present in both orbital and spin channels. A full-polarization study by measuring the polarizations of both incident and outgoing x-rays would help us to quantitatively determine the orbital and spin contributions to the peaks at~0.3 and 0.7 eV. Such full-polarization investigation is suggested for future more advanced studies. Nevertheless, it is possible to determine qualitatively the orbital and spin contributions by performing temperature-dependent measurements. Considering the experimental measurements integrated over all scattering polarizations, the contributions from these two channels are approximated at 2:1. Since the onset of orbital excitations can cause a rapid transfer of the spectral weight away from coherent peaks 41 , a rough estimation is that the intensity of these d − d excitations drops by 60% at high temperature due to the fluctuation of orbitals, consistent with the experimental observations for SC28K and SC42K. Moreover, we notice the suppression of the~0.7 eV peak involving the d xy orbital is consistent with ARPES measurements on monolayer FeSe that show the disappearance of d xy spectral weight at high temperature 12 [see Fig. 4g]. This agreement implies substantial orbital fluctuations and orbital redistribution in high-T c iron chalcogenides. We would like to stress that the two low-energy peaks at~0.3 and 0.7 eV are both orbital excitations no matter whether the spin-flip is involved. Therefore, the mixture of spin channels does not weaken our conclusion about orbital fluctuations.
In contrast, the orbital O(q = 0, ω) has almost no contribution to the high-energy regime. Thus, the 2.5 eV Raman feature originates primarily from local spin-flip d − d excitations, accounting for its weak temperature dependence. This mechanism can be further tested in the future by the temperature dependence of RIXS spectra employing a polarimeter for full-polarization measurements. In addition to the above three features at higher energies, another Raman feature at~0.1 eV is also dispersionless within our energy resolution (~80 meV), which is absent in the multiplet model calculation. ARPES and electron energy-loss spectroscopy (EELS) measurements on monolayer FeSe reported a phonon mode at~90-100 meV arising from the STO substrate 18,20,42 . However, our FeSe11111 thin-film has a different substrate LaAlO 3 and the thickness of the film is several hundreds of nanometers. A previous Raman study on FeSe11111 measured phonons only up to 400 cm −1 43 . Nevertheless, the presence of oxygen and hydrogen within (Li/Fe)OH layers may contribute to a 0.1 eV phonon excitation. Another possibility relates to an optical spinwave mode, which is typically much less dispersive compared to the acoustic spin-wave mode. In this case, the system exhibits short-or long-range magnetism involving block spins in an enlarged magnetic unit cell 44 . However, INS studies of FeSe11111 have not reported an optical magnon around 0.1 eV 45 , which renders this origin less likely.

DISCUSSION
The dispersionless excitations in FeSe11111 and monolayer FeSe contrast the dispersive spin excitations observed in bulk FeSe, reflecting their intrinsic differences. Nevertheless, the absence of orbital excitations in bulk FeSe might be due to these orbital excitations being broader and therefore harder to isolate from the large fluorescence background in bulk FeSe 28,29,31 , which remains to be explored in future studies. Our theory-experiment comparison suggests that a minimal multiplet model involving the Fe 3d orbitals can quantitatively address these dispersionless features. Informed by experimental results, our simulation ignores nonlocal interactions, which usually cause dispersive magnon excitations. That said, an orbital-dominated mechanism does not contradict and may coexist with the recently proposed scenario of hybridized magnons in a two-band Hubbard model 31 . We notice that previous INS studies reveal the spin excitations dispersing up tõ 140 meV near [π, π] in FeSe11111, but their spectral weight is much lower than that in bulk FeSe due to doping effects 45,46 , which would make it difficult to detect by RIXS within the accessible momenta near [0, 0]. Our observation of orbital excitations and strong temperature-dependent orbital fluctuations is still compatible with Hund's metal picture, which can lead to an orbital-selective Mott phase at zero temperature. The degrees of orbital fluctuations at elevated temperatures will be dependent on the competition between different energy scales, such as Hund's coupling and crystal-field splitting.
Orbital-selective Cooper pairing, manifesting as highly anisotropic superconducting energy gaps, has been discovered in both bulk FeSe and monolayer FeSe using Bogoliubov quasiparticle interference imaging 47,48 . By raising the temperature, ARPES studies have observed strong orbital-dependent correlation effects in monolayer FeSe where the d xy orbital completely loses its spectral weight 12 . However, this effect is not found in bulk FeSe 49 . Here we observe that FeSe11111 displays similar orbital excitations as those in monolayer FeSe 31 . Moreover, we find the spectral weights of orbital excitations at~0.3 and 0.7 eV decrease with increasing temperature, and the latter involving d xy orbitals shows stronger temperature-dependent behavior, consistent with ARPES results 12 . Since the evident change of orbital spectral weights occurs not only in FeSe11111 but also in monolayer FeSe, both exhibiting a relatively high transition temperature T c , it suggests that the orbital fluctuations may play important roles in superconductivity for these materials.
In our FeSe11111 material, the prominent orbital excitations at 0.3 and 0.7 eV are attributed to d − d transitions among the e g − e g and e g − t 2g orbitals. In iron pnictides, it was shown that the t 2g orbitals are necessary to lead to the ferro-orbital fluctuations 9 . The lattice-induced orbital fluctuation effect, in turn, can favor superconductivity. Moreover, when the electron-phonon coupling matrix elements for all five 3d orbitals are included, the critical value of phonon-mediated interactions that separates the phases of orbital fluctuation and orbital order is reduced to half compared to that in the three-orbital model. In these regards, all five 3d orbitals and their fluctuations can be important to superconductivity. This mechanism is consistent with the previous scanning tunneling microscopy (STM) studies in both monolayer FeSe and FeSe11111 50,51 , in favor of the phonon-mediated pairing enhancement scenario 18 or orbital fluctuation-mediated pairing mechanism 9 . The absence of the d xz /d yz hole Fermi surface may also suppress the nematic order that competes with superconductivity 52 , underscoring the crucial role of the orbitals in the pairing.
We stress that the orbital excitations referred to in this work involve both spin-flip and non-spin-flip excitations. Therefore, our conclusion does not exclude the role of spin fluctuations in the pairing mechanism, which has also been proposed in FeSe11111 45,53,54 . For example, the phase-sensitive quasiparticle Q. Xiao et al. interference technique has revealed a sign-changing order parameter in Zn-doped FeSe11111 53 and INS studies have observed the spin resonance below T c in FeSe11111 45,54 , indicating the important role of spin fluctuations in electron pairing. Thus, as a layered material without STO substrate, FeSe11111 and its excitations provide a unique opportunity to elucidate the evolution from bulk to monolayer FeSe materials, which suggests the importance of both orbital and spin degrees of freedom for high-T c iron selenides.

Sample preparation
The high-quality single-crystalline FeSe11111 films were synthesized by the matrix-assisted hydrothermal epitaxy (MHE) technique 55,56 . Their crystalline quality was characterized by x-ray diffraction (XRD) on a 9 kW Rigaku SmartLab x-ray diffractometer equipped with two Ge(220) monochromators at room temperature ( Supplementary Fig. 1). Both the XRD patterns exhibit a single preferred (00l) orientation and the full width at half maximum (FWHM) is around 0.3 ∘ of the x-ray rocking curves for the (006) reflection, showing the highly crystalline quality of these films (Supplementary Fig. 1). Their T c values are defined at the onset temperatures of the diamagnetism by magnetic measurements (Supplementary Fig. 1), performed with a 1 Oe magnetic field applied along the c-axis using a SQUID magnetometer (MPMS XL1, Quantum Design).

RIXS measurements
We performed high-resolution RIXS experiments by using the Super Advanced x-ray Spectrometer (SAXES) at the Advanced Resonant Spectroscopy (ADRESS) beamline of the Swiss Light Source, Paul Scherrer Institut, Switzerland 57,58 . The energy resolution was 80 meV (FWHM). We cleaved both samples in an ultra-high vacuum. We used the grazing-incident geometry, which provided stronger signals than the grazing-out geometry ( Supplementary Fig. 3). All samples were aligned with the surface normal (001) in the scattering plane. π polarized x-rays were used for the RIXS measurements since the signals were similar between πand σpolarizations ( Supplementary Fig. 3). The incident energy was tuned to the peak maximum of Fe L 3 -edge (709 eV), except for detuning measurements. The angle between incoming and outgoing x-rays was fixed to 130 degrees, resulting in a constant value of the total momentum transfer. We refer to reciprocal space coordinates in the conventional one-Fe unit cell with a = b = 2.680 Å, c = 9.318 Å for SC42K and c = 9.261 Å for SC28K 33 . The momentum transfer Q is defined as Q = Ha * + Kb * + Lc * where a * = 2π/a, b * = 2π/b, and c * = 2π/c. q is the projection of Q in the abplane. We measured ab-plane momentum transfer scans by varying the grazing-incident angle to the sample surface. Each RIXS spectrum was measured for 30 min except for the spectrum of incident energydependent RIXS maps, which was measured for 5 min.

Atomic multiplet theory
To model the experimental spectra, we make two assumptions about the Hamiltonian and the RIXS cross-section. First, we consider only multiplet interactions in the atomic limit 59 . In this case, our calculation can capture local spin or orbital flip excitations but not collective spin waves (magnon). Nevertheless, this assumption is appropriate in our system, as (i) the experiment exhibits little momentum dependence and (ii) we are not aware of any lattice model with short or long-range magnetism capable of exhibiting a complete flat spin-wave dispersion. The resulting multiplet Hamiltonian containing all Fe 3d orbitals reads Here, d ðαÞ σ (d ðαÞy σ ) annihilates (creates) a 3d α electron with spin σ (α = x 2 − y 2 , z 2 , xy, yz, and xz). ε α denotes the site energy of orbital α. In a distorted tetragonal CFE, ε α can be parameterized by three parameters 10Dq, δ, and μ 60 . We use 10Dq = 0.485 eV, δ = −0.12 eV, and μ = −0.32 eV. Here, 10Dq represents the site energy separation between e g and t 2g orbitals, and its typical value is~0.5 eV for a tetrahedral CFE in iron-based superconductors. δ and μ describe additional splitting between the e g and t 2g orbitals due to a lowered crystal symmetry. The intra-orbital U and interorbital U 0 interactions, as well as Hund's coupling and pair hopping J, can be parameterized by the Racah parameters A, B, and C. In the case of a single Fe atom, A only provides an overall energy shift, so it is simply set to zero. Without any fine-tuning, we use B = 0.12 eV and C = 0.40 eV, leading to a typical value of~0.8-0.9 eV for Fe 2+ Hund's coupling, independent of charge screening and the change of T c . The Hamiltonian is solved by exact diagonalization (ED) to obtain the energy eigenvalues and eigenstates. It is noted that within reasonable CFE and interaction parameters, a high spin S = 2 ground state is always favored, which agrees with our further ab initio CASPT2 simulations (see below). If the ground-state spin quantum number is reduced from S = 2, then the spectral intensity ratio between the 0.3 and 0.7 eV peaks also will change 61 . This can happen with charge or orbital fluctuation effects, which are always present to some extent in the experiment but not captured in the atomic multiplet calculation.
The RIXS cross-section is defined as where and the (momentum-resolved) dipole transition operator D k ¼ P i e ÀkÁri D i and E g is the ground state energy. In principle, the intermediate state of a RIXS process contains a core-hole interaction (with core-hole lifetime 1/Γ) H core , describing the interaction between Fe core 2p and valence 3d orbitals. To simplify the calculation and clearly assign the spectral features, we adopt the second commonly used assumption of ultrashort core-hole lifetime approximation 62,63 , i.e., Γ ≫ U α . In this case, the RIXS cross-section is reduced to the dynamical spin and orbital structure factors. The dynamical spin structure is defined as . These two structure factors approximate the spin-flip and non-spin-flip RIXS at the ultrashort core-hole lifetime limit. We note that the charge structure factor N(q, ω) also contributes to the nonspin-flip cross-section, but its inelastic response vanishes for q = 0 due to charge conservation. Using the atomic multiplet model, we only present the response function for q = 0 in this paper.

Ground-state configuration from ab initio quantum chemistry
To determine the ground-state spin configuration as a starting point for the ED calculation, we perform two different ab initio simulations for a single unit cell Fe(II)[Se 2− ] 4 complex. Due to the expectation of strong multi-reference effects, we complement the PBE0 DFT simulation with a CASPT2 simulation with (10e, 12o) active space. As shown in Table 1, both simulations conclude with a ground state with S = 2, justifying the selection of CFE in our atomic multiplet simulations. More specifically, the DFT simulations predict relatively large energy differences among all spin sectors compared to the CASPT2 results. This is because the lower-spin state exhibits a stronger multi-reference effect as a result of quantum fluctuations. Without accounting for this effect, DFT over-estimated the energy of these lower-spin states. In contrast, the energy differences obtained by CASPT2 are relatively small, accounting for the experimental observation of reduced magnetic moment with the presence of itinerancy.