Simple model of the ground state and spin-orbital excitations of free and adsorbed Fe(II) phthalocyanine molecules

We investigate the ground state and low-energy spin-orbital excitations of a single iron(II) phthalocyanine molecule in isolation and on an oxidized Cu(110) surface. Considering the subspace spanned by the three lowest spin-triplet states of 3 A 2 g and 3 E g symmetry, we diagonalize the Hamiltonian made of the anisotropic spin-orbit interaction and the ligand ﬁeld splitting (cid:2) , deﬁned as the energy difference between 3 E g and 3 A 2 g . We ﬁnd that the ground state switches from a 3 E g -like state with large orbital moment and out-of-plane easy axis for (cid:2) < − 60 meV to a 3 A 2 g -like singlet state with in-plane easy axis for (cid:2) > − 60 meV. The analysis of the ﬁrst excited states in the two regimes explains the zero-ﬁeld splitting data reported for β -FePc as well as for FePc molecules adsorbed on an oxidized Cu(110) surface [N. Tsukahara et al. , Phys. Rev. Lett. 102 , 167203 (2009)]. Importantly, the calculated magnetic susceptibility obtained with the ab initio value (cid:2) = 93 meV compares remarkably well with the experimental data of β -FePc


I. INTRODUCTION
Due to the importance of the Fe phthalocyanine (FePc) molecule for possible applications in nanodevices and spintronics, a quantum model of its magnetic ground state is of the utmost interest. There is, by now, converging theoretical evidence obtained with different methods that, neglecting spinorbit interaction, the ground state of an isolated FePc molecule is an orbital singlet of A 2g symmetry followed by two orbitally degenerate states of E g symmetry lying about 0.1 eV higher with spin S = 1 [1][2][3][4][5]. On the experimental side the analysis and interpretation of the Fe Kand L-edge x-ray magnetic circular dichroism (XMCD) in the case of a Au-supported FePc film carried out in Ref. [1] led to the conclusion that the above multiplets are indeed the constituents of the ground state of the isolated molecule, provided that the intermolecular interaction in the stack does not change the nature of the lowest multiplet levels in the single molecule.
It is therefore important to investigate the consequences of the model suggested above in comparison with experimental evidence available for the isolated molecule, namely, the temperature dependence of the paramagnetic susceptibility in the range of 1-300 K obtained by Dale et al. [6], Barraclough et al. [7] and Labarta et al. [8] in crystal powder samples of FePc in the β phase [9]. In this phase the molecular arrangement of the FePc molecules in a zigzag fashion with a large distance between metal ions results in weak mag-netic coupling between them, as proved by the fact that the sample remains paramagnetic down to the lowest measured temperature (1.5 K). The measured magnetic properties can therefore be ascribed to the single molecules with reasonable accuracy. They provide a clue to the nature of the ground state and give access to the first excited state of the molecule, the so-called zero-field splitting (ZFS). Low-energy excitations of FePc have also been reported from inelastic electron tunneling spectroscopy with scanning tunnel microscopy (STM) where the FePc molecules were adsorbed on an oxidized Cu(110) surface [10]. As the substrate perturbs the electronic structure of the molecule, it is interesting to see whether the model developed in Ref. [1] for a FePc film can be extended to this case. This is particularly important in view of technological applications, where FePc molecules will be subject to external perturbation, especially at interfaces.

II. THE MODEL FOR THE ISOLATED MOLECULE
As discussed in the Introduction, our starting assumption is that the ground and first excited states can be described as mixtures of the three lowest-lying multiplets, which are the orbital singlet A 2g and the orbital doublet E g , all being spin triplets, i.e., nine states altogether. All other multiplets are well separated in energy, typically by 0.5 eV or more according to Table III of Ref. [1]. The spin-orbit interaction (SOI), which is one order of magnitude smaller (∼50 meV), cannot significantly mix in these higher multiplets.
Here we determine the ground state, ZFS, and magnetic susceptibility by diagonalizing the Hamiltonian in the subspace of the three lowest multiplets ( 3 A 2g , 3 E g ) in the presence of SOI. This model was already considered in Sec. VI of Ref. [1], but only in the limit of saturated spin moments, which is appropriate for FePc films in the α phase with a large intermolecular exchange field. Here we focus on the isolated molecule and therefore solve the model without the constraint of spin saturation.
In general the Hamiltonian reads Here is the multiplet Hamiltonian, and E(A 2g ) and E(E g ) are the multiplet energies of the spin-triplet states. We define = E(E g ) − E(A 2g ) and take E(E g ) to be the energy zero, such that E(A 2g ) = − . The spin-orbit Hamiltonian H so has the usual form, H so = i ζ i l i · s i , where the sum runs over all the electrons i and s = (s x , s y , s z ) are the spin-1/2 matrices. For the SOI, the contribution of the light ligand atoms is negligible, and the calculation of the SOI parameter ζ i can be restricted to the Fe atom [11]. For a given molecular orbital i, ζ i is proportional to the fraction of charge λ i that lies inside of the muffin-tin sphere of Fe, and we find ζ (e g ) = 0.04 eV [1]. Finally, H ext describes the perturbation due to the substrate, which we model by introducing matrix elements between the A 2g and E g basis states (see Sec. III). For the isolated molecule, we have H ext = 0.

A. Analytic solution of the model at zero field
Indicating by d s the molecular orbital (MO) of symmetry s in the D 4h point group [2], the nine basis states are given by the following Slater determinants (see Ref. [1], Eq. (67)): The basis states are eigenstates of S z with eigenvalues M S = 1, 0, −1. They are ordered here as A 2g , E 1 g , E 2 g for fixed M S and then with decreasing M S . It is easy to see that all basis states have orbital momentum L z values between 1 and −1 and that it is possible to construct simultaneous eigenstates of S z and L z , with M L = 1, 0, −1. These transformed basis states are labeled |M L M S in Table I and classified according to When calculating the matrix elements of SOI between these states, one has to remember that the MOs are mixtures between Fe 3d and ligand orbitals. A crucial quantity for our analysis is the relative 3d weight in each MO, which we denote as λ s , where s labels the orbital symmetry. In multiple-scattering theory the wave functions are expanded in partial waves around each muffin-tin sphere, and the Fe d contribution can be written as s α s R d (r )Y 2s (r), from which the Fe d weights are obtained as where r 0 = 1.10 Å is the muffin-tin radius of the Fe atom. For an occupied MO, λ is the charge lying inside the Fe atomic sphere. Note that for free-atom Fe 3d orbitals, λ is as large as 0.98, so the λ values of the MOs accurately describe the reduction of 3d weight through hybridization. In particular, for a MO of E g symmetry λ e is found to be 0.59, whereas for a MO of A 2g symmetry λ a = 0.79. Therefore a matrix element of the SOI connecting two states of E g symmetry is proportional to λ e , whereas one connecting two states of respective symmetries E g and A 2g is proportional to √ λ e λ a times a factor of √ 3 coming from the application of l x,y [12]. Since the factor λ e is already incorporated in the definition of ζ , the ratio between the two matrix elements turns out to be r = √ 3λ a /λ e ≈ 2. The Hamiltonian matrix falls into blocks of dimensions 1, 2, 3, 2, 1 for M J = 2, 1, 0, −1, −2, which can be diagonalized analytically. Due to time-reversal symmetry, states with M J = 0 are doubly degenerate. Table II gives the matrix elements i|H so |j between the states in Table I Table III. Figure 1 (top panel) plots the three lowest eigenvalues, corresponding respectively to states |0 − , |1 − , and |2 as a function of the ligand field strength . We see that the three eigenvalues cross at x = 2 /ζ = −3 ( = −60 meV). At this point the ground state switches from an 3 E g -like, M J = ±2 doublet with large out-of-plane magnetic anisotropy to the 3 A g -like orbital and spin singlet |0 − with in-plane magnetic anisotropy (see the discussion of the magnetic susceptibility below). The bottom panel replots the three lowest eigenvalues, putting the ground state at E = 0 for all 's, visualizing the first and second zero-field splittings (ZFS1 and ZFS2; see below).

B. Analysis of magnetic susceptibility data
To establish the magnetic anisotropy for the singlet ground state |0 − , we note that 0 − |L α |0 − = 0 − |S α |0 − = 0 for α = x, y, z. Therefore one has to go to second-order perturbation when calculating the variation of the energy with magnetic field, i.e., the magnetic susceptibility tensor where |n is any excited state of the molecule and we have used the fact that in D 4h symmetry, χ αβ is diagonal and χ xx = χ yy = χ zz .
The relevant matrix elements different from zero are Therefore with the factor of 2 in the second equation coming from the double degeneracy of the excited state |1 ± . Figure 2 plots on a log scale χ z,z and χ x,x as a function of . For = 93 meV, the ratio χ x,x /χ z,z = 3.1 × 10 3 , indicating a strong in-plane anisotropy.
A quantity that can be compared with the experimental data is the ZFS, i.e., the energy difference between the ground state and the next excited state in zero magnetic field. The bottom panel of Fig. 1 plots the splittings between the ground state and the next two excited states (top panel of Fig. 1) as a function of , which we call ZFS1 and ZFS2 (ZFS1 . In Ref. [2], for the case of the isolated molecule, we obtained = 93 meV from ab initio molecular orbital calculations, so that from Fig. 1 we find ZFS1 = 8.8 meV and ZFS2 = 96.6 meV. However, rather than rely on the comparison of a single value, we now have all the ingredients to calculate the average magnetic susceptibility in the case of an isolated molecule and compare it with the experimental points taken by Dale et al. [6] in the temperature range between 1 and 300 K. We also compare it with the data by Barraclough et al. [7] and Labarta et al. [8]. To write down the magnetic susceptibility we need to expand each stationary energy state E α,j (H ) in Table III in terms of the magnetic field strength H as where j indicates the degeneracy. Indicating for brevity by n the pair α, j , the magnetization density associated with this state is In the limit H → 0 we can write n H/kT , (9) so that Eq. (8), using Eq. (7) and keeping only terms up to first order in H , becomes Ignoring saturation effects, one usually considers the part of the susceptibility χ = M (H, T )/H which is independent of the field strength. Denoting by N Avogadro's number and resuming again the pair of indices α, j , the molar susceptibility can be written as [13] provided or, in other words, that the system is in a paramagnetic state.
To calculate explicitly Eq. (11) we first observe that, since the excitation energy corresponding to the highest measured temperature (300 K) is around 25 meV, it is sufficient to consider the first excited state |1 − (ZFS1 = 8.8 meV) in the summation over α (apart from the ground state |0 − ) and to neglect the next excited state |2 (ZFS2 = 96.6 meV).
Turning now to the calculation of where β indicates the Bohr magneton. For E (2) α,j we need to consider two terms. The first one arises when α runs over the ground state |0 − , already calculated in Eq. (5) (without the degeneracy factor 2 and taking only the δ = − term of the sum) whereas the second one comes from identifying α with |1 − , Introducing the quantity d = D/kT = (ω − 1 − ω − 0 )/kT , we therefore obtain for the average susceptibility This expression should be compared with the one used by Dale et al. [6], to fit the parameter D to the experimental points. This is the usual expression (neglecting spin-independent terms) derived on the basis of a spin Hamiltonian for an orbitally nondegenerate state with effective spinŜ = 1 and magnetic field either parallel H z or perpendicular H x,y to the molecular axis [12], with respective eigenvalues Indeed, expansion of these eigenvalues in powers of the magnetic field provides E (1) α,j and E (2) α,j in Eq. (6), which inserted in Eq. (11) leads to Eq. (15).
According to the theory of the spin Hamiltonian, Dale et al. [6] used the relation 2D = λ(g ⊥ − g ), where λ = 18.6 meV was the spin-orbit coupling parameter λ L · S for the ground term of the free-ion value, reduced by 25%. Since S = 1, this quantity is related to our ζ by the relation λ = ζ /2 and is therefore very close to our value of ζ /2 = 20 meV. From the fit, they derived D = 8.67 meV, g = 1.93, and g ⊥ = 2.86. In their words, "any value taken from the range 17.3-21.0 meV will yield a reasonably close fit and acceptable values of g and g ⊥ . Moreover, we have found that the value obtained for D is not a sensitive function of the value chosen for λ." In our case, using our ab initio values = 93 meV and ζ = 40 meV, we obtain D = 8.83 meV, g = 1.88, and g ⊥ = 2.89. The following expressions for g and g ⊥ were derived from the comparison of the two expressions of χ m (T ) (14) and (15) and Eq. (5): With this identification, the two expressions for the susceptibility are identical, and in the high-temperature limit one finds a Curie-Weiss law [12], Inserting the ab initio values given above, one finds θ = 8.1 K, which should be compared with θ = 9 K reported by Lever [14], θ = 7 K given by Labarta et al. [8], and θ = 3 K reported by Dale et al. [6]. The discrepancy between these experimental values might be due to the different numbers of experimental points used by the various authors, 13, 23, and 4 in the range 90-300 K, respectively, to fit the Curie-Weiss law.
Note that the similarity between expressions (14) and (15) is not fortuitous. In both cases the low-energy excitation structure of the underlying models in the presence of an external magnetic field is the same and is described by Eq. (16), in one case with effective operatorŜ = 1 and in the other with effective operatorĴ =L +Ŝ = 1, with respective angular momentum projections M S and M J . The advantage of our approach lies in the fact that it gives an explicit expression for the spin-orbit coupled eigenstates, allowing the ab initio calculation of the various parameters of the effective Hamiltonian D, g , and g ⊥ . Figure 3 compares the theoretical curve, Eq. (14), with the experimental points by Dale et al. [6], Barraclough et al. [7], and Labarta et al. [8] for = 93 meV and two values of the anisotropy parameter, r = 2 and r = √ 3. The agreement between the theoretical curve for r = 2 and the experimental data is remarkable, considering that there is no adjustable parameter in the theory. All relevant quantities ( , ζ , r, λ e , λ a ) have been calculated ab initio. For comparison the curve with r = √ 3 is also shown, corresponding to neglecting the anisotropy of hybridization (λ a /λ e = 1). The agreement with the data is considerably worse. This shows that the common practice in ligand field multiplet calculations, which consists of treating hybridization effects through an isotropic reduction factor, is very questionable. In particular in the case of the FePc molecule the anisotropy of the hybridization parameters was found to be essential for establishing the correct ordering of the A 2g and E g multiplets in the isolated molecule and their energy gap . Finally, Barraclough et al. [7], besides measuring the average susceptibility, also provide measurements of the magnetic anisotropy of FePc from room temperature to about 90 K. Appendix A discusses these measurements in relation to the theoretical predictions.

III. ZERO-FIELD SPLITTING OF FEPC ADSORBED ON AN OXIDIZED CU SUBSTRATE
In the previous section on the isolated molecule, the external perturbation H ext in Eq. (1) was set to zero. Here we extend our model by introducing a perturbation due to a substrate and analyze the ground and first excited states of an FePc 195108-5 molecule adsorbed on an oxidized Cu substrate [10]. Single FePc molecules adsorb essentially flat on the Cu surface, but some tilting and/or distortion cannot be excluded [10]. On physical grounds, two kinds of effects can be expected from the interaction with a substrate. First, the ligand field parameter = E(E g ) − E(A 2g ) is modified and may even change sign. Second, the symmetry is reduced, which may lead to mixing between the D 4h molecular orbitals. Concerning the first effect, the a g (= d z 2 ) orbital, which is oriented normal to the molecular plane, hybridizes more strongly with the substrate atoms than the e g (d xz , d yz ) orbitals. Being an antibonding orbital [15], the a g level is pushed up with respect to the isolated molecule.
As a consequence, the |A 2g multiplet rises relative to |E g because |A 2g has one more electron in a 2g [see (2)]. This argument is confirmed by the constrained density functional theory calculations of Nakamura et al. [3], who found a ground state of type |A 2g for the isolated molecule and one of type |E g for a film of molecules. In other words, the ligand field parameter has changed sign.
Tsukahara et al. [10] measured the ZFS of the FePc molecules adsorbed on a Cu(110)(2 × 1)-O surface by inelastic electron tunneling spectroscopy with STM as a function of an applied external magnetic field. They observed a reduction of the ZFS value with respect to the isolated molecule and a change in magnetic anisotropy from easy plane to easy axis. Two FePc species, i.e., two different adsorption states, were found to coexist on the substrate: an α species with ZFS1 = 1.9 meV and ZFS2 = 4.7 meV and a β species with ZFS1 = 4.1 meV and ZFS2 = 9.0 meV [16]. These ZFSs can roughly be reproduced within the easy-axis phase of our free-molecule model by taking = −72.5 meV (ZFS1 = 2.0, ZFS2 = 3.5) and = −86.7 meV (ZFS1 = 4.2, ZFS2 = 7.4), respectively (see Fig. 1, bottom panel). However, as the easy-axis ground state is doubly degenerate due to timereversal symmetry (corresponding to M J = ±2 in Table I), one would expect each ZFS value to split into two branches upon application of a magnetic field, whereas experimentally, only one is observed in both species [10]. This finding is explained by the second effect, i.e., symmetry reduction upon adsorption. In the free-molecule model (1) with H ext = 0, i.e., the multiplet Hamiltonian in D 4h symmetry restricted to the { 3 A 2g , 3 E g } subspace, J z = L z + S z is conserved. The substrate interaction lowers the symmetry, which leads to mixing between multiplet states of different M J and thus to lifting of the ground-state degeneracy M J = ±2, as we are going to show.
The substrate perturbation H ext is assumed to be an effective one-particle operator. By extending the model used in Sec. VI of Ref. [1], it can be defined in the basis introduced in (2) via the matrix elements with the obvious meaning of the spin projection symbols u, 0, d. The rationale underlying this choice is that, because of the symmetry reduction upon adsorption, the MO |d xz and |d yz are no longer equivalent, and they mix with |d z 2 via the matrix elements d xz,yz |H ext |d z 2 . This hybridization is assumed to be spin dependent, which can be related to the antiferromagnetic coupling found between the Fe 3d and Cu 3d bands [17]. By noting that with similar relations for t σ 2 upon interchanging d σ yz with d σ xz , the definition of the perturbation in (20)  Guided by the findings in Sec. II for the free molecule, we keep < −60 meV in the vicinity of the transition point. We start with the β species. We have set all t σ 2 to zero for simplicity. By trial and error we have obtained a good fit of the ZFS and magnetic properties with the parameter values t u 1 = 35, t 0 1 = 25, t d 1 = 25, and = −65 meV. The three lowest eigenvalues are −35.15, −31.07, and −26.05 meV, corresponding to a ZFS1 of 4.08 meV and a ZFS2 of 9.1 meV, in excellent agreement with the experimental values of the β species (4.1 and 9.0 meV) [10]. The eigenvectors are listed in Table IV in Appendix B. Regarding the spin character, for the ground state (GS) we find GS|S z |GS = 0.49, and for the first excited state (EX) we find EX1|S z |EX1 = −0.40, in approximate agreement with the χ ± states of Fig. 2 of Ref. [10], but for the second excited state we have EX2|S z |EX2 = −0.30, at variance with the spin character of the upper state χ 0 . However, given the approximate character of the effective spin model used in Ref. [10], the correspondence is acceptable [18]. Note that the spin averages over a given state are those projected onto the Fe site. Those without projection are easily obtained from the mixing coefficients of the state in terms of the basis states. For the orbital moment of the GS we find GS|L z |GS = 0.46μ B , GS|L x |GS = 0, and GS|L y |GS = 0.45μ B , so that μ z /μ ⊥ = 2, indicating the easy-axis character of the GS, in agreement with the negative D value found for the molecule on the oxidized Cu surface.
The spin character of the states is similar to that of the β species. GS|S z |GS = 0.45 and EX1|S z |EX1 = −0.36, but EX2|S z |EX2 = −0.16. The orbital character of the GS is again easy axis, with μ z /μ ⊥ = 2.32 since GS|L z |GS = 0.52 μ B , GS|L x |GS = 0, and GS|L y |GS = 0.45μ B . The different sets of parameters needed to fit the excitation spectra of the two species reflect their different couplings with the same substrate, as observed experimentally. It is an indication of the internal consistency of the model that is the same for both species and the hybridization parameters t σ 1 of the α species are roughly half those of the β species.
Moreover, the need of a spin-correlated hopping to fit the data signals a kind of spin coupling between the molecule and substrate. The model is flexible enough to accommodate different situations but has no predictive character, although it can shed light on some aspects of the physics of the coupling. To make some progress, the relation between the hybridization parameters and the molecular structure still needs to be investigated, as well as their influence on the value of .
Finally, we have calculated the behavior of ZFS1 and ZFS2 when an external magnetic field B z (T ) is applied to the molecule, from 0 to 10 T (Fig. 5). All parameters are fixed by the above fits. A comparison with Fig. 2 of Ref. [10] is quantitatively good for the α2 and β1 branches but not for the α1 and β2 branches, which show, respectively, a marked upward and downward curvature not reproduced by the calculations. One reason for this discrepancy might be that under the external field the coupling of the molecule with the substrate, being spin dependent, slightly changes with the intensity of the applied field. As a last comment we observe that our model makes sense provided the perturbation brought about by the substrate is weak, in the sense that the electronic structure of the molecule is not substantially altered. This is not the case of the clean Cu(110) surface. As the work function of Cu ( 5eV) is much smaller than the binding energy (with respect to vacuum) of the free-FePc lowest unoccupied molecular orbital (LUMO) level (13 eV in our MS calculation [1]) the LUMO level will be pinned to the Cu Fermi level in the adsorbed molecule. This leads to charge transfer from the surface to the molecule which will partially fill the e g and a 1g holes and thus reduce the Fe spin and orbital moment, possibly quenching it completely. Indeed, some states (especially d z 2 states) become so strongly hybridized with the Cu band that they undergo a large broadening and spin/orbital polarization is lost.

IV. CONCLUSIONS
Based on the analysis of the XMCD spectroscopic data at the Fe K and L edges and recent theoretical investigations [1][2][3][4][5], we have presented a model of the magnetic ground state and low-lying excited states of the FePc molecule in which an orbital singlet of A 2g symmetry and two orbitally degenerate multiplets of E g symmetry separated by 93 meV are mixed by a spin-orbit interaction of the order of 50 meV. Relying on the level scheme summarized in Table III and the corresponding matrix elements of the magnetic operators L and S, we have calculated the paramagnetic susceptibility in the Van Vleck model. The remarkable agreement with the experimental data of Dale et al. [6] (complemented by the data by Labarta et al. [8]) in the range of 1-300 K obtained ab initio without adjustable parameters makes us confident that we have found the low-energy quantum structure of the free FePc molecule. We thus substantiate all the inferences contained in Ref. [6] regarding the ground and first excited states of the molecule. As stated in their abstract, "in the range 1.25-20 K the susceptibility is virtually independent of temperature and this result has been used to show that the ground state of the central iron atom is an orbital singlet, with a real spin triplet state. This is split by second-order spin-orbit coupling into a singlet ground state and a doublet state at 70 cm −1 , to give g = 1.9 and g ⊥ = 2.9." All this is reflected in our finding that the ground state |0 − is an orbital singlet made up of spin-triplet states combined in such a way that 0 − |S z |0 − = 0, followed by a doublet |1 − at 71.5 cm −1 , for which 1 − |S z |1 − = cos 2 θ − and 1 − |L z |1 − = sin 2 θ − .
Moreover, the disagreement of the theoretical curve for r = √ 3 with the paramagnetic susceptibility data provides us with the experimental proof that the usual practice of reducing the free-ion Coulomb parameters by a uniform quantity to account for hybridization in a molecular environment does not give accurate results. 195108-7 Using the eigenvalues listed in Table III, the eight excited states of the molecule are calculated as 8.8 meV (doublet), 97 meV (doublet), 132 meV (doublet), 137 meV (singlet), and 162 meV (singlet). A comparison with the same spectrum calculated in Ref. [19] as 2.4 meV (singlet), 3.6 meV (doublet), 11 meV (singlet), 47 meV (singlet), 130 meV (doublet), and 160 meV (singlet) provides further evidence, besides that given in Ref. [1], that their suggested ground state is not correct.
Finally, we have successfully reproduced the experimental data by Tsukahara et al. [10] by application of an extension of the model used in Ref. [1]. The results confirm that the type of magnetic anisotropy of the molecule on a substrate is mainly determined by the sign of the energy gap between the orbital doublet and the singlet, while the hybridization between them (the parameters t σ i ) is responsible for lifting any degeneracy left by SO coupling and for influencing the averages of the spin and orbital operators.
In Fig. 3 the data for the average susceptibility obtained by Barraclough et al. [7] are plotted against those measured by Dale et al. [6]. The slight difference between the two sets probably comes from their different ways of treating the data. Dale et al. found that the sum of the constant contributions to the molar susceptibility is very close to zero, indicating that contributions from diamagnetism and temperatureindependent paramagnetism cancel each other out to a large extent. Barraclough et al. instead corrected for diamagnetism using the value given by Lever [14].
At first glance, one might suppose that the fitted values for D, g ⊥ , g derived from the two sets are roughly the same. However, Barraclough et al. measure not only the average susceptibility χ m but also separately the susceptibilities perpendicular and parallel to the symmetry axis of the molecule in a crystalline sample. In practice, the anisotropy components (χ 1 − χ 2 ) and (χ 1 − χ 3 ) of the crystal tensor χ 1 , χ 2 , χ 3 with respect to the principal axes 1, 2, and 3 were measured, and from them the actual β-FePc molecule anisotropy components χ ⊥ m and χ m were derived. Clearly, χ m = (2χ ⊥ m + χ m )/3.

By making the identification
from the fit, they deduced the values D = 7.94 meV and g ⊥ = g = 2.74, without having to resort to the SO coupling constant λ. These values seem to be rather close to those derived by Dale et al. [6] (D = 8.67 meV and g ⊥ = 2.86, g = 1.93); however, they pose a problem. In fact, one of the tenets of the theory underlying the derivation of the spin Hamiltonian for an orbitally nondegenerate ground state, used to interpret the susceptibility data, is the relations where ⊥ and are the diagonal components in tetragonal symmetry of the second-order tensor ab associated with the product of the matrix elements of the components of the orbital moment in the basis functions of the atomic term [12]. Therefore g ⊥ = g is incompatible with a finite ZFS D.
We now consider the theoretical quantities corresponding to those in Eq. (A1), which are compared in Fig. 6 with the experimental points.
One can see that while the theoretical χ ⊥ m is acceptably close to the experimental points, χ m is roughly a factor of 2 off the measured quantities.
In trying to understand the possible origin of this discrepancy we have reanalyzed the data by Barraclough et al. [7]. In particular we have noted that they use the relation, deduced from their Table III, χ ⊥ m − χ m = (χ 1 − χ 2 ) c + (χ 1 − χ 3 ) c , which is at variance with that given by Gregson and Mitra [20], χ ⊥ m − χ m = 2 (χ 1 − χ 2 ) c − (χ 1 − χ 3 ) c . On the other hand, this latter relation is in keeping with that derived by Bleaney et al. [21] for an FePc molecule with D 4h symmetry in a monoclinic crystal (β phase).
Due to these inconsistencies we have decided not to take into account Barraclough et al.'s anisotropic magnetic data in our analysis.

Hamiltonian matrix elements
Here the matrix elements of the Hamiltonian (1), H 0 + H so + H ext , in the basis (2) are given in units of ζ /2. The matrix elements of H ext are listed in Eq. (20). The term H 0 , i.e., the ligand field splitting of the free molecule, is diagonal in this basis. We set E(E g ) = 0, so the only nonzero elements are The matrix elements of the spin-orbit interaction H so are given by (see Eq. (68) of Ref. [1] calculated at θ = 0) 1|H so |5 = 4|H so |8 = − 2|H so |4 = − 5|H so |7 = i r √ 2 exp (iφ), The φ dependence in Eq. (B2) stems from rotating by the same angle the reference frame around the z axis in view of averaging the in-plane orbital moment over the random positions of the molecule. In this case μ ⊥ = ( L x + L y )/2 would depend only on the combinations (t σ 1 ) 2 + (t σ 2 ) 2 for each spin direction. In the following, however, we shall set φ = 0 since in the experimental setting of Ref. [10] there is only a single molecule with a definite orientation with respect to the CuO chains of the substrate.

Eigenvector coefficients of low-lying molecular eigenstates
For the β species of the adsorbed FePc molecule, we obtained a good fit with the parameter values t u 1 = 35, t 0 1 = 25, t d 1 = 25, = −65 meV. The eigenvalues and eigenvectors of three lowest-lying states are given in Table IV. For the α species, we took t γ 2 = 0 and the best fit values t u 1 = 15, t 0 1 = 10, t d 1 = 10, and = −65 meV. The eigenvalues and eigenvectors are listed in Table V.