Many-Body Effects in FeN 4 Center Embedded in Graphene

: We introduce a computational approach to study porphyrin-like transition metal complexes, bridging density functional theory and exact many-body techniques, such as the density matrix renormalization group (DMRG). We ﬁrst derive a multi-orbital Anderson impurity Hamiltonian starting from ﬁrst principles considerations that qualitatively reproduce generalized gradient approximation (GGA)+U results when ignoring inter-orbital Coulomb repulsion U (cid:48) and Hund exchange J . An exact canonical transformation is used to reduce the dimensionality of the problem and make it amenable to DMRG calculations, including all many-body terms (both intra- and inter-orbital), which are treated in a numerically exact way. We apply this technique to FeN 4 centers in graphene and show that the inclusion of these terms has dramatic effects: as the iron orbitals become single occupied due to the Coulomb repulsion, the inter-orbital interaction further reduces the occupation, yielding a non-monotonic behavior of the magnetic moment as a function of the interactions, with maximum polarization only in a small window at intermediate values of the parameters. Furthermore, U (cid:48) changes the relative position of the peaks in the density of states, particularly on the iron d z 2 orbital, which is expected to affect the binding of ligands greatly.

Despite the relative simplicity of graphene and a transition-metal atom, understanding the combined electronic structure of a transition-metal center embedded in the graphene matrix remains a challenge as the transition metal atom typically involves an incomplete d shell that gives rise to confinement-induced correlations and magnetism. Studies of iron porphyrins and other heme-like molecules have been carried out within the framework of density functional theory (DFT) [44,[50][51][52][53][54][55][56][57], DFT+DMFT (dynamical mean-field theory) [35,58,59], quantum Monte Carlo [60,61], coupled-cluster [51], molecular dynamics [62], and configuration interaction techniques [63]. Here, we discuss a new approach in which we use the DFT to construct a material-specific effective model Hamiltonian, which is then solved via the density matrix renormalization group (DMRG) [64][65][66][67][68] technique along with a unitary transformation to account for the many-body physics in a numerically exact way. Due to the large number of geometries, transition-metal centers, and axial ligands involved, a versatile method such as that discussed here to account for strong correlation effects will be useful for exploring the properties and possible applications of these complexes.
Finding accurate approximations to the exchange-correlation functional represents the central problem of DFT. The simplest such functional is the local density approximation (LDA) [69]. LDA generally yields reasonable ground-state properties of 3d magnetic metals, including magnetic moments and Fermi surfaces [70]. LDA however does not correctly capture the phase diagrams of magnetic materials, including in particular the lowest energy crystal structure of elemental iron, where it predicts non-magnetic iron to be face-centered cubic and ferromagnetic iron to be body-centered cubic. The generalized gradient approximation (GGA) [71] provides a systematic improvement over the LDA with the added variational freedom by including the density gradient in the construction of the exchange-correlation functional. The GGA correctly predicts the relative stability of the ferromagnetic phase of pure Fe in addition to giving a good description of its ground state properties [72]. Further studies of iron porphyrin-type molecules based on the GGA+U scheme have revealed that interaction effects on the iron play an important role in determining the ground state magnetic moment [52,58,73].
In this paper, we consider a FeN 4 center embedded in a graphene lattice. For this purpose, we take a FeN 4 C 10 H 10 complex as the fundamental building block for this center, which was referred to as the D1 center in [47]. D1 centers have been realized experimentally in graphene and carbon nanotubes [40,44]. We adopt this configuration since it requires minimum structural modification of the graphene backbone. While GGA calculations with U = 0 capture the magnetic ground state quite well [55], U is needed to describe the energy splittings in the spectrum. Using a more physical value for the effective on-site Coulomb repulsion U of about 4 eV and the Hund interaction strength J of about 1 eV, the spin state changes to a high-spin state, as expected. However, these calculations still neglect the inter-orbital Coulomb effects, which will be shown below to be important.
The outline of the paper is as follows. The five-orbital Kanamori-Anderson effective Hamiltonian is derived in Section 2. The method used to solve our model Hamiltonian is discussed in Section 3. The solution for the molecule and graphene and the corresponding phase diagrams are described in Section 4. The concluding section provides a summary and discussion of our main results.

Model Hamiltonian
In order to derive a simplified model Hamiltonian, we focus on taking into account only the most relevant features that determine the electronic structure and occupation of the transition metal atom. Correlation effects in the carbon and nitrogen atoms are neglected, and the organic backbone of the molecule is constructed by using an LCAO-based tight-binding Hamiltonian. Exact numerical methods are used to account for the many-body physics introduced by the central atom, which we model as a multi-orbital Anderson-like impurity. Generalized Anderson impurity models have been applied previously to porphyrin-like molecules [32,34,37,38,74,75] to predict potential energy surfaces and electronic coupling factors [63,76] for various transition-metal complexes. The advantages of our model Hamiltonian approach are that in this way, we can account for the many-body correlation effects in a numerically exact way, and also, the computations can be scaled relatively easily to handle systems with multiple impurities and more complex geometries.
The transition-metal center is created by removing six carbon atoms from the graphene sheet and replacing these with four nitrogens and a single iron atom, as shown in Figure 1a. For simplicity, we assume the structure to be planar with point-group symmetry D 2h . In order to identify the active orbitals, we start by considering the carbon atoms to build a model for graphene. The carbon atoms are connected by σ-bonds of hybrid sp 2 orbitals, formed from linear combinations of carbon 2p x , 2p y , and 2s orbitals. The π-bonds are weaker than the σ-bonds, formed by the remaining p z orbitals. The graphene sheet is then functionalized by adding hydrogen atoms bonded with each dangling sp 2 orbital of the carbons at the edges. We approximated the values of the hopping parameters to approximate the results for bulk graphene. Modeling the p z orbitals then only requires a simple nearest-neighbor hybridization t, resulting in the well-known two-band model of graphene. Graphene's σ bands require a more sophisticated approach; see Appendix A for details. As an approximation, the N-C hopping is set equal to the C-C hopping. The nitrogens however have a remaining orbital that points toward iron.
Looking at the coordination of iron and nitrogen atoms, the complex is approximately square planar with D 4h symmetry, which is useful in determining the bonding orbitals and the crystal-field splitting. Using only the 3d-orbitals of iron, d x 2 −y 2 will form a σ bond with the dangling sp 2 orbital from the nitrogen, while π-bonds will form from the iron d xz , d yz and nitrogen 2p z orbitals, where we define the x and y axes to point along the lines connecting the Fe and the N atoms. The five dimensional (l = 2) irreducible representation of K (the continuous rotation group) becomes reduced to four irreducible representations when the symmetry is decreased to D 4h . If only electrostatic effects are considered, the crystal field splitting is predicted to be as seen in Figure 2a. However, this does not capture higher order effects or Jahn-Teller distortions [77]. For a more accurate description, the levels in a graphitic structure are arranged as in [47], as shown in Figure 2b. As discussed below, the positions of the energy levels are adjusted to match the occupations of the 3d orbitals approximately given by the DFT at U = 4.0eV. They are then held fixed in all other calculations. (b) Arrangement of energy levels used throughout this work following a similar convention adopted by Aoyama et al. [47]. The difference between the π and z 2 orbitals is small and has been reversed to match the occupation of the orbitals according to DFT calculations.
For iron, we include Coulomb and Hund interactions. The resulting general form of the Hamiltonian can be written as [78,79]: In the preceding expression, m labels the d-orbitals and V (r) is the screened Coulomb potential. This can be simplified to four matrix elements known as Kanamori parameters [80]: the intra-band Coulomb interaction U, the inter-band Coulomb interaction U , the inter-band exchange interaction J, and the pair-hopping amplitude J (Notably, in quantum chemistry codes such as Gaussian, Jaguar, and DMol, Coulomb integrals are denoted by the symbol J and not U. The symbol J is used here for magnetic exchange or Hund couplings, following common practice in the physics literature.). It can be shown [81] that J = J due to the symmetry of the orbitals involved and the fact that all coefficients are just integrals of the Coulomb term over the radial part of the wave functions. In order to ensure rotational invariance in orbital space, the condition U = U + 2J must be satisfied [78] for all combinations of orbitals. The effects of the crystal field on the Coulomb interactions are ignored and are assumed to be relatively small. The interaction part of the Hamiltonian now takes the form: The interactions U, U , and J can be expressed in terms of the Racah parameters A, B, and C, as shown in Table 1. The values of B and C were given in [82] for Fe 3+ and Fe 2+ and are restated in Table 2. While U and J take different values for different pairs of orbitals, the interaction U = A + 4B + 3C is orbital independent.  Our calculations assumed that there were approximately six electrons in the iron orbitals so that we were close to Fe 2+ . This leaves only one free parameter for the interactions, i.e., A, or equivalently U, since all U and U have the same dependence on A (note that J is independent of A). In general, it is the competition between U and J that determines if the complex is high or low spin.
In addition, spin-orbit coupling can be introduced as it becomes important for heavier elements that could substitute the iron in similar molecules. However, this interaction for iron is relatively small (≈0.05 eV) [83], and we have neglected it in our calculations, although it has been shown to be relevant in iron-based superconductors [84]. Notably, multi-orbital models for transition metals have also been used in dynamical-mean-field calculations of (bulk) correlated materials to discuss orbital-selective Mott transitions in iron-based superconductors, where the results can be sensitive to crystal-field splitting [35,[85][86][87]. Our focus in this study, however, is on understanding the correlation effects on a material-specific model Hamiltonian for the FeN4 cluster in graphene.

Method
Once all the parameters for the model are obtained, the problem can be recast and solved using the DMRG method. In order to do so, we mapped the problem onto an equivalent one-dimensional model by employing an exact canonical transformation, as presented in [88,89] and reviewed in detail in [90] and closely related to the chain mapping used in DMFT [35]. An outline of the method is as follows.
A conventional Hamiltonian for impurity problems will have the form: Here, H l denotes the single-particle tight-binding Hamiltonian for the lattice, which could include more than one band, and it can be obtained from DFT calculations. H imp and V describe the impurity and the coupling between the impurity and the lattice, respectively. Note that this method is applicable regardless of the geometry or dimensionality of the non-interacting Hamiltonian. The key is to map H l onto an equivalent one-dimensional chain. First, for simplicity, let us consider a single-impurity problem with one orbital per site. The more general case involving multiple orbitals and impurities will be discussed in the next section. The first step is to define the "seed" state to perform a Lanczos recursion as: where c † r 0 creates an electron at orbital r 0 and |0 is the vacuum state. For an Anderson-like impurity, which is the case here, we chose the seed to be the impurity orbital. The remaining states are then constructed with the following iterative procedure: The equations for a n and b n can be obtained by requiring the states to be orthogonal. Note, however, that at this stage, the states are not normalized.
After this transformation, H l assumes a tri-diagonal form: which corresponds to a one-dimensional Huckel Hamiltonian. Equivalently, in second quantization, it reads: wherec † i andc i are the normalized creation and destruction operators, respectively,ñ i =c † ic i is the particle number operator, and L is the total length of the chain. The diagonal a n terms are on-site potentials, while the b n 's are the new hopping along the chain. This is an exact canonical transformation. The remaining missing orbitals correspond to different symmetry sectors of the Hamiltonian, are completely decoupled from the impurity, and can be safely ignored, which highlights the power of this change of basis.
Coming back to the transition-metal complex, in our case, we had two orbitals (iron d xy and d yz ) coupled to different sites (nitrogens) that will generate two orthogonal chains along the lines very similar to those described above, except that we would now require two seeds for the π-bonded nitrogen p z orbitals. Labeling these as |α 0 and |β 0 , we chose these as: where the labeling refers to the nitrogen sites in Figure 1 (Note that due to symmetry, only these two -out of four-wave-functions couple to the d orbitals of the transition metal because the other degrees of freedom live in an orthogonal Hilbert space). After the Lanczos iterations are carried out, two chains represented by the red sites in Figure 3 are generated. The coupling Hamiltonian between the nitrogen and iron then becomes, For the σ-bands, the single seed mapping can be used starting from the d x 2 −y 2 orbital. The hopping integrals between the nitrogen and carbon 2s, 2p x , and 2p y orbitals are obtained from the tight-binding model described above. After this transformation, we obtain the green chain in Figure 3. As expected, due to the symmetries of the problem and the resulting dimensional reduction, the total number of orbitals in the equivalent system is smaller than that in the original system. The magnitude of the hopping between the iron and nitrogen sp 2 orbitals can be used as a fitting parameter, while the π-coupling between the nitrogen and iron can be estimated via comparison with DFT results to yield t ≈ 1.6 eV. In order to study the problem of iron embedded in bulk graphene, we extended the two sides of the chain in Figure 3 to the desired length as discussed in [88].

DFT is known to have difficulties in describing correlated electrons in open d and f shell systems.
In particular, DFT can fail in predicting if the ground state has low, intermediate, or high total spin polarization (this failure can however be mitigated if we allow spin contamination (spurious mixing of different spin-states) in the calculations [55]). Conventional DFT methods, such as LDA or GGA, also fail to account for the effects of the Coulomb interaction between localized electrons properly. For these reasons, the calculations presented here were performed using GGA+U within the VASP package [91,92]. One main source of error in the GGA+U scheme arises from the fact that the U and J terms in the Hamiltonian (or functional) are handled in a mean-field fashion, resulting in the single U e f f = U − J parameter. With this in mind, we first compared our exact DMRG results with DFT+U calculations, where we ignored all many-body terms except for the intra-orbital Coulomb repulsion and the Hund coupling of the spins. In other words, H int = U ∑ m n ↑m n ↓m − 2 ∑ m =m J mm S m S m .
In order to benchmark our approach, we began with the iron complex, FeC 10 N 4 , depicted in Figure 1. Here, we used level splittings between the d orbitals that matched the occupations obtained from GGA+U using the physically relevant value of U = 4 eV, as shown in Table 3 (Note that obtaining exact agreement between the DFT results and our model Hamiltonian was not possible due to the complex interplay between various terms in the Hamiltonian. The results shown were as close as we could get by searching within the multi-dimensional parameter space and fixing the value of U. Our formulation was SU(2)-invariant, so that any high-spin ground state would be a (2S + 1)-fold degenerate multiplet.). The ground-state occupation and magnetic moment of iron (total spin) as a function of U is shown in Figure 4. Calculations were carried out by varying the number of electrons and the value of the total spin S. The ground state was then obtained by minimizing the energy. The spin was seen to remain zero until the Coulomb interaction reached a value of U ≈ 3.0 eV. As the value of U increased beyond this point, the occupation of the iron levels decreased as the magnetic moment increased until it saturated at S = 2 ( N = 4). For large values of U, Coulomb repulsion prevented double-occupation of orbitals. Over the range of 3.5 ≤ U ≤ 4.5, the ground state has S = 1, which is the physically interesting range (Note that different Racah parameters should be used, in principle, for various occupations of the iron atom, for simplicity, we used a fixed value of these parameters, which corresponded to Fe 2+ ). When all the interaction terms were included in the Hamiltonian, the ground state changed drastically. In order to gain insight into this result, we considered a rigid shift, −V Fe ∑ m n m , in the position of the Fe energy levels. This ad hoc potential could be physically thought of as representing the screened interaction with the nucleus, so that it could be used as a parameter to mimic the occupation of the levels for various transition metals. In Figure 5, we show the overall occupation of the molecule, the Fe atom, and the total spin S, as a function of U and V Fe . The parameter regime of interest for Fe 2+ (S = 1) was seen to reside in a narrow band of values in the so-called C 231 electronic configuration (d 2 xy d 3 π d 1 z ) [93]. Other bands in the figure correspond to a different ionic state of iron or a different atomic species. Within each region, the occupation of various energy levels did not vary much. Figure 6, which focuses on the parameter regime corresponding to V Fe = 0, shows that the total charge and spin of the transition metal atom depended strongly on U. The physically-relevant region with S = 1 and N ∼ 6 occupied the small range 2.3 < U < 2.8. More importantly, increasing U further eventually led to a low-spin state (S = 1/2 or S = 0), as opposed to the high-spin state reached above without inter-Coulomb repulsion. The high-spin states in this model were thus limited to a finite window of U. This was due to the fact that as U was further increased, all orbitals became single-occupied, and the inter-orbital interactions began to dominate.  We turn now to discuss the effects of passivating our defect center via hydrogen atoms bonded to the dangling sp 2 bonds of the carbons. Calculations carried out with and without bonded hydrogens showed that there was not much overall change in the underlying physics, although occupation of the d x 2 −y 2 orbital was modified somewhat, and the effect was negligible when the U terms in the Hamiltonian were not included in the calculations. In all cases we considered, hydrogens did not significantly affect the overall spin state of iron. For these reasons, we only show the results without the hydrogens. Table 3. Occupations and magnetic moments of various iron orbitals. The first two columns are results from DFT+U calculations. The next three columns are obtained from the method described in the text with U = 0, while the last three columns refer to calculations in which all inter-orbital interactions are included.  Figures 7 and 8 show the projected density of states for the five iron orbitals calculated with dynamical DMRG [94,95] and DFT+U, respectively. The top panel in Figure 7 (without inter-orbital Coulomb interactions) shows qualitative agreement with DFT+U calculations (Figure 8) in the xy, z 2 , and x 2 − y 2 channels. States in the π channel in Figure 7, however, are more spread out than in Figure 8 due to greater hybridization with the conduction electrons. Note that the states associated with the z 2 and xy orbitals were sharp and did not experience substantial charge fluctuations, and for this reason, these states were reasonably described within the mean-field treatment, while many-body effects were reflected mainly in the π channel. Similar effects were also seen in the exact diagonalization calculations on a simplified five-orbital model [35] of a flat band with constant hybridization (our DMRG calculations reproduced the exact diagonalization results, but were not limited to small clusters). Parameters such as U, the energy levels of the d-orbitals, and the couplings were adjusted to obtain agreement with the occupation and magnetic moment given by DFT, since we expected DFT to be a reasonable approximation when inter-orbital interactions were neglected. The differences between DMRG and DFT reflected the effects of: (i) approximations inherent within the DFT, (ii) the parameters used in our DMRG modeling based on the DFT results, and (iii) the fact that we ignored Coulomb interactions between iron and nitrogen, as well as other parts of the defect complex even though these interactions were to some extent accounted for by the effective hopping parameters.

DFT+U
As expected, the inclusion of interactions led to modifications in the partial density of states. It was clear that the U terms shifted the energies slightly upward, with the exception of the d xy orbital, which was practically unchanged. Note however that the value of U was reduced substantially, from U = 4.0 to U = 2.6, indicating the importance of interaction terms. The d z 2 orbital was also affected greatly as seen from the splitting between the related peaks in Figure 7, which was controlled by U, while its relative position with respect to other orbitals was dictated by U and J. These results indicated that interaction effects would induce substantial changes in the binding energies of ligands, which normally involved the d z 2 orbital.
We also extended our modeling to include more carbon atoms in order to gain insight into the effects of embedding our defect center in a "bulk-like" flake of graphene. Here, our calculations included all interaction terms. As shown in Figure 9, for certain values of U, the spin actually increased compared to the molecular case. An intermediate S = 3/2 phase, which was not found before, appeared around U ≈ 2.6 eV, indicating that the surrounding material played a significant role in the physics of the defect center: a continuous density of states in bulk graphene compared to discrete "delta"-like peaks in the molecular case may allow for additional screening in the spirit of the Kondo effect.

Conclusions
We investigated heme-like iron centers in graphene and the FeC 10 N 4 molecule using an exact canonical transformation within the framework of the DMRG method. The DMRG technique has been used in quantum chemistry calculations as a solver for first-principles Hamiltonians in the same spirit as configuration interaction calculations [96][97][98][99]. Our approach took advantage of the weakly correlated nature of the carbon bond, which could be accounted for by DFT calculations, and recast the problem onto an LCAO model with an interacting transition metal center that was modeled as a five-orbital Kanamori-Anderson impurity. This allowed us to perform a unitary transformation that significantly simplified the Hamiltonian by reducing it to an equivalent Hamiltonian consisting of one-dimensional chains coupled to the iron d-orbitals. This made the problem amenable to efficient DMRG calculations in which all many-body terms were treated exactly. In this way, we obtained the occupation and magnetic moment of the iron atom as a function of the strength of the Coulomb interaction U and qualitatively recovered DFT results in the limit where the inter-orbital repulsion was neglected. Interaction effects were shown to play a crucial role in the physics underlying the defect center by shifting the relative positions of peaks in the density of states associated with various iron orbitals. These shifts could be expected to drive substantial changes in the binding energies of various ligands to the defect cluster. Moreover, our results demonstrated the importance of orbital-dependent interactions to account for correlation effects in an accurate manner, and thus question the common assumption that values of U ≈ 4 eV and J ≈ 1eV are sufficient to characterize the physics of the transition metal atom.
The technique described in this study can be combined with other quantum chemistry approaches such as CASPT2 [100,101], not only for benchmarking, but also to obtain realistic parameters to model bulk-embedded transition-metal complexes by mapping the system onto one-dimensional chains. Natural extensions of the method would include treatment of a variety of geometries of transition-metal complexes and the effects of spin-orbit interactions and correlated hybridizations [102]. Moreover, one could investigate two iron atoms embedded in a sheet of graphene to explore the possible emergence of indirect magnetic exchange interactions mediated by the conduction electrons. Finally, our approach can be extended readily to handle finite temperatures and adapted to study non-equilibrium phenomena such as transport and chemical reactions.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Tight-Binding Model for Graphene's σ-Bands
Before discussing the electronic structure of the σ bands, we briefly present the general approach for determining the vanishing and non-vanishing matrix elements between various orbitals. For the s and p orbitals, there are just four non-zero overlap integrals to consider: ssσ, spσ, ppσ, and ppπ. Due to the radial symmetry of the s-orbitals, the ssσ bond has no angular dependence. σ and π bonds are classified by whether the interatomic separation direction is parallel or perpendicular to the orbital axis. Since the p-orbitals can orient with any angle between them ( Figure A1), we project these orbitals into their normal and parallel (σ and π) components as follows. |p = cosθ |p σ + sinθ |p π .
Matrix elements between neighboring s and p states can then be written as: where the definition s| H |p σ = H spσ has been used, and s| H |p π = 0 by symmetry. The angle θ is defined in Figure A1. Similarly, the matrix elements between p states are: p 1 | H |p 2 = H ppσ cosθ 1 cosθ 2 + H ppπ sinθ 1 sinθ 2 . In general, sp 2 orbitals are planar and form angles of 120 deg. Since the unit cell of graphene has two atoms and each atom contributes three sp 2 states, we obtain six σ-bands, which are split into two groups of three bands each where one group lies above and the other below the Fermi energy. Following [103], the 6 × 6 matrix takes the following form: with all the remaining elements being zero. The numerical values for the hopping are (in eV): s = −8.7, p = 0, H ssσ = −6.7, H spσ = 5.5, H ppσ = 5.1, and H ppπ = −3.1.