Large effect of the metal substrate on the magnetic anisotropy of Co on hexagonal Boron Nitride

We combine x-ray absorption spectroscopy (XAS), x-ray magnetic circular dichroism (XMCD) and x-ray magnetic linear dichroism (XMLD) data with first principles density functional theory (DFT) calculations and a multiorbital many body Hamiltonian approach to understand the electronic and magnetic properties of Co atoms adsorbed on h-BN/Ru(0001) and h-BN/Ir(111). The XAS line shape reveals, for both substrates, an electronic configuration close to $3d^8$, corresponding to a spin $S = 1$. Magnetic field dependent XMCD data show large (14 meV) out-of-plane anisotropy on h-BN/Ru(0001), while it is almost isotropic (tens of $\mu$eV) on h-BN/Ir(111). XMLD data together with both DFT calculations and the results of the multiorbital Hubbard model suggest that the dissimilar magnetic anisotropy originates from different Co adsorption sites, namely atop N on h-BN/Ru(0001) and 6-fold hollow on h-BN/Ir(111).


Introduction
The stability of a magnetic moment against fluctuations critically depends on the magnetic anisotropy energy (MAE). At the microscopic scale, a large MAE may result from large spin-orbit coupling and suitable crystal field symmetry [1,2]. Therefore, although spin-orbit coupling in late 3d-transition metal atoms is relatively weak, under special circumstances the 3d adsorbates can show very large MAE values, as it is the case of Co single atoms and nanostructures on Pt single crystals [3,4]. However, large MAE do not guarantee stability, as it is the case of Co atoms on Pt(111), which have sub-nanosecond spin lifetime and do not show hysteresis even at sub-Kelvin temperatures [5]. These observations highlight that spin scattering with conduction electrons and substrate phonons is detrimental to the stability of the magnetic moment of isolated impurities. Furthermore, atomic-scale magnetic moments can be quenched when interacting with the itinerant electrons on a nearby metal, leading to the formation of Kondo singlets [6]. Thus, large MAE and weak electronic and phonon interactions with the surroundings are key ingredients to achieve stable nanomagnets [1,2,7]. These considerations explain the interest in 3d and 4f atoms deposited on surfaces like MgO/Ag(100) [1,8,9], graphene [10][11][12] or h-BN/Rh(111) [13][14][15]. In these systems, the decoupling layer is introduced to reduce the hybridization of the adsorbate with the metallic substrate, thus protecting the magnetic moments from destabilizing scattering processes, while still allowing for large values of the MAE. For instance, in the case of h-BN/Rh(111) [14], CoH magnetic impurities have been found to adsorb on atop N sites of the h-BN lattice and to show relatively large out-of-plane MAE close to 5 meV with spin quantum number S=1.
The MAE of single magnetic impurities adsorbed on surfaces can be measured mainly by two experimental techniques: x-ray magnetic circular dichroism (XMCD) and inelastic electron tunnelling spectroscopy (IETS). The interpretation of XMCD data relies on electronic multiplet calculations and the sum rules [16,17] analysis with some fitting parameters required to determine the crystal field and spin-orbit coupling strength. The outcome are the spin and orbital angular moments together with the energy spectrum of ground and excited states. The other technique, i. e., IETS [18], yields the spin excitation energies of a single or a few magnetic coupled atoms with atomic resolution. Then, by using a phenomenological spin Hamiltonian compatible with the symmetry of the crystal environment, one can usually describe the magnetic anisotropy in terms of some characteristic anisotropy parameters [1,8,14,19,20].
From the theoretical point of view, a recently developed approach [21,22] consists in using first principles density functional theory (DFT) calculations to obtain the crystal field and, later on, with this information at hand, to construct a crystal field term in a many-body multiorbital Hamiltonian that also includes spin-orbit and electron correlation effects. Additionally, DFT calculations including spin-orbit coupling permit an estimation of the MAE but, in general, the so obtained values for magnetic adatom impurities adsorbed on surfaces are significantly lower than measured values due to the overestimation of orbital momentum quenching [23][24][25].
In this work, we combine XAS, XMCD and XMLD measurements with a multiorbital Hamiltonian model and DFT calculations to disentangle the magnetic properties of Co atoms adsorbed on h-BN/Ir(111) and h-BN/Ru(0001) surfaces. h-BN is an insulating layer with a large band gap of about 6 eV, which can efficiently decouple the impurity from the metallic substrate. The substrate choice is dictated by the different binding energy between h-BN and the two single crystals [26], namely strong on Ru(0001) and weak on Ir(111). This is expected to promote different crystal fields (CF) and, thus, also different magnetic properties for the adsorbates. In agreement with these expectations, line shape and magnetic field dependence of XAS and XMCD measurements combined with multiplet calculations reveal a large out-of-plane anisotropy of about 14 meV for Co on h-BN/Ru(0001), contrary to the low value of a few tens of µeV for Co on h-BN/Ir(111). Additional XMLD measurements at high and low external magnetic fields suggest that the dissimilar magnetic anisotropy originates from different Co adsorption sites, namely atop N and hollow on h-BN/Ru(0001) and h-BN/Ir(111), respectively. The spin and orbital magnetic moment of the Co atom, as well as the MAE, obtained from DFT calculations confirm this picture at a qualitative level. Moreover, the use of a many-body multiorbital Hamiltonian permits the achievement of a better quantitative agreement with experiment and a deeper understanding of two systems, in particular of their response to an external magnetic field.

Sample Preparation
The Ir(111) and Ru(0001) single crystals were prepared in-situ by repeated cycles of Ar + sputtering and annealing. Subsequently, h-BN was grown by chemical vapor deposition (CVD) using borazine (125 Langmuir at 1030 K). The reaction is selflimited to one monolayer since the catalytic dissociation of the precursor molecule requires bare metal areas. Co was deposited from a high purity rod (99.998 %) using an e-beam evaporator in a background pressure of ≤ 3×10 −11 mbar. The Co coverage is expressed in monolayers (ML), where 1.0 ML is defined as one Co atom per h-BN unit cell.

DFT calculations
We have used the VASP package [27][28][29] to perform ab initio calculations based on DFT within the Projected Augmented Wave (PAW) method [30] and the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [31]. The calculations have been done using a 4x4 h-BN supercell, which matches very well with a rotation of 13.9 • of the Ir(111) and Ru(0001) surfaces. In this way, we avoid strain while having a computationally accessible unit cell. The effective vacuum region was larger than 7Å. A Γ-centred 3×3×1 k -point sampling and an high energy (600 eV ) cutoff were used. For the spin polarized calculations we also have included an on-site Coulomb interaction (U = 4 eV ) to account for electron correlations in the Co 3d shell [32].
Concerning the DFT calculations including spin orbit coupling (SOC), they have been done using the tetrahedron smearing method with Blöchl correction [33]. DFT+SOC calculations can be used to extract information about the magnetic anisotropy of the system. This requires a very strict convergence, which was obtained by using an energy convergence threshold smaller than 10 −6 eV. We achieve this convergence using Γ-centred 5×5×1, 7×7×1 and 9×9×1 k -point samplings.

XAS, XMCD and XMLD data
We investigate the electronic and magnetic properties of Co monomers created by low-temperature deposition on h-BN single layers grown by CVD on Ru(0001) and on Ir(111) (see methods for the sample preparation). Figure 1 compares the XAS spectra of isolated Co atoms deposited on the two surfaces with the spectra obtained from multiplet calculations with two different approaches. The experimental data show a fine multipeak structure, fingerprint of the Co electronic state. In particular we note that the L 3 edge splits in a main peak at 777.0 eV and in a small shoulder at 778.8 eV. This multipeak structure compared to the broad L 3 shape observed for Co atoms in bulk is the signature of an electronic state partially preserving an atomic like-character. In order to disentangle the Co electronic state, we performed multiplet calculations using the CTM4XAS6 [34] and the multiX [35] code. In both cases, a CF with C 3v symmetry is assumed, corresponding to Co adsorption on top of the N atom or in the hollow site. The CTM4XAS6 code describes the CF via Dq, Ds and Dt terms and allows for mixed configurations; the multiX code describes the CF via a point charge approach and it does not include partial charge transfer to the ligand, i.e., only pure 3d n configurations are considered. Both calculations represent the data well and reveal a ground state configuration with mainly 3d 8 character on both surfaces. This predisposition to acquire an extra electron in the 3d shell is frequently observed for Co atoms adsorbed on surfaces with low electron density at the Fermi level such as alkali metals [36], graphene [10,11,37] or MgO [1].
Insights into the magnetic ground state and magnetic anisotropy are obtained by measuring XMCD and XMLD spectra at normal (θ = 0 • ) and grazing (θ = 60 • ) incidence as shown in Figure 2. The quite similar XMCD shape and intensity observed at normal and grazing incidence for Co/h-BN/Ir(111) indicate negligible magnetic anisotropy. On the contrary, in the case of Co/h-BN/Ru(0001) we observe a strong angular dependence of the XMCD signal, larger at normal than at grazing incidence, indicating a strong out-of-plane magnetic anisotropy. These conclusions are also supported by the magnetic field dependent (independent) XMLD signal observed in Co/h-BN/Ir(111) (Co/h-BN/Ru(0001)). The microscopic origin of the magnetic anisotropy is the combined effect of the anisotropy in the atomic orbital moment dictated by the CF and the spin-orbit interaction [38,39]. In solids, the orientation of the orbital moment is defined by the CF symmetry and strength. However, in an external magnetic field, S and L tend to align to the field itself; thus, the resulting configuration depends on the competition between CF and magnetic field. Linear dichroism is a measure of the charge density involved in perpendicular versus in-plane bonds forming between Co atoms and supporting substrate. Its field dependence, i.e. the XMLD, thus provides a measure of this competition. In systems with strong CF and magnetic anisotropy, the application of an external field can only marginally change the orientation of L and S, resulting in a field independent XMLD as observed in Co/h-BN/Ru(0001); the opposite behaviour is observed in low magnetic anisotropy systems. The angular dependence of the magnetization curves shown in Figure 3 fully confirm the low and high magnetic anisotropy scenario for Co/h-BN/Ir(111) and Co/h-BN/Ru(0001), respectively.
Comparison of the experimental data with multiplet calculations allows us to provide a more quantitative analysis. Given the mainly 3d 8 character of the Co atoms on both substrates, we focused on calculations with the multiX code in which a pure 3d 8 electronic state has been assumed. These calculations include the effect of the external magnetic field, finite temperature, incidence angle of x-rays, and crystal field environment of the magnetic atom. The CF is defined by effective point charges which position and intensity were chosen in order to simultaneously fit the shape and intensity of XAS, XMCD, and XMLD spectra as well as the shape of the magnetization curves at the two x-ray incidence angles. The point charge distributions that best reproduce the two systems are sketched in Figure 2, the exact position and charge values are reported in the Appendix. They correspond to a hollow and top adsorption site for Co/h-BN/Ir(111) and Co/h-BN/Ru(0001), respectively. This is very similar to the behaviour of Co on graphene on both substrates [11]. Orbital (m L ) and effective spin magnetic moment (m S+D ), given by the sum of the spin and dipolar term, projected onto the x-ray incidence direction are evaluated by applying the sum rules to the experimental and calculated spectra (see Table 3.1). The orbital moment is relatively large on both samples with values close to 1 µ B suggesting an atomic like behaviour; however, it does not reach the close to free atom values observed for Co deposited on other decoupling layers such MgO [1] or graphene on Ru(0001) [11].
In addition, m L shows a strong angular dependence for Co/h-BN/Ru(0001), with the largest value observed at normal incidence, while there is only a fractional reduction by moving from normal to grazing incidence for Co/h-BN/Ir(111). The different angular dependence of the orbital momentum observed in the two systems explains the high (negligible) magnetic anisotropy observed in Co/h-BN/Ru(0001) (Co/h-BN/Ir(111)), as also highlighted by the angular dependence of the magnetization curves in the two systems (see Figure 3). We can quantify the strength of the MAE by the zero field splitting (ZFS), which is the energy difference between the ground and first excited state. The ZFS values, as deduced from the multiX calculations, are shown in Table 3.1. For Co/h-BN/Ru(0001) we find that the ground state consists of a doublet (S z = 0.98, L z = 1.42) separated by 13.7 meV from the singlet (S z = 0, L z = 0.02). The ZFS of 13.7 meV obtained for Co/h-BN/Ru(0001) is quite large, similar to the values observed for Co on graphene supported by different metallic substrates [10,11], and only lower than the ZFS of about 58 meV reported for Co/MgO/Ag(100) [1]. On the contrary, for Co/h-BN/Ir(111) we find that the doublet (S z = 0.92, L z = 1.13) and the singlet (S z = 0.003, L z = 0.022) are very close in energy with a ZFS of only 70 µeV. We note that on both surfaces the lowest energy states are not pure S z = 0, ±1 and L z states, since they contain admixtures of different spins and orbital moments from electronic levels belonging to multiplets higher in energy. This admixture is also responsible for the energy splitting of the ground state doublet by 96.3 (6.6) µeV for Co/h-BN/Ir(111) (Co/h-BN/Ru(0001)), permitting direct quantum tunnelling of the Large effect of the metal substrate on the magnetic anisotropy of Co on hexagonal Boron Nitride7   Table 1. Orbital (m L ) and effective spin (m S+D ) moments (in µ B ), as well as their ratios, for normal (0 • ) and grazing (60 • ) incidence evaluated by applying the sum rules to the experimental (calculated) spectra assuming a hole number n h = 2. The spin value and the ZFS as deduced from multiplet calculations are also reported.

Normal
Grazing   h-BN layer, rather than atop N. The corresponding adsorption energies differ by a few tenths of eV and both are lower than 1 eV, corresponding to weak chemisorption, where van der Waals forces are playing an important role. Indeed, a precise determination of the adsorption energies and the corresponding adsorption sites is not an easy task. The hollow adsorption site for Co on h-BN is in agreement with the result obtained by Yazyev and Pasquarello [40] who used similar computational methodologies. Therefore, as h-BN is weakly bound on Ir(111), we assume the same Co hollow adsorption site also on h-BN/Ir(111). However, the case of Co adsorption on h-BN/Ru(0001) corresponds to a more complex situation, which is at the limit of DFT reliability in giving the adsorption site. This is due to the fact that the hybridization between the h-BN layer and the Ru(0001) surface makes the h-BN/Ru(0001) surface more reactive than the h-BN/Ir(111) surface. A similar situation was found for molecular adsorption on metal surfaces in the weak chemisorption regime, like CO on Cu(111), where the correct adsorption site is only obtained when hybrid functionals, like, BLYP, are used [41]. In case of h-BN/Ru(0001), additionally, the large size of the Moiré supercell formed by h-BN on Ru(0001) makes the problem not tractable on the supercell of the real Moiré pattern. Therefore, we follow a strategy similar to the one used to study the adsorption energy variation in different Moiré domains for h-BN/Ni(111) with a lattice matched model [40] but using a rotated 4x4 h-BN surface unit cell on Ru(0001) with different N-Ru registries, as described in the Methods section. The calculated variation in the adsorption energy of Co on h-BN/Ru(0001) for different registries is of the order of a few tenths of eV and, therefore, comparable to the difference in adsorption energy calculated by moving the Co atom from atop to hollow adsorption on h-BN. Therefore, we prefer to use the information from XMCD and XMLD data to establish the adsorption site. Indeed, our multiorbital Hubbard model calculations described below are consistent with the observation of a large out-of-plane magnetic anisotropy for Co on h-BN/Ru(0001), when assuming atop N adsorption for the Co atom.
The calculated projected density of states (PDOS) onto 3d states of Co for the two adsorption sites on h-BN is shown in Figure 4. The energy integration of these PDOS curves onto 3d states up to the Fermi level (E = 0) gives values of about 7.8 electrons in the d shell for the Co adatom (see Table A3) in both cases, which approximately correspond to a spin S = 1 localized in the 3d shell of the Co atom, in good agreement with experiment. In both cases, the appearance of relatively sharp peaks is consistent with a weak hybridization between the 3d atom states and the h-BN states with practically a full occupation of the majority spin 3d states and of three minority spin states of Co. Both adsorption sites show weak hybridization of Co 3d states with 2p of first neighbour N and B atoms. This hybridization corresponds to d z 2 of Co with p z of N atom for atop N adsorption site and d yz and d xz of Co with p z of N and B atoms for 6-fold coordinate hollow adsorption site.
Our DFT calculations including spin-orbit interaction permit us to estimate the MAE from the total energy difference of two self-consistent calculations corresponding to two different magnetization directions, in-plane and out-of-plane. The high precision required to obtain well converged MAE values and, thus, reliable results, restricts us in practice to consider only the Co atom adsorbed atop N and on a hollow site of a 4 × 4 h-BN layer without the explicit inclusion of the metal surface atoms underneath. In this way, as described in the Methods section, we can do convergence checks with respect to both k-point sampling and energy cut-off in the plane wave expansion. This approximation is also justified by the fact that the magnetocrystalline anisotropy is predominantly determined by the nearest neighbour atoms that define the crystal field. This extreme is further confirmed by the multiorbital Hubbard model in the next section, where the crystal field term considers only four or six atoms for atop N or hollow adsorption respectively. The well-known overestimation of the orbital momentum quenching by DFT calculations for this type of systems usually translates into an underestimation of the calculated MAE and, thus, only the observed trend is reproduced: a significantly higher MAE value (1.5 meV) with out-of-plane anisotropy for Co adsorption on atop N sites compared to in-plane MAE for hollow site (0.4 meV). It is worth mentioning that the switching of the magnetization orientation from outof-plane to in-plane when going from Co adsorption atop N site to Co adsorption on the hollow site, as well as the reduction of its magnitude by about a factor of three, is also found in our multiorbital Hubbard model. However, the corresponding ZFS values are an order of magnitude larger and thereby in good agreement with XMCD data for the case of Co adsorption on atop N site (Co/h-BN/Ru(0001)) but not in the case of Co adsorption on the hollow site (Co/h-BN/Ir(0001)).

Spin excitation spectra from multiorbital Hubbard model
In order to calculate the low energy excitations of the Co atoms adsorbed on the surface we use a multiorbital Hubbard model. The purpose then is to build a many-body Hamiltonian derived from DFT calculations that accounts for the strong correlations in the system. This approach is adequate to describe spin excitations when charge redistribution and lattice deformation are negligible. Considering the information extracted from the fitting of the XMCD profile to the multiplet results of section 3.1, we concentrate on describing the correlations of the partially-filled Co d-shell electrons Table 2. Summary of spin-moment (m S ), orbital moment (m L ) and zero-field splitting (ZFS) of Co atoms on an h-BN layer obtained from DFT calculations including the spin-orbit interaction. Here HA and EA stand for hard axis and easy axis, respectively.
H Coul refers to the electron-electron Coulomb interaction, H CF is the crystal field contribution while H SO and H Zeem represents the spin-orbit and Zeeman interaction, respectively. The dimensionless λ SO ∈ [0−1] coupling constant is introduced to switch on and off the spin-orbit coupling when needed for the analysis of the results. Electronelectron interaction is accounted for only in the Co 3d-orbitals, while the strength of the Coulomb interaction is determined by the Slater integrals F n (3d) [42] or, similarly, in terms of an effective Hubbard repulsion U [22]. Here, we take the values F 0 = 3.41 eV, F 2 = 0.36 eV and F 4 = 0.13 eV. ‡ The spin-orbit and the Zeeman interaction in the presence of an external magnetic field are evaluated following Refs. [21,22]. The crystal field Hamiltonian is estimated using a point-charge model including only the first N and B neighbouring atoms. Both, charges and positions extracted from the DFT calculations are used for the computation, while the atomic-cloud dependent expectation values r 2 and r 4 are taken from the atomic values [43]. The single particle spin-orbit coupling constant (λ SO = 1) is taken as the atomic value ξ Co = 65.5 meV [43]. Although it is well known that these point charge models do not provide a good quantitative description, especially when covalent bonds are present, they yield the right qualitative behaviour and they correctly reproduce the symmetry of the environment. Quantitative aspects will be discussed hereafter. Figure 5 shows the results of our multiplet calculation with the effective multiorbital model for Co on pristine h-BN on two different adsorption sites: atop N (left panels) and on a hollow site (right panels). In both cases, the total angular momentum and spin quantum numbers are L = 3 and S = 1, as expected for an atomic 3d 8 configuration on the basis of Hund rules, but |L z | = 1 for the ground state. The crystal field splits the ground state multiplet, with degeneracy (2L + 1)(2S + 1) = 21, in a different way on the two sites. For the atop N site, the ground state multiplet is an orbital doublet, with total degeneracy 6. By contrast, the hollow site leads to an orbital singlet ground state. As a result, the spin-orbit coupling induces a qualitatively different behaviour in both systems.
As observed from the magnetic field dependence, for the atop N position, the lowest doublet is formed by the S z = ±1 states, while the excited state corresponds to S z ≈ 0, indicating an out-of-plane easy axis system, in good agreement with the data shown in Figure 3. This situation is reversed for the hollow absorption site, where ‡ Our results do not show a significant variation with the Slater integrals as long as the Hund's rules are satisfied. an out-of-plane hard axis is found. The change in the preferential magnetization axis between the two adsorption sites is also corroborated by the average magnetization along the applied magnetic field, see Figure 5(c) and (d). Interestingly, the orbital moment is drastically affected by the adsorption site: the atop N site leads to a significantly larger orbital moment, a situation reported before in similar systems with very high symmetry [1,44]. The origin of this unquenched orbital moment is the large and almost perfectly axial crystal field created by the underlying N atom. On the contrary, the hollow site corresponds to a much lower point symmetry in which both the spin and orbital components along the field direction are quenched compared to the atop N site.
The calculated average spin and orbital moments for the atop N site are in qualitative agreement with the experimental results for Co/h-BN/Ru(0001), showing a saturation for fields around 3 T, with a ratio between orbital (m L ) and spin moment (m S ), at both normal and grazing, of about 0.54, close to the experimetal values of 0.57 and 0.71 from Table 3.1. In contrast, the XMCD data for Co/h-BN/Ir(111) give a magnetic anisotropy much smaller than the one calculated for Co on the h-BN hollow site. Below we discuss some possible reasons. There are two additional key considerations to be done when comparing the experimental results with the multiorbital Hubbard model. The presence of the underlying metal surface partially reduces the symmetry of the Co environment. Furthermore, the adsorption on different regions of the Moiré pattern can also induce important changes in the crystal field felt by the Co impurities. For adsorption on the atop N site, since the dominant source of the crystal field is the underlying N atom, a corrugation of the h-BN layer is not expected to lead to qualitative changes, with only small changes of the MAE between different Co impurities on the h-BN/Ru(0001) surface, assuming the same N-Ru registry. In contrast, on the hollow site the contributions of the six first neighbours are of similar magnitude which, together with the symmetry and the opposite charges on the B and N atoms, tends to cancel the electrostatic potential on the Co site. This situation can also be anticipated from the smaller energy splittings between the dorbital peaks in the PDOS, see Figure 4. Hence, a small corrugation, charge transfer, or lattice strain can induce important qualitative changes for atoms adsorbed on hollow positions, which might partially explain the apparent discrepancy for the Co hollow site adsorption on h-BN/Ir(111).
A similar argument can also be applied to explain the deviations of the spin and orbital moments for Co/h-BN/Ir(111), especially for normal incidence. The moiré induced deviations from a perfect C 3v -symmetry yield a finite quantum tunnelling splitting between the states with S z = ±1 and, thus, a zero expectation value of S a and L a (a = x, y, z). The higher the deviations, the larger the splitting. The external field can partially recover the saturation values only if the induced Zeeman splitting is larger than the quantum spin tunnelling splitting.
Considering the approximations done and the simplicity of the model used in this section, the agreement with the experiment, see Table 3.1 and Figure 3, is remarkable for the atop N site. Notice that despite the similarities in the multiplet fitting of Sec. 3.1 and the multiorbital model employed here, in the latter case we are not using any fitting parameters, and the only free parameters are the Slater integrals that play a minor role, while all other variables entering into the model are taken as their atomic counterparts.  [43] , which depend only on the spin degrees of freedom.
In the case of a spin S = 1 system, the spin Hamiltonian can be written as whereŜ a is the a-component of the spin operator and D and E are the axial and transverse magnetic anisotropy parameters [(|D| =ZFS for a spin S = 1)]. At finite external field B, the spectrum changes due to the induced Zeeman splitting, with a response characterized by the Landé g-factor tensor. It should be pointed out that, although we may not expect the transverse term in (2) from the C 3v -symmetry of the two adsorption sites considered here, this term can appear due to any symmetry breaking, for instance, by the Ru(0001) underneath or by any surface stress. This situation has also been found for Co/h-BN/Rh(111) [14]. The usual approach is to estimate the parameters in (2) from fitting to experiments. Here, however, we take a radically different approach. We construct such a Hamiltonian from the many-body Hubbard Hamiltonian. If we denote by |E n and E n the eigenvectors and eigenvalues of the many-body Hamiltonian (1), and taking advantage that E n |Ŝ 2 |E n ≈ S(S + 1)δ nn , we can construct H S by projecting the (2S + 1) low energy states into the bases of eigenstates ofŜ 2 ,Ŝ z , i.e., H S = n=1,...,2S+1 E nPS |E n E n |P S , whereP S = m z |S, m z S, m z |. We notice that the quantization axis z may be different from the z-axis taken in the DFT calculations, see inset Figure 4, in which case some rotations may be needed to recover the simple form (2). Table 3 summarizes the parameters found from our analysis. In agreement with the previous results, the Co atop N can be described as an easy axis system with theẑ direction out-of-plane. By contrast, the Co on hollow site corresponds to a hard axis. The energy spectra of H S are also depicted in Figure  5 a) and b) with dashed lines. The accordance between these results and those of the multiorbital Hubbard model indicates the robustness of our assumptions in deriving H S .

Conclusions
XAS, XMLD and XMCD measurements reveal large out-of-plane magnetic anisotropy for Co individual atoms adsorbed on h-BN/Ru(0001), while Co atoms on h-BN/Ir(111) have basically no anisotropy. This surprising finding is explained using a combination of first principles DFT calculations and a multiorbital Hubbard Hamiltonian. The most important result of spin polarized DFT calculations is the determination of the spin S = 1 of Co in agreement with the XAS data. In addition, the adsorption sites for Co on h-BN/Ru(0001) and on h-BN/Ir(111) are identified as atop N and hollow sites, respectively, via the construction of a crystal field Hamiltonian that reproduces correctly the essential trends in magnetic anisotropy. The deduced Co adsorption sites for the two surfaces are further supported by comparison with the results of the multiorbital Hubbard model. This also leads to a large easy axis MAE for the atop N position, and a smaller in-plane MAE with hard axis for the hollow site. Our results illustrate how two different underlying metal surfaces can lead to dramatically different magnetic anisotropy energies. These differences are mainly attributed to the different crystal field of the different adsorption sites, while the retained orbital moments and d -shell filling are similar. as well as the magnetization curves, which display the magnetic field dependence of the XMCD peak intensity. Table A1 summarizes the spatial distribution and the values used for the effective charges on the two surfaces. The values of the spin-orbit coupling and Coulomb interactions were scaled to 95 % and 80 % of the Hartree-Fock values, respectively. The experimental line broadening due to the finite lifetime of the core-hole state was modelled by convolution with a Gaussian of FWHM = 0.5 eV. We describe the crystal field of the Co adatom using an effective point charge model, whose parameters were extracted from spin polarized DFT calculations. The charges and their positions are summarized in Table A3.