Ab initio Derivation of Correlated Superatom Model for Potassium Loaded Zeolite A

We derive an effective low-energy Hamiltonian for potassium loaded zeolite A, a unique ferromagnet from non-magnetic elements. We perform ab initio density functional calculations and construct maximally localized Wannier functions for low-energy states made from potassium s electrons. The resulting Wannier orbitals, spreading widely in the alminosilicate cage, are found to be the superatomic s and p orbitals in the confining potential formed by the host cage. We then make a tight-binding model for these superatomic orbitals and introduce interaction parameters such as the Hubbard U. After mean-field calculations for the effective model, we find that ab initio spin density functional results are well reproduced by choosing appropriate sets of the interaction parameters. The interaction parameters turn out to be as large as the band width, $\sim$ 0.5 eV, indicating the importance of electron correlation, and that the present system is an interesting analog of correlated multi-orbital transition metal oxides.


Introduction
Electron correlation in materials, which often causes various non-trivial phenomena, has been one of the central issues in condensed matter physics. While there are many kinds of fascinating strongly correlated systems such as organic molecular p-, transition metal d-, or heavy fermion f -electron systems, it is of great interest to consider the possibility of correlated s-electron systems. Since valence s-electrons are usually itinerant, we expect that correlations in s-electron systems cannot be so significant. However, recently, alkali-metal loaded zeolite is attracting broad interests as a unique exception; when clusters of alkali atoms are confined in the cage of alminosilicate, their kinetic energy is drastically suppressed and the system can exhibit a variety of correlation effects.
Although the unitcell of alkali-metal loaded zeolites is huge and complicated [typically it contains O(100) atoms], the low-energy electronic structure of zeolites is surprisingly simple. The host alminosilicate cage makes a gap of several eV, in which the guest alkali s-electrons form a few narrow bands around the Fermi level. This situation is described by the so-called superatom picture, [1][2][3] where each zeolite cage is regarded as a huge atom. Here, the confining potential formed by the host framework and the guest alkali-metal s electrons correspond to the atomic potential and valence electrons of the superatom, respectively. Recent ab initio calculations have elucidated that the band dispersion around the Fermi level is well represented by a very simple tightbinding model, [3][4][5][6][7][8] suggesting that correlation effects in low-energy physics of zeolites can be systematically described by the many-body superatom model.
Each cage has only one valence s electron, so that it can be viewed as a crystal of hydrogen-like superatoms. In ref. 4, a single-orbital extended Hubbard model for the superatomic s electrons has been constructed, and the correlation strength was found to be much larger than the band width.
On the other hand, we also have many intriguing superatomic multi-orbital systems in the family of zeolites. Potassium loaded zeolite A (K 16 Si 12 Al 12 O 48 , hereafter we call it K-LTA, where LTA is the code for the crystal structure of zeolite A) is a typical example, for which a spin-polarized ground state is realized. 1,2) While the mechanism of spin-polarization has not been fully understood, there are several experimental [9][10][11][12] and theoretical 3,13) indications that orbital degrees of freedom of the superatomic p orbital is important. In this paper, we derive an effective multi-orbital Hamiltonian of K-LTA from first principles. We estimate onsite potentials, transfer integrals, and interaction parameters in the Hamiltonian. We show that the energy scale of the electron correlation, estimated as ∼ 0.5eV, is as large as the band width, so that the system is strongly correlated and can be regarded as a unique "miniature" of transition metal oxides.
The structure of the paper is as follows. In § 2, we first briefly describe the atomic configuration of K-LTA and give the level diagram of the low-energy states, which is useful to discuss the band structure of K-LTA in terms of the superatom picture. In § 3, we describe the detail of the method to evaluate the parameters in the effective Hamiltonian. The basic idea is following: We first perform density functional calculation and make maximally localized Wannier functions 14) for several bands around the Fermi level. Using these Wannier functions, we then construct a tight-binding model and introduce interaction parameters such as the Hubbard U . Choosing appropriate sets of the interaction parameters, we perform a mean-field calculation for the effective model and fit to the result of ab initio spin density functional calculation. In § 4, we present our results which indicate that the system is indeed classified as a strongly correlated system. Section 5 is devoted to a summary and outlook.

Level Diagram of Superatomic p-orbitals in Zeolite A
Let us first look at the atomic configuration of K-LTA. The system has two kinds of alminosilicate cage, the socalled α and β cage (see Fig. 1 and related descriptions in ref. 13). While the diameter of the β cage is ∼ 7Å, the α cage is much larger (∼ 11Å), so that the guest potassium cluster is accommodated only in the α cage. In the unitcell, there are two inequivalent α and β cages, which are denoted as α 1 , α 2 , β 1 , and β 2 . Concerning the atomic positions of potassium atoms, neutron scattering measurements 15) revealed that there are four kinds of sites: K I , K II , K III , and K IV [ Fig. 1 (a)]. Here, the K I , K II , and K III atoms sit around the center of the six-, eight-, and four-membered rings in the alminosilicate network, respectively, and the K IV atom sits at the center of the α cage. The sites of K I and K II are fully occupied, while the K III and K IV sites are partially occupied. Following the argument in ref. 13, we assume that the system has three-fold symmetry.
The low energy bands around the Fermi level are mainly made from the superatomic p orbitals in the α cages. When we discuss the superatomic p levels, there are two key parameters. One is the difference between the potentials of the two α cages, defined as ∆I = min{ǫ pxα1 ,ǫ pyα1 ,ǫ pzα1 }−min{ǫ pxα2 ,ǫ pyα2 ,ǫ pz α2 }(1) with ǫ pµαi specifying the ionization potential of the superatomic p µ orbital in the α i cage. Since the α 1 cage accommodates more K + ions than the α 2 cage, we naively expect that the potential of α 1 is deeper than that of α 2 . However, it should be noted here that the K I atoms [denoted by solid (blue) circles in Fig. 1 (b)] are distributed sparsely in the α 1 cage than in the α 2 cage, and the α 1 cage does not have K IV at the cage center [solid light gray (pink) circles in Fig. 1 (c)]. These facts make, on the contrary to the above expectation, the α 1 -cage potential shallower than that of α 2 [the schematic diagram in Fig. 1 (b)].
The other important parameter is the crystal field splitting which is controlled by the atomic configuration of K III [solid dark gray (green) circles in Fig. 1 (c)]. For the K III configuration spreading in the xy plane [(i) and (iii) in the figure], the degenerated p x and p y levels become lower (lower schematic diagrams). In contrast, for the configuration which extends in the z direction [(ii) and (iv)], the crystal field lowers the p z level. Note that the α 1 cage contains 6 K III atoms and the α 2 cage does 3 K III atoms. Also note that the K IV atom is accommodated only in the α 2 cage. The crystal-field splitting ∆V αi in the α i cage is defined as , and (c) Crystal-field splitting (∆Vα i ) formed by clusters consisting of K III and K IV [denoted by solid dark gray (green) and light gray (pink) circles, respectively]. The z axis is taken to be perpendicular to one of the six-membered rings in the α cage, and we assume that the system has threefold symmetry along the z axis. Note that the sign and size of the crystal-field splitting depends on the atomic configuration of the cluster. Four types of the K III and K IV configurations considered in the present study are depicted.
In the previous study, 13) we performed spin density functional calculations for the four geometries in Fig. 2, which we call geometries I-IV. The figure also displays the schematic level diagrams with the superatomic pic-ture. It was found that the geometry I, whose total energy is the lowest among the four geometries, has the largest magnetic moment (∼ 2µ B ) in the α 2 cage. On the other hand, the magnetic moment is smallest for the geometry II. This geometry dependence was simply discussed in terms of ∆I and {∆V αi } in ref. 13. While these are key parameters for the formation of local spin in the α cage, transfers between the neighboring α cages is also important to discuss the magnetism in the bulk. We should also consider the degrees of freedom of the β cage neglected in the previous analysis, mediating the hopping between the superatomic states in the α cage. On top of that, careful analysis on the interplay between electron itinerancy and correlations are imperative. Thus in the present study, we evaluate the interaction parameters such as Hubbard U for the superatomic orbitals.

Effective Model Construction
Let us move on to the detail of the method to derive an effective model from first principles. For evaluation of interaction parameters such as the Hubbard U , in the previous study for alkali-metal loaded sodalite, 4) we used the constrained random phase approximation. 16,17) However, since the system size of K-LTA is much larger than that of sodalite and the numerical cost is so expensive, here we take a different approach.
First, we introduce the following model Hamiltonian, where c † imσ (c imσ ) creates (annihilates) an electron with spin σ in orbital m at site i, and n imσ ≡ c † imσ c imσ . The term H 0 describes the kinetic-energy part of the effective model, where ǫ im and t imjm ′ are the ionization potential and transfer, respectively. In this model, we take into account of not only the superatom p states in the α cage but also the superatom s states in the β cage, which lie just above the Fermi level as explained in the next section. The term H int is the interaction part, i.e., the intraorbital (U im ) and interorbital (U ′ imm ′ ) Coulomb interactions and the Ising component of Hund's coupling (J imm ′ ). Since the standard spin density functional calculation is not SU(2) symmetric, we ignore the spin-flip component of the Hund's coupling and the pair-hopping interaction. For simplicity, we also ignore orbital dependence of interaction parameters in the α cage; common U α is used for any m and common U ′ α and J α are used for any combination of m and m ′ . For the β cage, there is no orbital degree of freedom. Now, let us consider to derive the interaction parameters {U α , U β , U ′ α , and J α }. These parameters are determined by the fitting procedure in such that the Hartreelike mean field solution of the model reproduces both of ab initio non-magnetic and spin-polarized density functional results based on the generalized gradient approximation (GGA).
In order to relate the mean field in the model calculation and the density functional calculation, we introduce the following two terms. One is a term to scale the self interaction (SI) and the other is a term which guarantees that the non-magnetic model calculation gives the same result as the non-magnetic density functional calculation. For the first issue, as is well known, GGA includes the SI effect. 18) On the other hand, the Hubbard model defined in eq. (3) does not have SI; the intra-orbital interaction between electrons with the same spin are excluded in the model Hamiltonian. Therefore, to mimic the GGA situation in the model calculation, we need to introduce the SI term as where Thus the one-body part H 0 in eq. (4) is defined as where H GGA is the tight-binding representation of the non-magnetic GGA Kohn-Sham Hamiltonian with ǫ GGA im and t GGA imjm ′ in diagonal and off-diagonal parts, respectively. We use the maximally localized Wannier functions for the basis functions for those matrix elements. More precisely, the matrix elements of H 0 are calculated as and Here, ρ GGA imm is the matrix element of the non-magnetic GGA density matrix and N GGA i ≡ m ρ GGA imm . Note that transfers in H 0 is the same as those in H GGA , since H SI andH are local in real space. Now, the mean-field Hamiltonian H MF [ρ] for H reads as Note that the result of non-magnetic GGA is exactly reproduced as H MF [ρ GGA ] in the present framework.
We then consider to solve H MF by spin-polarized selfconsistent Hartree scheme. We can estimate the interac- tion and self-interaction parameters, U α , U β , U ′ α , J α , S α , and S β , by fitting the model band dispersion to the ab initio spin density functional results.

Results
Using the Tokyo Ab initio Program Package, 19) we first performed density functional calculations with the GGA exchange-correlation functional. 20) We used a plane-wave basis set with the ultrasoft pseudopotentials. 21) The energy cutoff was set to 36 Ry for wavefunction and 144 Ry for spin and charge densities and 4×4×4 k-mesh was adopted. Maximally localized Wannier function with the ultrasoft pseudopotential was constructed following the recipe of ref. 22.

Maximally localized Wannier functions
The GGA band structures for the geometries I-IV are shown in Fig. 3 as (red) dotted lines. We see that the low-energy band structures depend sensitively on the geometry. Especially, the bandwidth is rather different from each other. From the six (seven) bands around the Fermi level, we constructed the maximally localized Wannier functions for the geometry I and II (III and IV), which are employed as bases of the effective Hamiltonian in eq. (3). We estimated tight-binding parameters for the geometries I-IV by calculating matrix elements of H GGA in the Wannier basis. In Fig. 3, we plot interpolated band dispersions [(green) solid lines], from which we see that original ab initio bands are well reproduced by the Wannier interpolation. The band around −0.4∼−0.3 eV is formed by the superatomic s orbital in the α cage, which is not included in the effective model.
As shown in detail later (Fig. 6), the s orbital in the α cage is fully occupied and its exchange splitting is negligibly small. Thus, for simplicity, we do not consider explicitly this degree of freedom in the low-energy model. It should be noted that the s state in the α cage has small transfer integrals between the surrounding sites and tiny onsite hybridizations with the p states, because it is well localized in the cage and the confining potential is almost spherical. Moreover, the s state should have larger U and S than the p states, so that there should be a significant difference in the ionization potential, ǫ im [eq. (7)], for the s and p states. Therefore, the s state is expected to be irrelevant for the magnetic properties.
In Fig. 4, we display the isosurface contours of the plane-wave part of the Wannier functions for the geometry I. In this geometry, there are six bases (sβ 1 , sβ 2 , p z α 1 , p z α 2 , p y α 2 , and p x α 2 ), and the p x α 1 and p y α 1 states are not included in the effective model, because the latter two exist far above the Fermi level by, at least, ∼0.5 eV. Similarly, for the geometries II-IV, the high-energy p µ α 1 states were not included in the effective models. Table I lists our calculated spatial spreads of the Wannier functions, Ω mi , for all the geometries. As is indicated in the table, the size of the Wannier functions is as large as the diameter of the host cage (∼11Å for the α cage and ∼7 A for the β cage) and Ω β ≤ Ω mα , so that the onsite in-teractions in the β cages should be larger than those in the α cages.  Table I. Spatial spread of maximally localized Wannier functions for superatomic states, Ω, and self-consistently determined ionization potentials, ǫ defined by eq. (7). The difference in cage potentials ∆I [eq. (1)] and the crystal-field splitting in the α 2 cage ∆Vα 2 [eq. (2)] are also listed. The unit is given asÅ for Ω and eV for ǫ, ∆I, and ∆Vα 2 .  On the other hand, the transfers between neighboring α cages are |t pz αpz α | = 0.01 ∼ 0.02eV, |t pzαp x/y α | = 0.02 ∼ 0.03eV, and |t p x/y αp x/y α | = 0.03 ∼ 0.04eV. The distant transfers are negligibly small. We remark that, among the transfers between α cages, the p z -p z hopping is smaller than those of other p states. This is because the p z orbital is not pointing at the neighboring α cages. Thus, the p z states tend to be localized and to have large magnetic moments.
In order to estimate {U α , U β , U ′ α , J α , S α , S β }, we performed spin-polarized self-consistent calculations for H MF in eq. (9), and looked for the parameter set which reproduces the band dispersion of ab initio spin density functional calculations. First, we searched for the parameter set for which the difference between the GGA band dispersion of the geometry I and that of the model calculation along the symmetry line W→X→ Γ →L→ K in the Brillouin zone is minimized. We assumed that U α < U β < 0.6eV, S α < S β < 1, U ′ α /U α < 1, and J α /U α < 0.2 with the realistic parameter range. 4) To understand and explain this system by a minimum set of parameters, we then introduced a scaling factor x for {U α , U β , U ′ α , J α } and minimize the band-dispersion difference in GGA and model for the geometry II, III, and IV. We found that the GGA band structures for these three geometries are successfully reproduced for x = 0.6, S β = 0.7, and S α = 0.1. These parameters show that Geom. II-IV have smaller interaction parameters than Geom. I. This is consistent with the fact that the band widths in Geom. II-IV are broader than that in Geom. I as shown in Fig. 3.
Before presenting the resulting interaction parameters, we summarize in Table I our calculated self-consistent non-magnetic ionization potentials, ǫ in eq. (7). Using the ǫ data, we estimated the effective cage-potential differences ∆I in eq. (1) and the crystal field splitting in the α 2 cage ∆V α2 in eq. (2), which are given in the bottom two rows in Table I. 23) It should be noted here that these one-body quantities have appreciable geometrical dependence, suggesting that slight changes in the potassium configurations can cause substantial changes in the low-energy electronic structure of K-LTA. According to Table I, we draw in Fig. 5 the level diagram of the superatomic states in the α and β cages. We can see that the α 2 -cage potential is significantly deep compared to the other cage potentials and thus will form the trapping potentials for the low-energy electrons.
We next show in Fig. 6 our calculated spin-polarized model band dispersion of the majority spin [(red) solid lines] and minority spin [(green) dashed lines], comparing with those of the ab initio results [majority (+) and minority (×)]. Table II compares the resulting local magnetic moments in the model and ab initio calculations. In spite of the limited number of the free parameters, the ab initio band dispersion and moment size are well reproduced. 13) The magnetic moments are basically formed in the α 2 cage, being consistent with the level diagram in Fig. 5. The transfers cause the antiferromagnetic correlation between the α 1 and α 2 cages.
The interaction parameters estimated by the present fitting are summarized in Table III. We see that the re-sulting U is 0.24-0.5 eV and is always about ten times larger than the transfer integrals as 0.01-0.04 eV, thus indicating that the system resides in moderately/strongly correlation regime. We note that interaction parameters depend sensitively on the atomic configuration; they tend to be large (small) for the geometry I (II-IV). 24) It is also interesting to note that the self-interaction parameter S i tends to be larger when the interaction parameters are large. This is a reasonable trend, 18) because self-interaction is generally large for localized states, and localized states usually have strong onsite electron correlation. In the present study, the S i value is estimated by parameter fitting to the spin polarized GGA calculations. Therefore, the present S i parameter qualitatively measures the localization of the spin-polarized wavefunctions, but not that of the paramagnetic wavefunctions. Finally, let us comment on the controversial situation in the experiment. Experimentally, two kinds of possible magnetic structures have been proposed for ferromagnetism of K-LTA. One is the ferrimagnetic model, 9) and the other is the spin-canted antiferromagnetic model. 11) In the former, the system consists of two sublattices with large (∼ 2.8 µ B ) and small (almost 0.0 µ B ) magnetic moment, while, in the latter, every α cage has a moment of 1 µ B and the spins are canted in the low temperature by the spin-orbit interaction. These two magnetic structures seem to agree with the results for the geometries I and IV, respectively, 13) although the size of the moments is relatively smaller than the experimental values. It is interesting to note that the geometries I and IV are energetically stable (slightly lower in the geometry I by ∼0.1  eV). 13) The present results seem to support the former proposal, but we should emphasize that K-LTA resides in moderately or strongly correlation regime so that careful analyses with proper consideration of the correlation effect are required to clarify the mechanism of the spinpolarization in real K-LTA. Indeed, strongly correlated multi-orbital systems generally have rich phase diagrams. For example, a subtle competition between the interaction parameters and the tight-binding parameters can drastically change the magnetic properties of the system. Elaborated many-body analyses for K-LTA are indeed a fascinating future problem.

Summary and Outlook
We have derived an effective low-energy Hamiltonian for K-LTA. First, we constructed a tight-binding model from GGA using the maximally localized Wannier functions. We then introduced the interaction and self-interaction parameters, and solved the model selfconsistently. With appropriate sets of these parameters, we have succeeded in reproducing the result of spin density functional calculation. We have found that the system reside in the strong coupling regime in that the bandwidth and electron correlation have the same energy scale. We have also shown that the details of the potassium-cluster configurations crucially affects both the one-body and two-body parts of the effective Hamiltonian of K-LTA.
While the present method of evaluating the interaction parameters has a great advantage in that it can be applied to huge systems such as K-LTA, it does not work for systems for which spin density functional calculation gives a non-magnetic ground state. Since the size of magnetic moment tends to be underestimated in GGA, the present method is expected to give the lower bound of the interaction parameters. It is indeed an interesting future challenge to apply other ab initio methods to evaluate interaction parameters such as constrained density functional theory (ref. 25) or constrained randomphase approximation (refs. 16 and 17) to K-LTA. While formidable numerical cost will be required, comprehensive model construction by different methods is important to establish the low-energy correlation physics in the zeolite system.
We thank Professor Yasuo Nozue and Takehito Nakano for fruitful discussions. This work was supported by Grants-in-Aid for Scientific Research (No. 19051016, 22740215, 22104010, and 23110708) and the Next Generation Super Computing Project from MEXT, Japan, and "Funding Program for World-Leading Innovative R&D on Science and Technology (FIRST program)" JSPS, Japan and the JST PRESTO program, Japan. Numerical calculations were done at the ISSP supercomputer center, and the Supercomputing Division, Information Technology Center, University of Tokyo.