Determination of the ground state of an Au-supported FePc film based on the interpretation of Fe K - and L -edge x-ray magnetic circular dichroism measurements

J.B. and F.B. acknowledge the financial support of the Spanish MINEICO MAT2017-83468-R project, Aragonese E-12 17R RASMIA (co-funded by Fondo Social Europeo), and of the European Union FEDER (ES). We acknowledge support by the European Union COST MP1306 EUSPEC project.


I. INTRODUCTION
Since their discovery and the subsequent systematic studies of their molecular structure, phthalocyanines have been the subject of special research interest because of their multiple applications such as dyes, catalysts, and coatings. At present, they are one of the most studied organic materials for possible applications in nanodevices and spintronics.
Among the phthalocyanines, Fe-phthalocyanine (FePc) molecule is a promising candidate as new ideas and methods can be applied, due to its simple and robust planar structure. Moreover, the spin moment can be chemically manipulated by oxidation [1] and Li doping [2], while its strong magnetic anisotropy [3] can be switched from planar to perpendicular by the application of a small external electric field through the magnetoelectric effect [4]. It is therefore of extreme importance to have a model of the magnetic ground state of the molecule, an objective which has been elusive for decades (see Refs. [5][6][7][8]).
In an isolated molecule, the crystal field (CF) on the Fe site, with approximate D 4h symmetry, splits the Fe 3d states into three orbital singlets (d xy ,d z 2 ,d x 2 −y 2 ) and a doublet (d xz ,d yz ), which yield the basis for the formation of molecular orbitals (MO) of the corresponding symmetry. Ignoring the antibonding orbital of x 2 − y 2 symmetry lying at too high energy [5], out of the four remaining states one can construct six spin triplets, in Table I. where 3 E 1,2 g and 3 E 1,2 g are orbitally degenerate. Which one is the ground state in the molecule is * pkruger@chiba-u.jp still being questioned, and this uncertainty affects any model one can construct out of these building blocks for a molecule in an epitaxial film.
In density functional theory (DFT) calculations, the theoretically predicted ground state of the FePc molecule/film depends on computational details such as the exchange correlation potential and the way self-consistency is achieved (e.g., constrained DFT) [5]. Virtually all the states listed in Table I have been suggested as ground states of the isolated molecule.
Recently, measurements of the x-ray magnetic dichroism (XMCD) at Fe-L 23 absorption edge in a Au-supported FePc film showed a large, unquenched planar orbital moment of 0.532 ± 0.004 μ B , a perpendicular orbital moment of 0.297 ± 0.007 μ B , and a value of 2 S z equal to 0.640 ± 0.003 [3,9]. The FePc molecules in the film (with a thickness of 133 nm) lie parallel to the Au/sapphire substrate and are stacked in chains (see Fig. 1).
While accounting for the large orbital moment anisotropy of around a factor of two, the interpretation of these data could not answer the crucial question of the ground electronic configuration. Quite naturally, these findings prompted a further investigation of the Fe K edge XMCD, since this signal is known to be sensitive only to the orbital moment [10]. Based on these measurements and on ab initio molecular orbital multiplet calculations of the isolated FePc molecule, we have proposed a model for the magnetic ground state of the FePc film that explains the Fe K edge XMCD data and reproduces the observed values of the orbital moments in the perpendicular and planar directions [11]. However, the expression of the ground state was introduced as an Ansatz and no attempt was given to derive it in a quantum mechanical way. The main objective of this paper is to present the quantum mechanical derivation of the ground state of FePc molecules in the Au supported thin film as a result of the diagonalization of a model Hamiltonian, accounting in this way for the experimentally observed orbital anisotropy values and XMCD measurements both at the Fe K and L 23 edges [3,11]. At the same time, we give arguments to rule out possible alternatives based on the states listed in Table I, in particular, the ground state recently proposed in Ref. [8].

II. EXPRESSION FOR THE FePc GROUND STATE FOR THE ISOLATED MOLECULE AND THE THIN FILM
The molecular orbitals of FePc have been obtained using the multiple scattering code by Smith and Johnson [12] in a spin-restricted calculation. We focus on the antibonding states with mainly Fe 3d character, which we label by the d-orbital symmetry as 3z 2 − r 2 , x 2 − y 2 , xy, xz, and yz, corresponding to the irreps a 1g , b 1g , b 2g , and e g of the D 4h point group, respectively. The computed energy levels are E = −13.34eV + E, where E is given in Table II. The x 2 − y 2 orbital, which forms direct bonds with the ligand N atoms, is 2.5-3 eV higher than the other four orbitals. Among them, xy is the lowest. These level splittings are in good agreement with those calculated with a projector-augmented wave code (VASP), given in parenthesis in Table II.
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. The calculated λ values are listed in Table II and agree within 5% with those obtained using the VASP code. 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. As the splittings between the four lowest orbitals in Table II are very small (<0.35 eV), it is essential to take account of Coulomb interactions beyond DFT for establishing the nature of the ground state and the first excited states. From the five MOs, times a spin function, denoted by ± subscript, we build Slater determinants (SD) as basis states for the multiplet calculation and the spectroscopic analysis. In the SD, the single-particle states are ordered by energy with spin up before spin down. We consider the space of all C 10 6 = 210 Slater determinants (SD) in which six out of the ten spin orbitals are occupied, corresponding to a nominal Fe 2+ configuration. In this sub-space, the many-body Hamiltonian containing the MO levels and the Coulomb interaction is diagonalized. The Coulomb matrix elements are calculated ab initio from the MO wave functions, omitting the monopole term of the multipole expansion, which is already included in the MO energies through the Hartree potential in the multiple scattering calculation. An advantage of such a procedure over the usual one of taking the Coulomb matrix elements from an atomic calculation (iron in the case of FePc) with a uniform reduction factor of about 30% to simulate the ligand effect is that in our approach the anisotropy of the various molecular orbitals is automatically taken into account in the calculation of the Coulomb matrix elements. Since such anisotropy factors vary from 0.84 to 0.6 (see Table II) their effect can be substantial on the ordering of the various molecular multiplets. Table III shows the lowest lying eigenstates and excitation energies (eV) referred to the ground state so obtained. Spaces between spin-orbitals in a SD have no physical meaning but FIG. 1. Left panel: scheme of the experimental molecular structure of a thick film (TF) of α-FePc on a Au substrate [3]. Fe atoms are indicated in red. Right panel: top view of the molecular stack with Fe in gold, N in silver, and C in brown (H is omitted). In agreement with Ref. [5], we find for an isolated molecule, by the method of molecular multiplets, that the ground state is an orbital singlet 3 A 2g , followed by an orbitally degenerate state of 3 E g symmetry 0.093 eV higher. However, the first state of 3 B 2g symmetry (nonorbitally degenerate), which involves the excitation of a spin-down orbital d − xy into an orbital of E g symmetry, lies 0.82 eV higher, at variance with Ref. [5], which reports a ten times lower energy. Our order is in keeping with the fact that the molecular orbital energy of the d xy orbital is lower than d xz,yz by 0.32 eV. The points stand for three singlet and four quintet states (S = 2), with excitation energies between 0.82 and 1.45 eV. The second 3 E g is mentioned here because it was suggested as a possible ground state in Ref. [6] in a parametric ligand field approach and in Ref. [8] in a ligand field multiplet model. On the other hand, the possibility of 3 E g and 3 B 2g states lying very close in energy and being mixed by spin-orbit coupling has been suggested in Refs. [2,8]. These assignments will be discussed in Sec. V in comparison with our findings.

155139-2
We now assume that, due to the perturbation of the nearest layers in the stack and the spin-orbit interaction, the two lowest molecular states A 2g and E 1 g , E 2 g mix together in the new ground state. Even though in the stack the symmetry D 4h is broken and the molecular orbitals broaden into unidimensional bands along the axis of the stack, in keeping with Ref. [5] we continue to label the band states according to the local molecular symmetry D 4h , in consideration of the fact that the magnetic behavior originates from the Fe atom, so that the main band amplitudes bear the old symmetry. In our approach, the band effects in the film are neglected, in the sense that the dispersion of the molecular orbitals that constitute the Slater determinants representing the various multiplets in the isolated molecule is not taken into account. Consequently, we represent band states of a given symmetry dispersing in a certain energy interval by a single state located at the centroid of the band and coinciding with the molecular orbital of the same symmetry in the isolated molecule. This is a very reasonable approximation, because the bands most affected by the stacking (of symmetry d xz,yz and d z 2 ) disperse no more than 0.4 eV, as shown by Fig. 4 of Ref. [5], whereas the bands of symmetry d xy are almost dispersionless. A bandwidth of 0.4 eV in a unidimensional band means a hopping integral of 0.1 eV between two adjacent molecules, which can be considered a weak perturbation on the molecule under consideration. Therefore we expect that in general the ground state of the molecule in the stack can be represented as where the coefficients c i (i = 0,1,2) are in general complex and such that i |c i | 2 = 1. In Ref. [11], we have proposed a model for the magnetic ground state of the FePc film that explains the quadrupole XMCD data at the Fe K edge and reproduces the observed values of the orbital moments in the perpendicular and planar directions. The model was based on the following expression for the ground state: which is a particular form of Eq. (1). As anticipated in Introduction, this expression was introduced as an Ansatz and no attempt was given to derive it as a result of the diagonalization of a model Hamiltonian. In this paper we shall fill this gap.

III. DIAGONALIZATION OF THE MODEL HAMILTONIAN. ANGULAR DEPENDENCE OF THE ORBITAL MOMENT
As a preliminary step for the diagonalization of the model Hamiltonian, we compute the matrix elements of the angular momentum operator L α , α = x,y,z, between the three nearly degenerate multiplets, Here, the d xy orbital, which is always double occupied, has been omitted. Using |l y |d xz , and the angular momentum matrix elements in a real spherical harmonics basis [13] (see also Table I of Sec. I in Ref. [14]), we find the results given in Table IV, where indices λ a and λ e correspond to a 1g and e g orbitals, respectively. The parameter r ≡ √ 3λ a /λ e = 2.0 is a sort of anisotropy factor and shall play an important role in the following. We have r = | A 2g |L x,y |E 1,2 g / E 1 g |L z |E 2 g |. The factor of √ 3 is intrinsic and originates from the application of l x,y to the orbital d z 2 (see Sec. I of Ref. [14]).
From Table IV. it can be seen that the groundstate orbital moment is roughly proportional to λ e , the charge of e g symmetry in the Fe atomic sphere. We choose the overall phase of the ground-state wave function (1) such that c 0 is real. For the wave function (2), we then find GS|L x |GS = λ e r/(1 + β 2 ), GS|L y |GS = λ e rβ 2 /(1 + β 2 ), and GS|L z |GS = λ e β/(1 + β 2 ) in agreement with Ref. [11].
Our task is now to diagonalize the total Hamiltonian in the subspace spanned by the three states |A 2g , |E 1 g , and |E 2 g . We can write where H 0 is the multiplet Hamiltonian given by where E(A 2g ) and E(E g ) are the multiplet levels. Defining = E(E g ) − E(A 2g ) and taking the reference energy E(E g ) = 0, we have E(A 2g ) = − . The spin-orbit Hamiltonian is where the sum over i runs over all the electrons and Here, α is the fine structure constant, V is the Fe atomic potential, and λ is the charge inside the Fe atomic sphere for a given molecular orbital. Because a 2g corresponds to l z = 0, we only need ζ for the e g orbital, for which the numerical calculation gives ζ = 0.04 eV. Finally, we introduce an external perturbation H ext , which models the effect of the molecular stacking in the film. Stacking gives rise to weak band formation [5] and, because the stack has low symmetry, to a hybridization between the D 4h multiplets.
The parameters b and c are chosen to be real without loss of generality because all wave functions may be taken real before spin-orbit coupling is introduced.
When diagonalizing the Hamiltonian (3) we have to take into account that the sample was magnetized at saturation with a magnetic field (μ 0 H = 5 T) parallel to light incidence, along a direction z rotated by θ , φ from the z axis of the molecular frame (coinciding with the surface normal to the film). We define a spin frame with axis x ,y ,z such that the two frames are congruent when z ≡ z . To describe this situation, we should take the ± spin states in Table III as eigenstates of the spin operator S z in the spin frame, given by (see Sec. III B of Ref. [14]) so that where r α is the component of unit vector of the spin direction (i.e., cos φ sin θ , sin φ sin θ , cos θ ). The diagonalization of the total Hamiltonian (3) then leads to the following eigenvalue problem: where ω is the energy eigenvalue and r the anisotropy factor introduced in Table IV. Here and henceforth, all energies are in units of ζ /2 = 0.02 eV. The eigenvalue equation is where we have put To find the eigenvector coefficients, we solve for the last two rows in Eq. (9), using Eq. (10) to equate ω 2 − cos 2 θ to ω (t 2 + r 2 sin 2 θ )/(ω + ). Then, indicating by ω g (θ ) the lowest (negative) eigenvalue, we find for the ground state where is the normalization factor. We have added the subscript (θ,φ) to the ket |GS (θ,φ) in order to stress the angular dependence of the mixing coefficients, which have been indicated by c i (θ,φ). Notice, however, that ω g (θ ) depends only on θ . This procedure introduces an angle dependence into the expression of the ground state, since different directions of the applied field give rise to different states of the sample with slightly different energies. All these angle dependent states are very similar to the state in Eq. (2) hypothesized in Ref. [11] as the magnetic ground state of the FePc molecule in the film. Indeed, the GS given in Eq. (2) was expressed in the molecular frame, ignoring the small energy differences introduced by the external field applied at different angles and can be formally recovered from Eq. (11) by putting θ = φ = 0. We defer the discussion on their relation until the end of this section.
In the following, for typographical simplicity, we shall indicate the angular dependence by only the subscript θ , having in mind that a final average over the φ variable will always be performed before comparison with experiments, to account for the random orientation of the FePc molecules in the plane of the film. When necessary for the argument at hand, the φ dependence will be explicitly resumed.
Using Table IV with the angle dependent coefficients of Eq. (11), we can now write down the averages GS|L α |GS θ (α = x,y,z). Putting for short we find from which we obtain the wanted projection Since ω g (θ ) is negative, this projection, averaged over the angle φ, is positive as a function of θ and depends only on t 2 and . The relation (15) is to be compared with the usual expression used to analyze the angular dependence of XMCD data when light impinges at an angle θ on the sample, being parallel (or antiparallel) to the applied magnetic field, where μ ⊥ = 1 2 ( GS|L x |GS + GS|L y |GS ) θ=π/2 (calculated at θ = π/2) and μ z = GS|L z |GS θ=0 (calculated at θ = 0). Macroscopically, assuming a diagonal susceptibility tensor appropriate to our case, Eq. (16) follows from the relation = μ x sin 2 θ cos 2 φ + μ y sin 2 θ sin 2 φ+ μ z cos 2 θ, (17) where μ α is the momentum measured when the magnetic field H lies along α = x,y,z. By averaging over the angle φ one obtains Eq. (16). In contrast, Eq. (15) is the correct quantum mechanical expression for the class of ground states considered here.
The comparison of the experimental values [9] μ z = 0.297 ± 0.007 and μ ⊥ = 0.532 ± 0.004 μ B with the corresponding expressions obtained via Eq. (14), determines the best fit values t = 70 meV and = −60 meV. As anticipated, the calculated values for the axial and planar magnetization depend only on t 2 and . Figure 2 shows the behavior of GS|L ·Ĥ |GS θ given by Eq. (15) as a function of the angle θ in the range [0 : π/2], for t = 70 meV and = −60 meV (labeled gs) in contrast with the same function calculated for the value = 60 meV, after performing an azimuthal φ average over the range [0 : 2π ]. In both cases, we give also the curve described by Eq. (16) 155139-5  (16) ("std") and experimental data ("exp") from Ref. [3].
Here the experimental error bars hardly exceed the symbol size. Parameter values are given in meV.
(labeled std). Also shown for comparison are the experimental points with their error bars from Ref. [3]. We note that there seems to be a slightly better agreement with the experimental data of the curve labeled "gs" for = −60 meV than for the one labeled "std." We defer, however, the discussion on this point until Sec. VI, in the framework of a more realistic treatment of the spin-orbit interaction.
Before finishing this section, we discuss the relation between the ground-state (2) used in Ref. [11] to the one given by Eq. (11) calculated at θ = φ = 0. By equating the real and imaginary part of the corresponding coefficients one is able to get the same result and phase, provided ω g = ±1, which is clearly incompatible with the eigenvalue equation (10) when t 2 = b 2 + c 2 = 0. This result is not surprising, since the state (2) represents a sort of average of all the angle dependent states in Eq. (11). In this perspective, the best way to relate the state (2) to the "average" state in Eq. (11) is to equate the value of the magnetic anisotropy μ ⊥ /μ z = (1 + β 2 )/β to the same quantity calculated in the framework of the present model [see Eq. (18)] and find out for which values of the parameters we obtain the value β = 0.5 used in Ref. [11]. Table V  nature of the old approach and that β was the only parameter used, the correspondence with the new set of parameters is remarkable.

Dependence of the magnetization anisotropy on the model parameter
Equations (18) determine the magnetization anisotropy of an FePc molecule under an external perturbation and a saturating magnetic field. According to the last column of Table IV and Eq. (14), each component of the magnetization is proportional to the charge λ e of the e g orbital inside the Fe atomic sphere, while the ratio of the in-plane to the axial magnetization is proportional to r = 2. However, this ratio is also affected substantially by the value of the energy gap and the hybridization parameter t 2 .
We illustrate its dependence on for representative values of the hybridization parameter t. Table V gives the values of μ z and μ ⊥ in units of μ B for various values of for t = 70 meV. It also shows the value of the ratio I − tot /I + tot for future reference in Sec. IV A. We see that μ z increases in going from positive to negative values of , reaching the experimental value at = −60 meV as found above. The type of anisotropy in the range −80 80 meV is always planar. For = −160 meV (not shown in Table), one encounters a transition to an axial type with μ z = 0.58 and μ ⊥ = 0.50 μ B .
The same trend is found for t = 35 meV and t = 0 in the same range of variation of . This latter case is discussed analytically in Sec. II of Ref. [14].
This behavior is an indication that the external perturbation in the film reverses the order of the |A 2g and |E g states. This finding is reasonable, since the external perturbation, which schematizes the effect of the molecular stack onto a single molecule, besides coupling the states, might affect differently their energy position. In the film in fact, the 3d z 2 molecular state, being oriented perpendicularly to the stack [15], hybridizes more strongly than the 3d x(y)z state with the ligand atoms above and below the FePc molecule in consideration. Being an antibonding state, it is pushed up in energy with respect to the isolated molecule, so that the |A 2g state rises in energy in comparison with the |E g state, due to the different occupation of the two molecular orbitals (see Table III). This argument is confirmed by the constrained DFT calculations of Nakamura et al. [5], who find a ground state of type |A 2g for the isolated molecule and of type |E g for the linear chain. Finally, we note that the best fit value of the parameter t (70 meV) is not far from the value of 100 meV that we estimated for the hopping integral between two adjacent molecules in the film just before Eq. (1), indicating an internal consistency of the model. However, the precise relation between the hybridization parameters (b,c) and the environment-molecule interaction still needs to be investigated, but this goes beyond the scope of the present paper, which is focused on the interpretation of the XMCD core-level spectra at the Fe K and L edges and their pertinence to the determination of the ground state of Au-supported FePc films. Nonetheless, the model introduced here is not limited to this system. As clear from its derivation, it can also describe the case of a single monolayer of molecules (or a single molecule) on a substrate, provided again that the perturbation brought about by the substrate could be considered as weak. This case will be the subject of future investigation.

IV. DERIVATION OF Fe K -AND L-EDGE XMCD AND SUM RULES
The possibility of three states lying very close in energy and being mixed by spin-orbit coupling and an external perturbation in the ground state of the FePc molecule is not new and has been suggested in previous publications [8,16,17]. It is possible to discriminate among the various suggestions by analyzing their prediction with respect to the XMCD measurements at the Fe K and L edge [3,11], obtained by taking the ground-state averages of the orbital moment and the effective spin operator and comparing them with the sum rule data of the experiment. However this would imply the determination of the ground state for each new combination of component states along the lines of Sec. III, a difficult task in the absence of the relevant data. Moreover, the proof of the sum rules in the case of anisotropic systems has never been given. In the following we provide explicitly such a proof.
We start by analyzing the spectroscopic predictions of our ground state (1), deferring the discussion on the alternative suggestions to Sec. V. The possible final states are obtained by exciting the relevant core electron (with both spin values) up to an empty orbital state. For example, in calculating quadrupole absorption and dichroism from the Fe K shell, we need to consider the following final states: where we have omitted the d ± xy states, since they are always occupied. Similarly for a 2p hole. These states are assumed to be good approximate eigenstates of the FePc molecule in the stack, due to its strong ligand-field regime.

A. Quadrupole absorption at the Fe K edge and angular dependence of the XMCD signal
Due to the experimental setting, in the spin frame, we have Therefore, using complex spherical harmonics for the transition operator, we have Here, the notation (1,1) refers to the orbital projection of the transition operator Y 21 (r ) in the two matrix elements of Eq. (21). Notice that also the molecular states are referred to the spin frame. Since the state of the molecule is expressed with reference to the laboratory frame x,y,z, we have to transform everything into this frame. For the transition operator Y 21 (r ), this is accomplished by the transformation [18,19] where D 2 μμ (α,β,γ ) = e −μα d 2 μμ (β) e −μ γ is the Wigner rotation matrix for the rotation specified by the Euler angles α,β,γ that takes the frame x ,y ,z into x,y,z. Therefore where now the transitions I (μμ ) are referred to the molecular frame. Similarly for I (−1,−1). We now specify α = φ, β = θ, γ = 0 and observe that in the expression (23) Remembering that d 2 11 (cos θ ) = (cos 2 θ + cos θ − 1)/2, we easily find that This is the sought relation between dichroism expressed in the two frames. The angular dependence of quadrupole XMCD is obtained by using the expression (11) for the ground state. Writing the complex spherical harmonics Y 21 (r) in terms of real ones we find for absorption where A = 4π 2h ωαk 2 /3 and R 02 = r 4 R 1s (r) R 3d (r) dr is the quadrupole radial matrix element. Therefore the quadrupole dichroic signal in the experimental setting is the last step following from the last column of Table IV. The signal is negative for all θ 's, in agreement with experiments. From Eq. (A7) of Appendix, we derive that the total absorption at the K edge is given by where n h (θ ) is the effective number of the 3d holes, so that Eq. (27) can be rewritten as Figure 3 plots the negative of Eq. (29), averaged over φ, as a function of θ in the interval [0 : π/2] compared with the standard expression GS|L z |GS θ=0 cos 2 θ , where GS|L z |GS θ=0 is given by the first Eq. (14) calculated at θ = 0. Also plotted are the measured integrated intensities of the dichroic signal with their error bars. Again, there seems to be a better agreement with the data of the curve labeled "gs." We show below that the right-hand side (rhs) of Eq. (29) calculated at θ = 0 is also the negative of the ground-state expectation value of the operator 3L z + 2O zzz , introduced by Carra et al. for quadrupole transitions [10].
Using the above expressions, we are now in a position to calculate the intensity ratio for normal absorption between the two helicities and compare with the experimental value. From Eqs. (27) and (11) taken at θ = 0, we find Likewise, remembering Eqs. (13) and (12) calculated at θ = 0, for the other two possible transitions showing only absorption, we find Collecting all these terms we find within a common factor so that the total intensity for left circularly polarized light I − tot is greater than for right circularly polarized light I + tot , in agreement with the experimental finding [11], if one remembers that ω g is negative. From Table V, we see that for the optimal values of t = 70 and = −60 meV the ratio I − tot /I + tot is 2.25, higher than the measured value 1.7(1) [11] by 30%. We ascribe this discrepancy to our neglect of the band formation in the film, which emphasizes the localized character of the excitations (see also the discussion following Fig. 4).

Quadrupole sum rules
According to the quadrupole sum rules established by Carra et al. [10], the total intensity of the quadrupole XMCD at the K edge for a final l = 2 transition in normal incidence is given by 2 5l(2l + 1) GS|L z |GS + 2 2 5l(l − 1)(2l − 1)(2l + 1) The operator O zzz is defined as where the sum over i is over all electrons in the system. Using Table I of Sec. I of Ref. [14], we find that, when l z acts on the manifold spanned by the states |A 2g ,|E 1 g , and |E 2 g with orbital occupation given in Table III [ so that Consequently, in agreement with the dichroic result (27). We now consider absorption and dichroism in grazing incidence. In order to take into account the random orientation of the various grains of the film, we need to calculate the matrix element Having a look at the possible final states in (19) we see that no transition is possible involving both the real and imaginary part of the transition operator (38). Therefore there is absorption but no dichroism, in agreement with the experimental finding. This result is in keeping with the fact that the expectation value of Carra's operator both in the x and y directions is zero. In fact, in grazing incidence we need to evaluate the averages of O ααα for α = x,y. Using again Table I of Sec. I of Ref. [14] we find (see note [20] adapted to this case) so that Consequently, as it should be, since in grazing incidence there is absorption but no dichroism.

B. Absorption from L edges
We now use expression (11) for the ground state to calculate the absorption coefficient and the dipolar XMCD at the Fe L edges. In analogy with Eq. (21), we should calculate where the sum over j z is over the components of the appropriate manifold of L 2 or L 3 states and |F i is one of the possible final states listed in (19) with an 1s hole replaced by a 2p hole. In the molecular frame, this expression becomes where D 1 mm (α,β,γ ) is the usual Wigner rotation matrix. In terms of polar coordinates (α = φ, β = θ, γ = 0), the corresponding dichroic signal is given by [21] We should therefore calculate the quantities (46)

Orbital sum rule
According to the first sum rule, integration over all possible final-state transitions and summation over the two L edges should provide a measure of the orbital moment of the system along the direction θ . We then take advantage of the fact that summing over the two edges is equivalent to summing over all initial p α σ states (α = x,y,z) and replace in Eq. (45) j z with α σ . Moreover, in performing the sum over spin components, we can use the spin representation of the molecular frame since χ m |χ m = χ m |χ m = δ mm . We already used this property for the quadrupole transitions.
The calculation of the various quantities I L (m,m ) is straightforward, but long and tedious, and is given in Sec. III A of Ref. [14]. We state here only the results.
The axial part I L (1,1) − I L (−1,−1) of Eq. (46) gives whereas the contribution of the transverse part, proportional to the second step following from the last column of Table IV. Putting together Eqs. (46)-(48), we obtain for the dichroic signal along the applied field: Remembering that, from Eq. (A9) of Appendix, the total absorption from both L edges is where n h is the effective number of 3d holes, we find which is clearly the projection of the orbital moment along the applied magnetic field. We notice, en passant, that the expression for n h found in Appendix, when calculated with the actual values of the effective muffintin charges given in Table II, provides the value of ≈2.66, almost constant in the whole θ range [0 − π/2], in very good agreement with the value used by the authors of Ref. [3] in their analysis of XMCD at the Fe L edge and taken from Ref. [22].

Spin sum rule
In order to find the expression for the spin sum rule, we need to calculate the integrated dichroism from one of the two spin-orbit split L edges. For the sake of convenience, we choose the L 2 edge. Therefore the sum in Eq. (45) should be performed from the 2p 1 2 core state and the corresponding dichroic signal I L 2 is obtained from Eq. (46). The calculation is laborious and is carried out in Sec. III B of Ref. [14]. Clearly, I L 3 = I L 3 +L 2 − I L 2 , where I L 3 +L 2 is obtained from the orbital sum rule (49). Since we need the combination I L 3 − 2 I L 2 to calculate the spin sum rule, we have to evaluate the quantity Therefore omitting for simplicity the θ dependence of the coefficients c(θ ) in the GS, using Eq. (49) and the result for I L 2 derived in Eq. (61) of Sec. III B of Ref. [14], we obtain According to the spin sum rule, the right-hand side of this equation should represent the quantity GS| S eff ·Ĥ |GS| >, where S eff = 2 S + 7 T . Here, S = i s i is the total spin of the system and T = i t i = i ( s i − 3 r i ( r i · s i )/r 2 i ) the total magnetic dipole. Section IV of Ref. [14] shows that this is indeed the case.
By averaging over the angle φ the expression (53), we get [23] GS| S eff ·Ĥ |GS θ exploiting the normalization of the ground state to write |c 1 | 2 + |c 2 | 2 = 1 − |c 0 | 2 . Looking at Eqs. (64) and (65) of Sec. IV of Ref. [14], we can easily make the identification Since all the model parameters have been fixed by the orbital sum rule, we can compare the prediction of Eq. (54) directly with the experimental results (cf. Eq. (6), Fig. 6 and Table I of Ref. [3]). The only quantity needed is |c 0 | 2 , the square of the coefficient of the |A 2g state in the ground state (11). Despite its angular dependence, its value ranges from 0.26 at θ = 0 to 0.33 at θ = π/2. For simplicity of discussion, we fix its value to 0.33 ≈ 1/3 in the whole angular range, without affecting the following conclusions. Notice the change of sign in the definition of the components of m T with respect to those in Ref. [3]. Figure 4 plots the function 2S + 7(m z T cos 2 θ + m (x,y) T sin 2 θ ) for both sets (57) and (56), except that the calculated curve labeled by S eff (gs) has been rescaled by a factor of 0.48 for the reasons discussed below.
There are two main discrepancies between the two results. The first one is the magnitude of the three parameters. This difference can be reduced if one observes that in the spin sum rule we probe the spin structure of the ground state, which was assumed to point along the z axis with the full (non projected) value S z = 1. Since all values in Eq. (56) are proportional to S p z , it makes sense to reduce them according to the ratio between the measured and the calculated value for 2S, that is, 0.64/1.33 = 0.48. In Sec. VI, we shall study the consequences of diagonalizing the spin-orbit Hamiltonian in the whole spin subspace. A further difference concerns the magnitude of the components 7T (z,xy) of the magnetic dipole operator, originating the second discrepancy. We note in fact that the calculated value for the perpendicular dichroic signal is negative, contrary to the experimental finding. The origin of this negative sign can be traced back to the factor of 2 that multiplies the term proportional to λ a of 7T z as opposed to 1 in the similar term of 2S in Eq. (56). This high value of 7T z is related to the particular flat structure of the molecule, which in a crystal field approach has been assumed to model the thin film. As discussed in Sec. II, our model neglects completely the band formation of the thin film along the stacking direction and its spin structure (see, for example, Fig. 4 of Ref. [5] depicting the density of states of the 3d band in the film). In other words, the three-dimensionality of the film might tend to reduce the ratio 7T (z,xy) /2S. We intend to pursue the investigation on this point using a more realistic model for the film.
As a last consideration, we would like to comment on the fact that in our calculation of the XMCD at the L 23 edges, we have neglected the effect of the 2p core hole on the reorganization of the final excited states. For the orbital sum rule, the reorganization is ineffective since we sum over all possible 3d transitions. There might be a difference for the spin sum rule. However, we have shown that the calculated XMCD signal for the particular combination (53) of the individual L 23 edges is equal to the ground-state expectation value of the S eff operator. This equality implies that the final states used in the L 23 transitions are a good approximation of the excited states including the effect of the core hole, due to the strong ligand field regime of the FePc molecule.

V. EXAMINATION OF ALTERNATIVE CANDIDATES FOR THE FePc GROUND STATE
We examine in this section the consequences of the absence of quadrupole dichroism in grazing incidence to discriminate against alternative candidates. In the GS (11), which we now label |GS 1 , this absence is mainly due to the fact that all its SD's have the molecular orbital d xy doubly occupied. This is reasonable, since according to Table II, its orbital energy lies ≈0.3 eV lower than the next one, in other words, it is quite stable. This is also the reason why the 3 B 2g state was found to lie 0.821 eV higher than the 3 A 2g state. We mentioned earlier that this finding is in contrast with the constrained DFT calculation of Ref. [5], where the 3 B 2g state was found to lie at roughly the same energy as the 3 E 1,2 g states. The absence of dichroism in grazing incidence gives us an argument to exclude this possibility.
In fact, if the four states above were within ≈0.1 eV of each other, the perturbation of the stack would also mix in the GS the 3 B 2g state, so that we would have interactions. By rotating the CF parameters they fit the angular dependence of the linear dichroism of the x-ray absorption spectra (XAS) of Ref. [3] and calculate the XMCD at the L 2,3 edges. The CF parameters so obtained determine the relative energies of the molecular orbitals, indicated in square brackets (in eV): Note that these energies are very different from ours and from what is generally found in DFT calculations [5]. The diagonalization of the full many-body Hamiltonian without the spin-orbit interaction gives the E g doublet (a 2 1g ,e 3 g ,b 1 2g ) as the ground state, followed by the B 2g singlet (a 1 1g ,e 4 g ,b 1 2g ) at ∼0.08 eV. These are the same states that in our Table III are labeled E g and B 2g , lying respectively at 1.451 and 0.821 eV. Inclusion of the spin-orbit interaction leads to a ground state (62).
The three states above are to be compared with ours: the singlet A 2g (a 2 1g ,e 2 g ,b 2 2g ) and the doublet E g (a 1 1g ,e 3 g ,b 2 2g ), keeping the same ordering. As one can see, within each set the doublet and the singlet differ by the transfer of one e g orbital to another of type a 1g in the respective Slater determinants. Therefore the matrix elements of L α (α = x,y,z) are the same in both sets. As a consequence all the considerations of Sec. III remain valid and cannot be used to discriminate between the two types of ground state. Moreover, we see that the coefficients c 0 ,c 1 ,c 2 in (62) are the same as those of the ground state in Eq. (11), taking into account that the excitation energy of the singlet with respect to the doublet is almost the same in the two cases [remember that we changed sign to the difference = E(E g ) − E(A 2g ) due to the interaction of the molecule with the surrounding].
In order to discriminate between the two sets, we now consider the spectroscopic predictions of the new ground-state (62) concerning the quadrupole dichroism from the Fe K edge and the dipolar dichroism from the L edge. We can immediately see that it is not possible to discriminate on the basis of the quadrupole dichroic absorption between the two ground states, because their average of Carra's operators are the same, as a consequence of the fact that the matrix elements of the orbital moment are the same. For the same reason the orbital sum rule does not provide any new indication. We find a difference though in the spin sum rule, due to the presence of the newly occupied orbital d xy ∝ Y s 22 and the marked difference between the two Gaunt coefficients (see Sec. VI of Ref. [14]): Instead of calculating the transitions to the excited states of (62), we calculate directly GS n | S eff ·Ĥ |GS n θ . Following the same procedure as in Sec. IV B 2, we find (n = 4) GS n | S eff ·Ĥ |GS n = λ b2 + λ a |c 0 | 2 + λ e (1 − |c 0 | 2 ) where we have introduced λ b2 = 0.84, the charge of the orbital d xy (see Table II). Therefore the quantities corresponding to Eq. (55) and (56) are 2S = λ b2 + λ a |c 0 | 2 + λ e (1 − |c 0 | 2 ) = 1.51, using again |c 0 | 2 ≈ 1 3 . The curve labeled S eff (alt) in Fig. 4 represents the function in Eq. (65) rescaled by the factor 0.42 = 0.64/1.51. Its behavior is opposite to the experimental one. The perpendicular dichroism is very strong and even putting 7T z = 0 it cannot be less than the value 2S = 0.64. The experimental result instead seems to indicate a cancellation between the 2S and 7T z terms, pointing to a negative value of the perpendicular component of the magnetic dipole. For the ground state (62), 7T z is dominated by the contribution of the d xy orbital, which is strong and positive [see Eq. (66)] and the only way to change this situation, even in the film, would be an almost equal occupancy of the two spin components of the orbital.
The list of |GS i (i = 1,2,3,4) considered above exhausts all the possible ground states of the FePc film that one can form out of the list (I), due to the fact that on the basis of the large unquenched orbital moment found in the XMCD of the Fe L 23 edges [3] the presence of two degenerate state of E g symmetry is mandatory in the ground state. |GS 2 and |GS 3 are ruled out on the basis of the quadrupole XMCD, whereas |GS 4 leads to a spin sum rule in contrast with the experimental observations [3].
The relevant, nonzero, matrix elements of the spin-orbit Hamiltonian (5), assuming that the quantization axis is along the direction of the applied field at angles θ,φ, again in units of ζ /2, are 1|H so |2 = − 7|H so |8 = i r cos φ sin θ, To these, we have to add the diagonal matrix elements 3,5,6,8,9) = E E 1 g ,E 2 g = 0.0, (69) and those of the external perturbation H ext in the form of a hopping integral between states A 2g , E 1 g , and E 2 g . Due to its crystal field nature, we assume that they are independent of the spin. Therefore we put 1|H ext |2 = 4|H ext |5 = 7|H ext |8 = −t 1 , Moreover, we add a Zeeman term (ZT) −μ B H eff gS z of 14 meV to account for the effect of the external magnetic field and the internal effective exchange field acting on the Fe atom. Its value will not be subject to variation in the following, because it corresponds to the physically measured quantity [24,25]. The resulting matrix is diagonalized for each θ and φ in a suitable mesh. The ground state |GS θφ is given by the following combination: where the basis states |i are those specified in Eq. (67) and the coefficients d i (θ,φ) are those of the eigenvector corresponding to the lowest eigenvalue resulting from the diagonalization of the angle dependent Hamiltonian matrix given by Eq. (68). Section V of Ref. [14] details how the various averages of the relevant operators over the GS (71) can be calculated in terms of the coefficients d i (θ,φ).
We are now able to fit the angular dependence of the projection of the orbital moment along the applied magnetic field at an angle θ from the molecular axis, after averaging over the azimuthal angle φ. Again, the φ-averaged quantities depend only on the combination (t ) 2 = t 2 1 + t 2 2 . The best fit (in units of meV) occurs for t = 63 ; = −50 (72) very near to the values already found in Sec. III (t ≡ t = 70 meV and = −60 meV). Figure 5 pictures the angular dependence of the projection of the orbital moment along the magnetic field together with the standard dependence (16) against the experimental points with their error bars. It reproduces the behavior already found in Fig. 2 of Sec. III. Notice, however, the slight reduction of the value of μ ⊥ from 0.56 to 0.53 in better agreement with experiments. The distribution of the experimental points seems to favor the quantum behavior of L ·Ĥ with respect to the standard one. Figure 6 shows the angular dependence of GS| S eff · H |GS θ labeled "gs," together with the standard expression GS| S eff ·Ĥ |GS θ=0 cos 2 θ + 1 2 GS|S eff x |GS θ=π/2 + GS|S eff y |GS θ=π/2 ) sin 2 θ (73) labeled "std." Also shown is the φ averaged quantity GS|2 S p · H |GS (θφ) , which is almost constant, around 1.20, to be compared with the value of 1.33 found in Eq. (56). The superscript p in S p refers to the fact that, in taking the average of the spin components over |GS , the integration volume is restricted to the central atom, as dictated by the presence of the core hole.
In keeping with this slight reduction of the value of 2S p , we find that the terminal values of GS| S eff ·Ĥ |GS θ for θ = 0,π/2 differ only by 10% from those found in Sec. IV B 2. In fact, mainly S z = 0 components of the basis states (67) appear in the ground state, so that they can reduce the spin projection on the z axis, but do not affect the value of GS|7 T · H |GS .
Therefore the findings summarized in Figs. 5 and 6 show that the completely spin-polarized model of Sec. III is valid, as long as we look only at properties of the orbital moment, whereas to find agreement with the spin properties (like in the spin sum rule) requires a more realistic approach based not only on the single FePc molecule but also on the unidimensional band formation that induces spin down components in the ground state of the film.
We feel, however, that the results of the present model and the constrained DFT calculations of Nakamura et al. [5], who find a ground state of type A 2g for the isolated molecule and of type E g for the linear chain, are sufficient indication for a ground state of the type (11), in view also of the fact that the other alternatives contradict in a way or another the experimental data, as discussed above.

VII. SUMMARY
We have presented a quantum mechanical derivation of the magnetic ground state of FePc molecules in an Au supported thin film, starting from the ground and low-lying excited states of the isolated molecule. In order to identify the Hilbert subspace in which to diagonalize the spin-orbit interaction, we have performed a molecular orbital multiplet calculation and used ab initio calculated orbital hybrization coefficients throughout the analysis of the ground state and spectra. Thus we have relied neither on DFT predictions nor on standard ligand field multiplet calculations, where free atom Coulomb matrix elements must be reduced uniformly by some 30% to account for hybridization effects [8,17].
Both methods have indeed disadvantages. DFT calculations predict a variety of ground states, depending on computational details such as the exchange-correlation potential and the way self-consistency is achieved. Ligand field multiplet calculations ignore instead the anisotropy factors of the various molecular orbitals that influence the Coulomb matrix elements in a different way. In a situation in which many multiplets lie within 1 eV, disregarding such anisotropy may lead to the wrong ordering of the multiplets, as we have numerically checked.
Our method of calculating directly molecular multiplets, using molecular orbitals as basis functions, avoids this problem and can be expected to give a more reliable answer. It is rewarding that the DFT constrained calculations by Nakamura et al. [5], based on charge symmetry constraints, lead to the same answer, indicating the 3 A 2g multiplet as the ground state, followed by the E g multiplet at about 0.09 eV. This energy order has very recently also been confirmed by many-body diffusion Monte Carlo calculations for the free FePc molecule [26].
We have therefore diagonalized the spin-orbit interaction in the subspace of these multiplets, assuming that the sample is saturated along the external magnetic field the direction of which does not coincide with the molecular axis. Moreover, an external perturbation is introduced to mimic the effects of molecular stack in the film. In this way, the ground state becomes angle dependent, as expected, since different directions of the applied field give rise to different states of the sample with slightly different energies. The ensuing eigenvalue problem can be dealt with almost analytically, providing a better insight into the physics of the magnetic ground state compared to a numerical calculation. For example, it becomes clear why one should reverse the order of the 3 A 2g and E g multiplets found in the isolated molecule, due to an external perturbation, in order to obtain the correct value of the measured L z , or which parameters to vary in order to increase or reduce the magnetic anisotropy as described by Table V. This same table shows that is the main parameter for the control of the magnetic anisotropy.
Based on this ground state, we have also been able to calculate analytically the angular dependence of the quadrupole XMCD at the Fe K edge and of the orbital and spin sum rules at the L 23 edges, showing explicitly that they are equal to the expectation value of the orbital moment and the effective spin moment over the assumed ground state. As a result, we justify quantum mechanically the classical angular dependence of the magnetization, relation (16), usually used to analyze and interpret XMCD data. We note, however, that the quantum and classical predictions of the angular dependence are slightly different and that the experimental points seem to favor the quantum behavior.
This analysis highlights the following conclusions. (a) The ground state of the molecule in the film is not the ground state of the isolated molecule. The high value of μ z = 0.3 μ B found experimentally via the orbital sum rule at the Fe L edge requires the inversion of the energies of the singlet with respect to the doublet as compared to the isolated molecule. In fact, the ground state of the isolated molecule has S z = L z = 0 as can be seen by diagonalizing the 9 × 9 problem of Sec. VI in the molecular frame (θ = φ = 0), putting t i = 0 (i = 1,2), = 93 meV and the Zeeman term equal to zero. (b) The various multiplets entering the composition of the molecule in the film all have the d xy molecular orbital doubly occupied. States with a hole in this orbital give predictions in contrast with the spectroscopic data. (c) The interpretation of the XMCD spectroscopic data allows the identification of these multiplets. If the perturbation on the single molecule due to the stack does not disrupt its electronic configuration, this is also the configuration of the ground state of the isolated molecule, obviously with different mixing coefficients. ACKNOWLEDGMENTS C.R.N. would like to thank S. Di Matteo and Ondřej Šipr for helpful discussions. J.B. and F.B. acknowledge the financial support of the Spanish MINEICO MAT2017-83468-R project, Aragonese E-12 17R RASMIA (co-funded by Fondo Social Europeo), and of the European Union FEDER (ES). We acknowledge support by the European Union COST MP1306 EUSPEC project.

APPENDIX: CALCULATION OF TOTAL ABSORPTION
In order to eliminate the radial matrix element from the expression of the various dichroic formulas, we need to calculate the total absorption coefficient for the transition under consideration. For a quadrupole 1s → 3d (E2) transition, we should calculate where α = |m|,c(s) labels the type of transition operator and the states |F i are those listed in Eq. (19). In an atomic system it would be easy to calculate Eq. (A1) by passing to complex spherical harmonics by a unitary tranformation and making use of the Wigner-Eckart theorem. In this case, the result is where n h is the number of 3d holes, C (2) = 4π 5 Y 2m is the Racah transition operator, and (0||C (2) ||2) = 1 is the reduced angular matrix element. We have n h = 10 − n , where n is the expectation value of the number operator in the l = 2 shell. Moreover, R 02 = ∞ 0 r 4 R 0 (r) R 2 (r) dr with obvious meaning of the symbols.
In the case of an anisotropic system, like the FePc molecule or film that we are considering in this paper, we cannot close the sum over the square of the Gaunt coefficients since the radial matrix elements depend in general on the final state through the transition operator. We have therefore to consider all the possible transitions individually and sum up all the contributions.
If in Eq. (A7), we assume that λ b1 = λ e = λ a = 1, then n h (θ ) = 4, as it should be. From the value c 2 0 (0) = 0.33 and the numerical values of the effective charges given in Table II, we calculate n h (0) = 2.66, in very good agreement with the number of holes used in Ref. [3]. From Eq. (A8), we see that n h (θ ) is angle dependent, but in the interval (0 : π/2) it is almost constant, varying only between 2.66 and 2.64.
A similar, but less straightforward, calculation for the dipolar (E1) transition 2p → 3d leads to the result