Magnetic and Transport Properties of Fe4 Single-Molecule Magnets: A Theoretical Insight

Here, methods of density functional theory (DFT) were employed to study the magnetic and transport properties of a star-shaped single-molecule magnet Fe4 S=5 complex deposited on a gold surface. The study devoted to the magnetic properties focused on changes in the exchange coupling constants and magnetic anisotropy (zero-field splitting parameters) of the isolated and deposited molecules. Molecule-surface interactions induced significant changes in the antiferromagnetic exchange coupling constants because these depend closely on the geometry of the metal complex. Meanwhile, the magnetic anisotropy remained almost constant. Transport properties were analysed using two different approaches. First, we studied the change in the magnetic anisotropy by reducing and oxidizing the Fe4 complex as in a Coulomb blockade mechanism. Then we studied the coherent tunnelling using DFT methods combined with Green functions. Spin filter behaviour was found because of the different numbers of alpha and beta electrons, due to the S=5 ground state.


Introduction
The field of molecular electronics has been developed with the goal of providing miniaturized electronic devices beyond the limits of silicon technology. It has been suggested that the new technologies "More than Moore" and "beyond CMOS" 1 could achieve this in the future and it is clear that heterogeneous devices could replace some silicon components. The use of molecules in such heterogeneous devices could open up new functionalities, especially if magnetic molecules are employed. 2,3 Initially, the combination of molecules and metal surfaces was considered a drawback because in some cases, such as the functionalized single-molecule magnet (SMM) consisting of a Mn 12 molecule, deposition on the surface was followed by the loss of SMM behaviour due to electron transfer from the surface. 4,5 SMMs are compounds with a relatively large magnetic anisotropy that allows us to fix the spin direction even in the absence of an external magnetic field. 6,7 Thus, they have been proposed for the storage of information at the molecular level. However, recently, in the field of Molecular Spintronics, a new challenge has arisen: to control the electronic structure of the "spinterface" that can provide the combination of molecule+surface system with unexpected new magnetic properties. 2,8 Among others, it is worth to mentioning the magnetic properties of non-magnetic copper or manganese surfaces covered by fullerene molecules under a threshold magnetic field. 9 Also, noteworthy is the spin polarization of the current when injecting electrons from a gold surface to deposited chiral molecules [10][11][12] or small magnetic molecules 13 at room temperature in single-molecule devices. If the goal is to make devices using molecules showing magnetic anisotropy into SMMs, 14 one of the logical options is to employ molecules that are robust from a chemical point of view; 15 that is, molecules that are more electrochemically stable than the mixed-valence Mn 12 systems. The two most commonly employed systems are Fe 4 complexes [15][16][17][18] and Tb III phthalocyaninate. [19][20][21][22] The use of such molecules is relatively new and most of the experiments have been performed far from room conditions; such as in ultrahigh vacuum and at low temperatures. The use of magnetic molecules opens up the possibility of creating new devices for emerging technologies such as those based on spin-transfer torque mechanisms, 23,24 where the key property is the exchange interaction between the carriers of a spin-polarized current and the spin of the molecule.
Sessoli and coworkers have performed experiments with Fe 4 complexes grafted onto gold surfaces. 15,16,[25][26][27][28][29][30] The magnetic properties of the Fe 4 complexes can be summarized as an S=5 ground state, due to antiferromagnetic coupling between central and terminal Fe III centres, and also as exhibiting small magnetic anisotropy, due to the isotropic d 5 electron configuration of the Fe III cations. 31 Such deposited molecules behave as SMMs, showing hysteresis loops at low temperatures in X-ray Magnetic Circular Dichroism measurements. 15,16,26 More recently, they studied the influence of scanning tunnelling microscopy electrodes on the exchange interactions and found considerable strengthening when the Fe 4 complex is placed in the two-contact junction. 28 The same type of complexes have been studied by the van der Zant and Cornia groups using junction devices to analyse the influence of the electric field on their magnetic properties 18,32,33 and the role of vibron-electron coupling in their transport properties. 17 The analysis of the anisotropy of a single Fe 4 captured between the electrodes indicated that its transverse anisotropy is larger than the bulk-phase value. 34 Our goal here is to provide a complete theoretical study of the electronic structure, and magnetic and transport properties of SMM Fe 4 complexed on surfaces in order to rationalize our understanding of the properties of single-molecule devices based on such systems (see Fig. 1). The magnitude of the magnetic anisotropy is the key physical parameter for an SMM system, because it determines the height of the spin-flip barrier D·S 2 , with D being the axial zero-field splitting parameter and S the total spin of the ground state, which depends on polynuclear complexes of the ferromagnetic or antiferromagnetic couplings between the paramagnetic centres. Thus, our study will initially focus on the magnetic properties (exchange interactions and zero-field splitting) of the isolated molecules. The next step will consist of analysis of how such magnetic properties change when the Fe 4 molecules are deposited on gold surfaces. Finally, we will analyse the redox processes involved in a Coulomb blockade regime and their influence on the magnetic anisotropy by considering just the isolated Fe 4 molecule and the transport properties of the coherent tunnelling regime of single-molecule devices with the Fe 4 complexes sandwiched between two gold electrodes (see Fig. 1).

Computational details
Exchange coupling calculations for isolated and deposited Fe 4 complexes were carried out using the SIESTA (Spanish Initiative for Electronic Simulations with Thousands of Atoms) code. 35,36 The code is highly efficient when dealing with systems that contain a large number of atoms and also periodic structures, and therefore for those studied here. The generalized-gradient approximation (GGA) functional expression of Perdew, Burke and Ernzerhof (PBE) 37 was employed and pseudopotentials were generated according to the method suggested by Trouiller and Martins. 38 For the interaction of the Fe 4 molecule with terminal S· radical groups, an S=11 state was considered; due to the long distance between the central Fe 4 and the radical, the ferromagnetic S=11 and antiferromagnetic S=9 cases are degenerate. For the Fe 4 complex with no ligand functionalization to covalently anchor the gold surface, the DRSLL functional including a dispersion term was employed to optimize the molecule-surface interaction. 39 A triple-ζ numerical basis set with polarization functions for the metal atoms was employed, while a double-ζ basis with polarization functions was used for the other elements. Values of 50 meV for the energy shift and 250 Ry for the mesh cut-off provide a good compromise between accuracy and the computational cost of estimating the exchange coupling constants, according to previous studies 40, 41 A detailed description of the procedure used to calculate the exchange coupling constants can be found in previous papers. [42][43][44][45] A phenomenological Heisenberg Hamiltonian was used, excluding the terms relating to magnetic anisotropy (D and E zero-field splitting parameters) to describe the exchange coupling in the polynuclear complex: where Ŝ a and Ŝ b are the spin operators of the different paramagnetic centres, and the J ab parameters are the pairwise coupling constants between the paramagnetic centres of the molecule.
In general terms, we needed to calculate the energy of n+1 spin distributions for a system with n different exchange coupling constants. In our particular case, only first neighbour exchange interactions of the Fe 4 complexes were considered, thus, depending on the symmetry, there are complexes with from only one to three J values. Thus, for each system, the high spin S=10 spin configuration was calculated and the spin configurations corresponding to the spin-flipping of each terminal Fe III centre directly provided the J value between the terminal and the central Fe III centres.
Common exchange-correlations functionals provide good agreement with the experimental data avoiding the double-counting problem of some electronic correlation contribution generated by the spin-projection method. Due to the large spin of the Fe III centers, the difference in the calculated J value for the inclusion of the spin-projection is only 20%.
The calculation of the zero-field splitting parameters using DFT methods for molecular systems was introduced by Pederson and Khanna. 46 The method introduces the spin-orbit effect through second-order perturbation theory with the following spin-orbit operator: where S is the spin moment operator, L is the angular moment operator, r is the distance and F(r) the Coulomb potential operator. V x matrix elements are then defined as: and second-order perturbative energy can be expressed as: where σ considers all the spins, while i and j correspond to the three spatial directions. The matrix elements are defined as: for the full and empty orbitals, respectively, with e s and e s the corresponding energies. Due to the similarity of Eq. 4 with the spin Hamiltonian (Eq. 6), it is possible to compute such a zero-field D value: Thus, for a diagonal form of the tensor, the following expression can be obtained: which allows us to obtain each component independently: Using the most common form of the Hamiltonian: we can obtain the final expression for the D and E zero-field splitting parameters from the diagonal terms of the D tensor: This method was implemented in the NRLMOL code, using the PBE functional and its own basis set. As the calculations are restricted to molecular systems, they were performed on the isolated Fe 4 complexes using the experimental structures. [46][47][48] For the deposited systems, the molecular structure was taken from the structural optimization with the gold surface, but the zero-field splitting parameters were determined for the molecule alone, without the surface.
Transport calculations were performed using version 3.0 of the SIESTA code 35,36 using the PBE functional 37 with and double-zeta basis set for all the elements, with the exception of gold; for which a single-zeta basis was used combined with a 1-electron pseudopotential. 49 This 1-electron pseudopotential gives incorrect structures if it is employed for geometry optimization but reasonable transport properties in single-point calculations. 50  The most usual approach used to analyse the phase coherent transport in molecular transport junctions is based on the work of Landauer, Imry and Buttiker. 52 The expression for the conductance (G) for a given system with current I and voltage V is:

Exchange interactions in isolated Fe 4 complexes
In order to study the exchange interactions in the family of Fe 4 complexes, we selected six different complexes (see Figure 2

Exchange interactions of Fe 4 complexes deposited on gold surfaces
For the analysis of the variation in the exchange interactions of the Fe 4 complexes with the gold surfaces, two systems were selected: the complexes XUBWIB 60 and ICOCIN 59 (hereafter, Fe 4 C9 and Fe 4 ph, respectively). There is an important difference between the systems: in Fe 4 C9, the ligand functionalization with a thioester group provides efficient anchoring but with a long distance between the central Fe 4 moiety and the surface; while for the Fe4ph system, there are just van der Waals interactions. Due to the long length of the ligands of the Fe 4 C9 system, two possible interaction modes with the surface are possible ("standing-up" and "lying-down", see Fig. 3). The third molecule also employed in the molecular electronics experiments is an equivalent system but with a shorter ligand than Fe 4 C9: with only five methylene units (Fe 4 C5); thus, only standing-up is possible. 16 For the study of the interaction of the Fe 4 C9 complex with the Au(111) surface, we considered two different chemical structures of the anchoring group: (i) the thioester group maintains its structure during the interaction with the surface; and (ii) homolytic cleavage of the S-C bond occurs with the consequent generation of two radicals: Fe4-S· and an acetyl CH 3 -C=O.
Thus, two free acetyl radicals will form a diacetyl molecule; this is the mechanism that is usually proposed for thiol groups leading to the formation of hydrogen molecules from hydrogen radicals.
Optimized DFT structures (see Computational details section) gave an interaction energy of -6. . Hence, the most stable case corresponds to the lying-down structure (5.7 kcal/mol more stable than standing-up, with radical contacts) and the interaction through the radical groups (87 kcal/mol more stable lying-down with radical contact than with the thioester groups).  Table 2 for the structures with a radical S· anchoring group and the equivalent values for the unchanged thioester ligand. Comparison of the J values with those in Table 1 Figure 3 as the most stable. The phenyl rings are parallel to the gold surface which maximizes the contact between equatorial ligands and the surface, which is 10.1 kcal/mol more stable than the case with one phenyl ring oriented towards the surface. Intermediate tilted orientations of the molecule are much less stable. There is an important change in the structure of the complex due to the interaction with the molecule, but also the optimization of the geometry reduces the bridging angles in comparison with the X-ray data used for the values in Table 1. Such structural modifications are especially important in the case of a considerable decrease in the Fe-O-Fe angle or the antiferromagnetic coupling.

Zero-field splitting parameters of isolated Fe 4 complexes and those deposited on gold surfaces
The magnetic anisotropy of the isolated Fe 4 molecules was studied using DFT methods, including spin-orbit effects (see Table 3). Such an approach usually gives reasonable results for polynuclear complexes with low degeneracy, which allow a realistic description of the ground state with a single Slater determinant. 66 The D value that is often employed to quantify the magnetic anisotropy is correlated with the average of the γ angles between the planes involving the Fe 3 unit (the three external Fe centres) and Fe 2 O 2 unit (between the central and the three terminal Fe cations) planes (see Figure 4). This correlation with the γ angle as proposed by Cornia and Sessoli 67,68 is clearly reflected in the D values calculated using DFT. Thus, for the two isolated molecules, the larger angle in NIPJEC is also followed by a slight increase in the D parameter that is in reasonably good agreement with the experimental data in both cases. The E value also increases with the γ angle, but in the optimized DFT structures of the molecule deposited on the gold surface this parameter becomes larger. This fact should indicate larger tunnelling effects in the spin relaxation for the deposited molecules.

Transport properties of single-molecule Fe 4 devices
The transport properties 69,70 of the Fe 4 complexes were studied via two approaches: (i) one assuming a Coulomb blockade mechanism with a molecular redox process associated with the transport that was analysed by comparing DFT results for isolated neutral and charged molecules; and (ii) coherent tunnelling transport was studied using DFT+EGF (see Computational details section) and a model with the Fe 4 complex sandwiched between two periodic gold electrodes.
Magnetic anisotropy controlled the electric field as proposed by Zyazin et al. 18,33,48 is based on the change of the oxidation state of the ICOCIN Fe 4 complex from the S=5 neutral state to a reduced S=11/2 state with a change in the D parameter from -0.06 to -0.09 meV (-0.48 to -0.27 cm -1 ). We performed DFT calculations using the Gaussian code with the B3LYP functional (see Computational details) to compare the energies of the optimized neutral, reduced and oxidized NIPJEC Fe 4 complex. 62 At first glance, reduction of a Fe 4 complex is expected at the Fe III centres while the S values could be 9/2 or 11/2; while from a chemical point of view, it is expected that the additional electron will cause the reduction of the Fe III cations to give a high-spin S=2 Fe II centre.
These predictions were confirmed by the DFT calculations which also showed greater stability of the S=9/2 configuration (around 4843 cm -1 ) than the S=11/2 case. In the case of the oxidation process, due to the high oxidation state of the Fe III centre, the DFT results indicated that the process occurs on the oxygen atoms of the bridging ligands with a small energy preference for the S=11/2 case (the Fe 4 centres remain unchanged and an additional unpaired alpha electron is delocalized over the ligands). Furthermore, we calculated the D (E) zero-field splitting values for the three cases: the neutral S=5, the reduced S=9/2 and the oxidized S=11/2 NIPJEC Fe4 complex, using the optimized DFT structures; we obtained: -0.40 (0.034), -0.75 (0.007) and +0.47 (0.14) cm -1 , respectively. These results clearly corroborate the idea that the electric field can modified the redox state of the Fe 4 complex, as suggested by Zyazin et al., 18,33 not only by increasing the D value but also the sign of the D value could be changed.
In the second part of this section, the coherent transport properties of the Fe 4 complexes will be analysed with the help of DFT methods. The DFT study was performed using a GGA+U approach using the PBE functional (see Computational details section). The first system that we studied was the Fe 4 C9 system with a similar structure to that depicted in Figure 1. However, the capacity of such a system to transport carriers is rather limited due to the long aliphatic chain of the anchoring ligand. For instance, the current calculated for a voltage of 2 V is only 2 pA; thus, we focused on the Fe 4 ph system. Hao and coworkers reported a theoretical study of the transport properties using DFT+NEGFT but reducing the length of the Fe 4 C9 ligand chain. 71 The transport properties of Fe 4 ph system were already studied by adding some thiol groups in the phenyl rings to increase the strength of the molecule-electrode contacts. 72 In our case, we performed the study with the original non-functionalized structure (see Fig. 5) as employed in the experiments.  The results for Fe4ph are shown in Figure 6. The GGA+U provides more accurate results than usual GGA functionals (the geometry was optimised using a functional including dispersion contributions 39 ), because it corrects the self-interaction error that causes an underestimation of the energy difference between occupied and empty orbitals. Thus, in Figure 6, the occupied molecular orbitals with a large contribution from the Fe d orbitals are below -5 eV (see DOS) and consequently they are not very important for transport. Empty Fe d orbitals (at +1.5 eV) are closer to the Fermi level but are not relevant for transport, especially at low voltages, when it is mostly through occupied ligand channels. Thus, we should expect transport mediated mainly by holes and due to the proximity of this level to the Fermi level, relatively low bias will considerably increase the current between the electrodes. The current is slightly spin polarized, with alpha carriers predominant, as can be seen in the transmission spectra where such levels (in green) are closer to the Fermi level than the beta ones. Also, it is worth noting that the current at 2 V is around 1000 nA: several orders of magnitude larger than that found for the Fe 4 C9 system; and that despite the almost negligible influence of the Fe d "magnetic" orbitals on the transport properties, a slightly spin-polarized current is obtained due to the influence of such metal orbitals on the ligands.

Concluding Remarks
New spintronic devices based on magnetic molecules face the challenge of the problem of the electronic structure of the spinterface created between the anchored molecule and the surface. In The second magnetic property that we examined was the magnetic anisotropy by calculating D (and E) zero-field splitting parameters. Magnetic anisotropy is less sensitive to structural changes than exchange coupling constants are, thus, molecule-surface interactions cause only small variations in the calculated D values. The D value correlated, as previously reported, with the average γ angle between the planes involving the Fe 3 unit (the three external Fe centres) and Fe 2 O 2 unit (between the central and the three terminal Fe cations) planes. It is worth noting that despite the small E values calculated, as this magnitude is related to the strength of the spin relaxation tunnelling effects, these should be more important when the molecule is deposited on the surface. This result is in agreement with the common experimentally observed tendency of SMMs to reduce their blocking temperature (or lose the hysteresis loop) when deposited. Furthermore, some spin-phonon relaxation due to the vibrations of the surface will also decrease SMM behaviour.
The transport properties of the Fe 4 complexes were also analysed: firstly, by considering a Coulomb blockade regime where the charge carriers are "captured" by the molecule for a long time period that will imply a change in the oxidation state of the molecule and geometrical rearrangement. The reduction and oxidation processes of one S=5 Fe 4 complex were also studied using DFT methods and showed that in the reduction process, the extra electron is placed at the iron centres (reducing their oxidation state) and as the alpha Fe III d orbitals are completely filled, the electron must be in a beta orbital, resulting in a most stable S=9/2 state. In the oxidation process, as the Fe III centres are already in a high oxidation state, the bridging oxygen atoms lose one electron and for such a case, the S=11/2 state is the most stable. The magnetic anisotropy (D zero-field splitting) was calculated for the neutral, oxidized and reduced complexes showing an important variation of the D values not only in magnitude but even in the sign (changing from easy-axis to easy-plane for the oxidized structure). This result is quite important because it confirms that an electric field employed to reduce or oxidize the molecule can control the magnetic anisotropy of the system.
Finally, the coherent transport properties for the Fe 4 ph complex between gold electrodes were analysed in detail using periodic DFT methods combined with Green functions. The results indicated a very low participation of the transport channels involving a large contribution of the Fe III d orbitals, which are not close to the Fermi level. Despite such small contributions of the "magnetic orbitals" of the molecule, the transport (I/V characteristics) showed a spin polarization due to the influence of these orbitals on those of the ligands. The occupied alpha orbitals of the ligands are mostly responsible for the current; thus, positive carriers should be expected and they could be experimentally corroborated by measuring the thermoelectric effect.