# Correlated electron behavior of metalorganic molecules: insights from density functional theory combined with many-body effects using exact diagonalization

###### Abstract

A proper theoretical description of electronic structure of the 3d orbitals in the metal centers of functional metalorganics is a challenging problem. In this letter, we apply density functional theory and an exact diagonalization method in a many body approach to study the ground state electronic configuration of an iron porphyrin (FeP) molecule. Our study reveals that dynamical correlation effects are important, and FeP is a potential candidate for realizing a spin crossover due to a subtle balance of crystal field effects, on-site Coulomb repulsion and hybridization between the Fe d-orbitals and ligand N p-states. The mechanism of switching between two close lying electronic configurations of Fe-d orbitals is shown. We discuss the generality of the suggested approach and the possibility to properly describe the electronic structure and related low energy physics of the whole class of correlated metal centered organometallic molecules.

Molecular magnets combine low dimensionality and inherent confinement effects with strong electron correlations and hold prospects in the context of spintronics. An important molecular property is bistability, i.e., the possibility of realizing two different spin states, which can in principle be accessed and manipulated externally. This is important as the switching of spin has a pronounced effect on measurable quantities, like magnetic anisotropy and spin dipole moment contribution sumantaprl (). Finding ways for the efficient manipulation of the magnetic state spinswitch1 (); sumantachemi (); sumantaprl (); nsr () of transition metal (TM) centered porphyrin (TM-P) and phthalocyanine (TM-Pc) molecules have critical consequences in this regard. A crucial interplay between molecular ligand field and spin pairing energy makes only a subspace of this class of materials to respond to spin crossover.

The magnetic properties in TM-P/TM-Pc are largely governed by the metal center, which features sizable local Coulomb interactions ( eV and eV) and is simultaneously subjected to crystal fields, spin-orbit coupling and orbitally dependent hybridization with the ligands. Electronic correlations are expected to arise Fazekas () and the description by local density approximation (LDA) or generalized gradient approximation (GGA) thus potentially becomes inadequate, for example, leading to underestimated or even vanishing HOMO-LUMO (HOMO=Highest Occupied Molecular Orbital and LUMO=Lowest Unoccupied Molecular Orbital) gap. Hence, the treatment of the molecular electronic structure in terms of correlated electron theories like ligand field theories or Anderson impurity models become crucial. These model based approaches can be very helpful to trace the physical origin of phenomena like spin-state switching Fazekas (), emergence of magnetic anisotropies bi2se3 () or many body resonances hui () as soon as solid links between model and the realistic structure can be established.

In this letter, we have adapted a hybrid approach dft++ () (DFT++), which links Density Functional Theory (DFT) and Anderson’s impurity model to study physical properties of FeP and FePc. We demonstrate how the interplay of Coulomb interactions, crystal fields and hybridization with the ligands, which are fully captured in our theory, lead to correlated electron physics, and how this theory describes the spin-crossover in the Fe metal center of FeP. Furthermore, the crossover between different close lying ground state electronic configurations within the subspace of the correlated Fe-d orbitals is analyzed.

The problem can be cast in the following way. The delocalized orbitals in the organic ring in FeP are described by LDA or GGA but the Fe center is considered as an impurity embedded in the organic host and is described by Anderson’s impurity model anderson ().

The model Hamiltonian for the impurity problem can be expressed as :

(1) |

where describes the onsite energies and represents combined orbital and spin indices. is the annihilation operator while represents the local Coulomb interaction. is parametrized by Slater parameters. We have chosen U=4 eV and J=1 eV for (3d-orbitals) in our calculations. In the above equation, the first two terms represent electrons in Fe 3d orbitals. The third term describes the interaction with the surrounding atoms while the fourth term is for the delocalized ligand electrons with energies . The hopping matrix element appears in the hybridization function, which is represented as :

(2) |

The energy dependent hybridization function can be obtained from first principles calculations. We follow the approach considered by Karolak et al.dft++ () to construct the local Green function from DFT. The Kohn-Sham Green function can be calculated from the Lehmann representation economou () using

(3) |

where ’s and ’s are the Kohn-Sham eigenstates and eigenvalues for band and reciprocal space point . The projection of the full Green’s function to an atom centred local Green’s function for localized orbitals is needed, which in our case are cubic harmonics ().

(4) |

where and . The hybridization function is calculated from the local impurity Green’s function from the expression :

(5) |

In the above expression, is the projected Green’s function on local orbitals and combines the hybridization function and static crystal field . If the bath orbitals (defined by and in Eqn. (2)) are limited to a small number of discrete orbitals only, the many body problem defined in Eqn. (1) can be solved by means of exact diagonalization (ED), which will be used here.

Electrons in Fe-3d orbitals hybridize mostly with the orbitals of the four surrounding N-atoms, that provide a square planar ligand field. The effect of the outer C-ring is rather indirect as that mainly shifts N-p levels by both in-plane and interaction. The energy dependence of that interaction with surrounding N-atoms can be described by an energy-dependent hybridization function . We have employed non-spin polarized density functional calculations within local density and generalized gradient approximations to extract the hybridization functions. The DFT calculations were performed using the VASP code vasp () that employs a plane wave basis and projector augmented wave method. In Fig. 1, real and imaginary parts of are shown. The imaginary part of Im quantifies the density of bath states coupling to each impurity orbital weighted by the hybridization matrix elements V. As seen from Fig. 1, the most dominant peak in Im is observed for the Fe-d orbital at 2.04 eV below the Fermi energy. The formation of in-plane bonds with axial N-ligands explains this pronounced peak. It should be noted that the other in-plane orbital d shows almost no hybridization apart from a small peak at 4.8 eV below the Fermi energy. The out-of-plane orbitals have relatively small peaks in Im, among which the one closest to the Fermi energy is of d character at -2 eV. Appearance of this peak reflects a interaction of Fe d (d, d) orbitals with N-p orbitals, which is expected in the square planar ligand field of the FeP molecule. The other out-of-plane and in-plane contributions are present in the -4.5 to -10 eV energy range. Taken together, the hybridization function reveals a strong in-plane interaction between Fe-d orbitals and the N-p ligand states along with a much weaker interaction amongst all other orbitals.

The real part of the hybridization function Re describes the energy dependent ligand field, combined with a static crystal field, which can be obtained at the limit. The strong resonance with the host and d orbital manifests itself also around -2.04 eV in Re. This strong ligand field pushes d high in energy causing almost no occupancy in either spin channel. In the gas phase, six electrons in Fe, thus, are distributed in the remaining four d-orbitals, with four and two electrons in majority and minority spin channels, respectively, giving rise to an intermediate spin state (S=1).

An intriguing phenomenon of strain induced spin state switching is observed theoretically in FeP sumantachemi (); sumantaprl (); liao (), which is relatively difficult to obtain for other molecules with other TM centers kahn (); gutlich () or different structures, e.g., transition metal phthalocyanines. Iron Phthalocyanine (FePc), for example, has an Fe center but along with a larger organic ring. This results in a stronger hybridization of Fe with the neighboring N atoms, which is reflected in the appearance of a d -peak in Im at -1.85 eV and a stronger hybridization, as shown in the Supplementary Information (SI). The shift of the bath energy to -1.85 eV for FePc compared to -2.04 eV in FeP is due to the presence of four extra N atoms connected to the pyrolle rings in FePc (See Fig. 2 in SI). These four N atoms are not directly bonded to Fe but are connected to N atoms in the Fe-N block via C atoms.

To quantify the conditions required for spin state switching of FeP, we have varied crystal field and hybridization strength within DFT++ method. The variation of these parameters mimics the strain effect on the molecule. Also as mentioned in the previous section we found that the crystal field splitting of d is larger compared to other orbitals and as will be discussed below, this splitting is responsible for the spin switching. For this part of the calculation we have kept a common reference level of the rest of the orbitals separated from d by V, which was varied. An independent variation of the hopping matrix elements, , is done along with the variation of V. One needs to keep in mind that for our model calculation, V and are independent, while in DFT calculation (discussed in the latter part of this section), these parameters are implicitly related. The ligand field variation due to the strain effect simultaneously changes V and . For the remaining part of the discussion we will refer to as V, as this the only bath-site coupling considered for our model calculations.

Fig. 2 depicts the spin phase diagram of the central atom in FeP. The phases are defined by the characteristic energy contributions and demonstrated with RGB color code (see supplementary information). For V 0 and with sufficiently high crystal field (V 2.6 eV), 6 electrons occupy the degenerate d level. Hence, the ground state will be a low spin state with an energy gain from crystal field (red) but with the cost of exchange energy (blue). The ground state attains a high spin state for low crystal fields. A strong V pushes the molecular orbital containing predominantly d and hence, required crystal field for spin crossover is small. The intermediate spin state, in this situation is governed by the Fe-N bonding (green). Crystal field can be tuned to occupy anti-bonding molecular orbital and hence, S=2 spin state can be achieved (AS). The dependence of V with V is follows a quadratic behaviour. The fitted curve in Fig. 2 at the phase boundary is obtained from a mean field model (see supplementary information) containing one particle Hamiltonian but with a renormalized onsite energy . The discrepancy between the phase boundary and the fitted curve lies in proper renormalization of , that depends linearly with , by electron correlation energy, which has a dependence of .

As the spectrum of the FeP molecule is gapped, we placed chemical potential in the middle of the gap. But it is also possible to vary the average onsite d-energy , while keeping the total number of particles in the system fixed. In addition, the average -electron on-site energy is not exactly known from the DFT simulations — a problem usually referred to as ”double counting problem” in DFT++, DFT+U, or DFT+DMFT approaches. In our case, this means that the average on-site energy or more precisely the average energy difference between the bath level and the Fe -block carries some uncertainty. We thus vary in a range, while keeping the total number of particles in our system constant. In this way, the high-spin to low-spin transition line has error bars associated with some of the calculated points shown in Fig. 2.

At a stronger coupling V, the shifting of the onsite energy will allow a variation of V, where a spin crossover (SCO) can happen. The range of V due to the variation of is presented by the red error bars in Fig. 2. In the regime of weak hybridization, the system is described predominantly by crystal fields. There are in particular, no charge fluctuations to or Fe impurity states, which are generally affected by changing . As there is no mixing with impurity states in the limit of , the error bar is vanishingly small in this limit. For V beyond 2.8 eV, this particular scenario of crystal field will not be able to switch the spin state. An orbital reversal, i.e., the d energy becoming lowered compared to other orbitals is needed in this case, which requires a different kind of charge distribution in the molecule. Thus a transition between S=1 and S=2 states could be realized for V 2.8 eV while S=1 should be obtained generally for V 2.8 eV.

Weakening of the ligand field leading to a spin state change, however, needs a nontrivial chemical or physical procedure. The existence of high spin (S=2) porphyrin complexes with configuration has so far been observed in non-planar molecules with five or six coordination of the central Fe atoms2 (). Higher coordination leads to out of plane shift of the central Fe atom (five coordination) or symmetrical Fe-N block expansion, resulting in a weaker ligand field. For four coordination, however, the S=2 state is yet to be observed experimentally in gas phase or on a surface. Fig. 2 establishes the parameter space for when this state is to be expected.

The variation of V is also studied with PBE by varying the Fe-N bond length in the molecule. The highest value of V (3.16 eV) is obtained for a Fe-N bond length of 1.97 Å. The four data points, shown in Fig. 2 in filled orange circles, are for bond lengths of 2.00, 2.04, 2.07 and 2.11 Å, respectively. As mentioned in Refs. sumantaprl, ; sumantachemi, for FeP either physisorbed or chemisorbed on surfaces, the required bond length of Fe-N in FeP for spin switching is beyond 2.03 Å within PBE+U approximation with U=3 eV. marzari () This is in agreement with the data shown in Fig. 2, where both DFT++ and PBE suggest a spin crossover beyond 2.04 Å. A detailed comparison between the bond lengths required for spin-crossover in PBE+U and DFT++ is shown in Fig. 3. A comparison of PBE and LDA calculations yield values of the static crystal field, bath energy and hybridisation as 0.62(0.65) eV, -2.04(-1.9) eV and 3.16(3.39) eV for PBE(LDA). These values indicate that the two approximations for exchange-correlation functional do not show pronounced difference.

In Fig. 3, calculated magnetic anisotropy energies (MAE) within DFT++ are shown for different strengths of V. The MAE is defined as the energy difference between the lowest many-body eigenstate with out-of-plane magnetic moment and the ground state which has in-plane oriented moment. It is observed that for low values of V (large Fe-N bond lengths), the MAE decreases with decreasing V whereas for large V, a more or less constant value of MAE is obtained. In all the cases, we have found the easy axis of magnetization to lie in the plane of the molecule. In the inset of Fig. 3, the spin-crossover properties of FeP as predicted by different flavors of PBE+U and DFT++ is compared. In the PBE+U calculations, we stabilize the non-favorable spin solutions by constraining the spin moments. It is seen from the figure that there is a notable difference between two PBE+U methods (Dudarev and Lichtenstein) in the Fe-N bond length required for LS-HS transition, while the results obtained from Lichtenstein PBE+U and DFT++ methods are close to each other. This is understandable since both, the Lichtenstein variant of PBE+U and DFT++, employ the full four fermion Coulomb matrix as defined through the Slater parameters, whereas Dudarev assumes a simplified Coulomb vertex.

The free molecule spin state (S=1) is particularly interesting from the point of view of ground state configuration. A strong ligand field in the free molecule leaves d nearly unoccupied and an intermediate spin state (S=1) can have ground state electronic configuration with 2 electrons in d, 3 electrons in d and 1 electron in d. We will refer to this configuration a (ddd) configuration. The other possible configurations within the S=1 ground state multiplet are (ddd), (ddd) and (ddd) among which the appears to be very close in energy to the energy of . From Re, a splitting can also be seen among , d and d orbitals but in a relatively small energy scale. To acquire a clear view of the ground state configuration, we varied the crystal field in presence and absence of coupling V. At the first step we only considered a crystal field splitting between d and the averaged position of remaining orbitals, as shown in the left most part of Fig. 4. In this pure crystal field situation, the S=1 state leads to the ground state configuration. As the coupling V is switched on, the ground state configuration becomes onsite energy dependent. This happens because the inclusion of strong hybridization will admix N=5 and N=7 impurity occupancies with the pure N=6 ground state configuration of the pure crystal field situation. Varying on site energy, the ground state configuration can be modified. The energy scale associated to this change is of the order of meV. A relative splitting among the remaining orbitals has even more pronounced effects in determining the ground state configuration. As shown in Fig. 4, with additional splitting, if d stays above d, is stabilized. However, in the reversed situation, is obtained. This change occurs in at least one or two order higher energy scale compared to the many-body effects induced configuration change, revealing the crystal field to be the dominant factor.

In summary, we have presented ground state electronic properties of Fe in FeP molecule with a hybrid approach of DFT combined with a many body treatment, using exact diagonalization. We have demonstrated that a delicate interplay of the static crystal field, ligand hybridization and Coulomb interactions promotes iron porphyrin to be a potential candidate for realizing spin-crossover behavior. In general, our calculated phase diagram indicates the possibility of tuning electronic and magnetic properties of organometallic molecules to serve the purpose of molecular electronics and storage devices. Moreover, the long-standing debate regarding the electronic ground state configuration, e.g., vs. has been solved by identifying the proper parameters required to switch one to the other, which will have important consequences for the spin dipole moments and magnetic anisotropies, where the energy positions of d-orbitals with specific symmetries are important.

B. S. acknowledges Carl Tryggers Stiftelse for financial support. O.E. Acknowledges support from VR, the KAW foundation, eSSENCE. We acknowledge supercomputing allocation from Swedish National Infrastructure for Computing. T.O.W. and M. S. acknowledge support from the Central Research Development Fund of the University of Bremen and the DFG via FOR 1346.

## References

- (1) H. Wende et al., Nat. Mater. 6, 516 (2007); V. A. Dediu, L. E. Hueso and C. Taliani, Nat. Mater. 8, 707 (2009).
- (2) S. Bhandary et al., Phys. Rev. B 88, 024401 (2013).
- (3) S. Bhandary et al., Phys. Rev. Lett. 107, 257202 (2011).
- (4) S. Bhandary, O. Eriksson and B. Sanyal, Nat. Sci. Rep. 3, 3405 (2013).
- (5) P. Fazekas, Lecture Notes on Electron Correlation and Magnetism, Series in Modern Condensed Matter Physics- Vol.5, World Scientific (1999).
- (6) J. Honolka et al. Phys. Rev. Lett. 108, 256811 (2012)
- (7) H. Zhang et al. Phys. Rev. Lett. 107, 026801(2011)
- (8) P. W. Anderson, Phys. rev., 124 41 (1961).
- (9) M-S. Liao and S. Scheiner, J. Chem. Phys. 116, 3635 (2002); W. R. Scheidt and C. A. Reed, Chem. Rev. 81, 543 (1981); D. C. Rawlings, M. Gouterman, E. R. Davidson and D. Feller, Int. J. Quant. Chem. xxviii, 773 (1985).
- (10) M. Karolak, T. O. Wehling, F. Lechermann and A. I. Lichtenstein, J. Phys. Condens. Mat. 23, 085601 (2011).
- (11) E. N. Economou, Green’s Functions in Quantum Physics, Springer Series in Solid-State Sciences (2006).
- (12) G. Kresse G. and J. Hafner, Phys. Rev. B 47, R558 (1993); G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
- (13) O. Kahn and C. J. Martinez, Science, 279 44-48 (1998).
- (14) P. Gütlich and H. A. Goodwin, Spin Crossover in Transition Metal Compounds I, 233 1-47 (2004).
- (15) C. A. Reed et al., J. Am. Chem. Soc., 102, 2302 (1980); J. P. Collman et al., J. Am. Chem. Soc., 97, 2676 (1975); G. B. Jameson et al., J. Am. Chem. Soc., 102, 3224 (1980)
- (16) D. A. Scherlis, M. Cococcioni, P. Sit and N. Marzari, J. Phys. Chem. B 111, 7384 (2007).