Engineering magnetic anisotropy in two-dimensional magnetic materials

Abstract Although magnetism is one of the oldest branches of solid-state physics, studies of nanomagnetism are extremely vigorous in recent years, because of the accelerating miniaturization of magnetic units in spintronics devices, which drives the sizes of the magnetic units down to nanometer scale. In this realm, the magnetic anisotropy is the critical factor because it prevents the random spin reorientation induced by thermal fluctuation. Extensive theoretical and experimental efforts have been made to enhance the magnetic anisotropy of the magnetic nanostructures to promote the stability of the magnetization, for the potential applications at high temperature. In this review, we will take a series of examples to address how the magnetic properties including the magnetic anisotropy can be manipulated, as well as the underlying mechanism associated with the manipulation. Thorough understanding of the magnetism of magnetic nanostructures not only provides guidance for engineering the magnetic properties in experiment, but also predicts promising candidates for applications in spintronics devices.


Introduction
The application of innovative experimental tools and the desire to grow nanostructures and novel materials in a controlled manner have produced an urgent need for highly reliable, robust, and predictive theoretical models and tools to achieve a quantitative understanding of their growth dynamics and the size, shape and process dependences of their physical properties. It has been recognized in many fields of materials science that state-of-the-art ab initio electronic structure calculations based on the density-functional theory [1,2] have been enormously successful, in both explaining existing phenomena and, more importantly, in predicting the properties of new systems. For example, the prediction of enhanced magnetic moments with lowered coordination numbers at clean metal surfaces has stimulated both theoretical and experimental explorations for new magnetic systems and phenomena in man-made nanostructures [3,4]. In addition, the theoretical exploration of the symmetry of the interface structure in the Fe/MgO/ Fe greatly facilitated the realization of tunneling magnetoresistance (TMR) junctions composed of the Fe/MgO/Fe(001) sandwich structure in experiments [5][6][7]. Synergistic applications of theory and experiment, as demonstrated repeatedly in many areas of materials science, become a 'must' to further advance our understanding of magnetism in reduced dimensions.
TMR junctions have been widely used in spintronics devices such as nonvolatile magnetic random access memories. As the spintronics devices keep miniaturizing rapidly in the recent years, TMR junctions with sizes scaling down to a few nanometers are desired [8]. To this end, the thickness of the ferromagnetic thin films in TMR junctions should be reduced to a few atomic layers. It is known that the magnetism of ultra small entities is unstable due to super paramagnetism [9]. Thermal fluctuation is the major obstacle for further size reduction of magnetic devices and has become an urgent issue for fundamental research. Innovative strategies to block thermal fluctuation with high magnetic anisotropy energy (MAE) have attracted enormous attention in recent years. Moreover, it is preferred to have the easy axis perpendicular to the surface of the ferromagnetic thin films, a stipulation for the next generation high-density magnetic recording media [10][11][12][13][14].
For ferromagnetic ultrathin films, MAEs are mainly governed by the magnetocrystalline anisotropy (MCA) which derives from the spin-orbit coupling (SOC). However, the strong magnetic materials such as 3d transition metals (TMs) which are commonly used in TMR junctions usually display weak SOC. Therefore, heavy rare-earth and noble TM elements (such as Gd, Tb, Pd and Pt) were alloyed with 3d TMs (especially Fe and Co), in order to enhance the SOC through the interaction between the light and heavy TM atoms [15][16][17][18]. Interestingly, large perpendicular magnetic anisotropy (PMA) was successfully achieved in a series of alloys as expected. For example, the MAEs in L1 0 FePt and CoPt crystals are, respectively, enhanced to 2.7 and 1.0 meV per formula unit, hundreds of times of the intrinsic MAEs (~10 −3 meV per atom) in bulk Fe and Co [19][20][21]. The MAEs are even larger in ultrathin films, due to the more localized surface states and quantum confinement. It was found that the MAE can be as large as ~5.2 meV/ Fe atom in the sandwich structure Pt/Fe/Pt(001) where the capping Pt layer and ferromagnetic Fe layer are both only one-atom thick [22,23]. For the sake of applications in TMR junctions, it is more attractive to use less expensive but more abundant elements such as Mo and W as substrate. Unfortunately, the Fe/Mo (110) and Fe/W (110) systems show only in-plane MA [24,25]. Meanwhile, alloys and multilayers with only 3d TM elements were also investigated, but their MAEs are much smaller than that of Pt/Fe/Pt(001) [26][27][28]. Besides, many other tactics have been proposed to tailor the magnetic anisotropy of ferromagnetic ultrathin films by manipulating lattice strains, compositions, capping layers, growth procedures, surface adsorbates, and external electric field (EEF) [29][30][31][32][33][34][35][36][37][38][39]. The underlying electronic mechanism has been extensively explored through density functional theory (DFT) calculations [40][41][42][43].
Parallel to the studies on ferromagnetic ultrathin films, magnetic nanostructures, which typically contain tens to thousands of magnetic atoms, have spurred increasing interest, due to both the crucial role that they play in advanced magnetic information-storage devices and the fundamental understanding from the investigation of magnetism at the nanometer scale [44,45]. Generally, the magnetic moments per atom in nanostructures of the 3d transition-metal elements such as Fe, Co, and Ni are typically enhanced compared with those of the corresponding bulk metals [46], whereas 4d elements such as Rh and Pd, which are not magnetic in bulk, can exhibit a magnetic moment in small nanostructures. Research on magnetic nanostructures has progressed rapidly in the past decade, partially because of experimental advances in synthesis and characterization [47,48]. Usually, magnetic nanostructures include magnetic clusters that are artificially fabricated and organic single-molecule magnets that widely exist in nature [49][50][51]. To build spintronics devises, the magnetic nanostructures will be placed on proper substrate, such as Ag(001) or Pt (111) [52,53]. However, conducting substrates have a series of disadvantages, such as the nanostructure is highly mobile, the itinerant electrons compensate the magnetic moment and mediate long-range interaction between the clusters, and strong screening effect. Semiconducting substrates are more favorable to achieve stable and controllable magnetic nanostructures [44]. Great progress on fabricating two-dimensional (2D) semiconducting materials has been made in recent years [54], since the discovery of graphene [55]. Nowadays, it is easy to obtain high quality 2D materials such as graphene, hexagonal boron nitride, and transition metal dichalcogenides, and the thickness of these materials are only one-atom to a few-atom thick. Therefore, they are good platform to investigate the properties especially the magnetic anisotropy of magnetic nanostructures.
In this paper, we will review the current status of theoretical determination for magnetic anisotropy of transition metal thin films and magnetic nanostructures.
Examples used here, mainly from our own work, will elucidate the main challenges, solutions and prospects of these vigorous research fields.

Theoretical background
Generally, magnetic anisotropy depends on two factors: the spin-orbit coupling which contributs magnetocrystalline anisotropy and the magnetostatic dipoledipole interaction which contributes shape anisotropy. The shape anisotropy energy (E shape ) can be evaluated through summing up over discrete lattice Here A is a unit-dependent constant; ⇀ m and ⇀ r are local magnetic moments and position vectors, respectively. For magnetic thin films, the shape anisotropy is important, but for magnetic nanostructures, the shape anisotropy is negligible. The determination of magnetocrystalline anisotropy energy (E MCA ), on the other hand, was a major challenge for theoretical calculations in the last two decades. The reliable determination of E MCA for a given material requires highly accurate electronic structures and proper treatment of the weak SOC Hamiltonian which can be expressed as In practical calculations, the SOC Hamiltonian is expressed in a matrix format in the spin space (the two states are denoted as ↑ and ↓, respectively) as, Here θ and ϕ are for polar and azimuthal angles of magnetization and A + = e i∅ (L x + iL y ) while A − = e i∅ (L x − iL y ). For 3d TMs, it has been shown that involvement of H SOC alters charge density, spin density and spin moment negligibly [41]. In contrast, H SOC has to be invoked in the self-consistence loop for the solution of the Kohn-Sham equations when heavy elements are involved, especially in nanoentities in which 4d and 5d elements are also magnetic.
In low dimensions, the leading term in E MCA is the uniaxial anisotropy, E MCA = Ksin 2 θ + Osin 4 θ, sin 4 ϕ). For the determination of K, the magnetocrystalline anisotropy force theorem [56,57] (1) (2) was adopted in most previous ab initio calculations. Here ɛ i (0 0 ) and ɛ i (90 0 ) are the eigenvalues of the Kohn-Sham equations with the magnetization perpendicular and parallel to the surface, respectively. Here positive and negative E MCA indicate perpendicular and in-plane MCA, respectively. Because the sets of occupied states, i.e. {occ′} and {occ″}, were determined through the Fermi filling scheme, which relies on very limited information from ε i alone, a huge number of k-points (>10,000 in the two dimensional Brillouin zone for thin films) are needed to restrain the numerical fluctuation [56,58]. This hurdle was overcome in the last decade using various broadening techniques [59,60], the state tracking [61], and the torque [62] schemes. Particularly, the torque method circumvents the need to differentiate energies. Instead, E MCA is evaluated through the expectation values of the angular derivative of H SOC with wave function at a magic angle ψ′(θ m , ∅ m ): Here ξ is the strength of SOC; ε u,α and ε o,β are the energy levels of the unoccupied states with spin α (�u, ⟩) and occupied states with spin β (�o, ⟩), respectively. Therefore, the keys to achieve large MCA are: (i) large SOC constant ξ which exists in heavy elements such as 5d TM elements; (ii) specific energy diagram to reduce the denominator in Equation (7) that can be realized by appropriate ligand field.

Transition metal thin films
Magnetic TM thin films were widely investigated in the last two decades, in order to find out good candidates with large perpendicular MCA for their potential applications in spintronics devices. In most of magnetic thin films, large uniaxial MCA is induced not only by the surface and interface effects but also by lattice strains due to lattice mismatches and the presence of step edges [30,70]. The surface and interface contributions sensitively depend on the atomic relaxations, chemisorptions and growth morphologies. Well-optimized structures are hence prerequisites for the determination of E MCA through DFT calculations. Fortunately, most of the modern DFT packages can directly calculate atomic forces to expedite the structural optimization procedures. Nevertheless, one should note that while generalized gradient approximation (GGA) significantly improves results of atomic structures for 3d TMs over local spin density approximation (LSDA), it overestimates the volume of 4d and 5d elements as well as oxides. Extra care is thus needed to handle systems with mixing building blocks.
To understand complex magnetic films, it is instructive to first analyze the electronic origin of the E MCA for a freestanding magnetic monolayer. The k distributions of E MCA and the band structures of a Co monolayer [41], for example, are plotted in Figure 1 along the high-symmetry directions in the 2D Brillouin zone. Obvious correlation between the abrupt changes in E MCA and the locations where bands cross the Fermi level discloses simple physics. Clearly, the large negative E MCA for the Co monolayer (−1.34 meV/atom) mainly stems from contributions around the M point. This can be further traced to SOC interaction(s) between the occupied d xz,yz states (m = ±1, denoted as state 5 in Figure 1) and the unoccupied   Figure 1). Note only the pair across E F with the same (different) magnetic quantum number(s) leads to a positive (negative) contribution to E MCA . By knowing the origin of E MCA of simple systems, one can tailor magnetic anisotropy by engineering the bands of the magnetic layer. For example, owing to the Co-Cu d-band hybridization, the Co-d xz,yz states are lowered in energy in Co/Cu(001) and Cu/Co/Cu(001). As a result, E MCA becomes less negative in Co/Cu(001)and is positive in Cu/Co/Cu(001) [71], in good accordance with experimental data [72]. Many calculations have been done in the last decade to study the effects of metal substrates or capping layers on the E MCA of ultrathin magnetic thin films [41,42]. It is now well established that values of K for magnetic films can be determined through high-quality DFT calculations with satisfactory accuracy; some of the calculations also give atom-resolved contributions from different layers [73].
Interestingly, E MCA of Fe, Co and Ni films can be tuned by surface chemisorptions of O, CO and H [35,[74][75][76][77][78], and the underlying electronic mechanism has been actively explored through theory and experiment interplays. Figure 2 illustrates theoretical and experimental results of MAEs for O/Ni/Cu(001) from different regions, namely, the y interception of K(1/d) line is for bulk and its slope is for surface/interface [43]. For O/Ni/Cu(001), the bulk contribution is 27 meV/ atom, primarily caused by the tetragonal distortion in the Ni films grown on Cu(001) [34,60,79].  Figure 2(a). This brings the optimum record up to date for SRT in Ni/Cu(001), as compared to cases using Cu, hydrogen or CO. Theoretical calculations reproduced the trend of experimental data very well, indicating the appropriateness of models. Analyses in electronic structures furthermore attributed the O-induced change in E MCA to the new surface state with the d xz feature caused by the O adlayer. O-and H-induced relaxation and buckling in the second layer also play a role. Surprisingly, the presence of O does not purge the magnetization in the surface region [43]. Instead, the Ni surface magnetic moments are significantly enhanced (by more than 5% for Ni alone), especially when the sizeable spin polarization of oxygen (m O = 0.18 μ B ) is included. However, the magnetic moments of both O and Ni are significantly reduced if the atomic structure is optimized with LDA that produces better agreement with experimental X-ray magnetic circular dichroism data [80].

Transition metal dimers supported by 2D materials
Dimers represent the small-size end point in the magnetic nanostructures. Recently, it was predicted that several freestanding transition metal dimers (TMDs) such as Rh-Rh, Ir-Ir and Pt-Pt might possess MAEs of 40-70 meV [81][82][83]. However, actual realization of such huge MAEs is extremely difficult without a careful search of suitable supporting materials. Graphene can be a benign substrate since it is easy to grow and highly stable. In fact, Xiao et al. recently reported giant MAEs of Co-Co and Co-Ir dimers on benzene and graphene, indicating the possibility of making usable magnetic nanostructures on graphene [84,85]. Nonetheless, most TM adatoms are highly mobile on pristine graphene [86] and various geometrical configurations of Co-Co dimers are almost degenerate [87]. Hence, the structural instability of TMDs on graphene is a major challenge for room temperature applications, except for the large MAEs.

Transition metal dimers on defected graphene
It has been revealed that point defects such as single vacancies (SVs), divacancies (DVs) or nitrogenized divacancies (NDVs) can be created in graphene controllably in experiments [88][89][90][91]. By placing a TMD on a defect site of graphene in a vertical geometry as depicted in Figure 3(a), one of the TM atoms (A-type) turns away from graphene so it can retain its magnetization, whereas the other one (B-type) is embedded in the carbon network and hence becomes immobile. A series of TMDs were considered [92], using 5d elements for the A-type atom (i.e. A = Os, Ir or Pt) to take the advantage of their strong SOC effect and Fe group and Co group for the B-type atom (i.e. B = Fe, Ru, Os, Co, Rh or Ir) to ensure substantial magnetic moments. Figure 4 shows the calculated MAEs. There are nine cases with positive MAEs larger than 30 meV, namely, Ir-Co@SV, Pt-(Co, Rh, Ir)@SV, Os-(Ru, Os)@NDV, Ir-(Co, Rh)@NDV and Pt-Co@NDV, out of 36 A-B@SV and A-B@NDV structures.
For the practical uses, another major concern is regarding their structural stability. There are several mechanisms that may trigger the instability of the vertical geometry of TMDs on graphene. First, the A and B atoms may switch their positions and lose their large MAE. The probability of such position exchange should scale with the energy difference Table 1 lists the ΔE pos-ex for all cases with large positive MAEs. Clearly, Pt strongly prefers the A site, with large positive ΔE pos-ex (>1 eV), whereas Ir and Os tend to switch with other magnetic atoms on both SV and NDV. This is reasonable since both Ir and Os are more active than 3d and 4d elements toward defected graphene, whereas Pt is an 'inert' element. It is more difficult for the position exchange in  A-B@NDV, due to the strong binding between NDV and the B-type atoms, even for cases that have negative ΔE pos-ex . The other structural instability associates with the drift of the A-type atoms away from the vertical geometry, depending on the competition between metalmetal bonds and metal-carbon bonds. As shown in Figure 5, four possible pathways for the diffusion of A-type atom near SV or NDV were considered. From the energy profiles along the diffusion pathways, there are three cases (Pt-Ir@ SV, Os-Ru@NDV and Os-Os@NDV) which exhibit good structural stability. The energy barrier for Pt diffusion along the path 1 is 0.87 eV, while the energy barrier for Os to diffuse away from Ru/Os through all pathways is larger than 1.6 eV. Therefore, both Pt-Ir@SV and Os-Ru@NDV are promising candidates that have both good structural stability (with large energy barrier for diffusion) and giant positive MAEs over 60 meV.
In order to understand the driving forces for the development of giant MAEs and hence find out possible approach the manipulate the MAE, it is instructive  [92]. notes: For all cases, the total energy of the vertical geometry is set to zero. Figure reproduced from Ref. [92].
to the electronic structures such as the projected density of states (PDOS) as shown in Figure 6. Combined with Equation (7), it was revealed that the couplings between states in the minority spin channel, through ⟨d u xy∕x 2 −y 2 �L z �d o xy∕x 2 −y 2 ⟩ of Pt and ⟨d u x 2 −y 2 �L z �d o xy ⟩ of Os, play the dominant role in the giant positive MAEs of Pr-Ir@SV and Os-Ru@NDV. It is also found that the Coulomb correlations between the 5d orbitals are rather moderate, because they are localized, and well separated in energy and space.
The dependence of the MAEs on the position of the Fermi level can be obtained from Equation (5) by varying the highest level of occupied states. As shown in Figure 7(a), the competition between three spin-resolved components, i.e. uu (from majority spin states), dd (from minority spin states), and ud/du (from across spin states) produces a narrow peak of MAE near the actual Fermi level, E 0 F , for Pt-Ir@SV. This mainly associates with the narrow peak of the Pt − d xy∕x 2 −y 2 states right near the Fermi level. In contrast, the plateau of MAE near E 0 F for Os-Ru@ NDV appears to be mainly from the minority spin channel (dd), particularly from coupling between the d xy and d x 2 −y 2 orbitals of Os. Interestingly, the total MAE of Pt-Ir@SV rapidly decreases to zero when the E F shifts upward, which corresponds to more electron occupancy on the d xy∕x 2 −y 2 orbital of Pt. This provides an easy way to operate the magnetic unit by adding more electrons or by applying an EEF (ε) [92]. In fact, the presence of an electric field toward the graphene sheet noticeably reduces the magnetic moment and MAE of Pt-Ir@SV, as shown in Figure 7(c). For instance, its MAE becomes less than 10 meV at ε > 0.8 V/Å, easy for the writing operation. As shown in Figure 7(d), the main effect of the electric field is to induce charge redistribution from d o z 2 orbital to d u xy∕x 2 −y 2 within the Pt atom. As a result of the additional occupation of the d u xy∕x 2 −y 2 states in the minority spin channel, their contributions to M S and MAE(dd) are reduced. In contrast, the actual Fermi  The torque method in Equation (5) treats the SOC term in the Hamiltonian as a perturbation [61,69]. So the validation should be checked for the heavy 5d elements here, by comparing the MAEs from the torque method to those from self-consistent calculations. In principle, both spin and orbital magnetic moments of heavy atoms should depend on the orientation of magnetization and the easy axis is usually along the direction that maximizes spin and orbital moments [82,83,93]. As listed in Table 2, the spin moment of Pt-Ir@SV changes from 0.60 to 0.44 μ B and its orbital magnetic moment changes from 0.96 to 0.34 μ B , as the magnetization is switched from the easy axis to the hard axis. Large changes of spin (from 1.87 to 1.28 μ B ) and orbital moments (from 1.29 to 0.21 μ B ) are also observed Table 2. Total and atom-resolved spin and orbital magnetic moments (in μ B ) of Pt-ir@sv and Os-Ru@ndv with the direction of magnetization along the vertical (⊥) or in-plane (↔) axis.
notes: in both cases, '⊥' direction is the easy axis. Figure reproduced from Ref. [92]. for Os-Ru@NDV, as expected from its huge positive MAE. Self-consistent calculations give MAEs of 73 meV for Pt-Ir@SV and 54 meV for Os-Ru@NDV. These values are rather close to the corresponding values from the perturbative torque method, indicating the applicability of both approaches for the determination of MAEs. In particular, the torque approach is more efficient and practical for studies of large systems and also for broad searches as in the present work. Finally, it is feasible to fabricate the vertical Pt-Ir@SV and Os-Ru@NDV structures on graphene. Embedding TM atoms in SVs and NDVs has been realized [91,94,95] and the binding energies for Ir to SV and Ru or Os to NDV of graphene are very large (5.3 and 7.3 eV, respectively). Therefore, it should be possible to first trap the B-type atoms (Ir and Ru) on SV and NDV. One extra step to obtain the Pt-Ir@SV and Os-Ru@NDV structures is the deposition of the A-type (Pt and Os) atoms. Pt and Os atoms are reasonably mobile on the perfect graphene, and Ir@SV and Ru@NDV act as attractive centers for them as indicated in Figure 5. Therefore, the formation of vertical Pt-Ir and Os-Ru dimers should be straightforward in a well-controlled lab condition. The actual application of this prediction requires more fundamental and engineering studies, including finding appropriate protecting materials and fabrication conditions. Nevertheless, an area density of 190 TB/in 2 is conceivable for magnetic recording and quantum computing operations, providing that the 8 × 8 graphene cells can be fabricated as the smallest magnetic units.

Transition metal dimers on hexagonal boron nitride
Hexagonal boron nitride (h-BN) has almost the same atomic structure as graphene, so the prediction in Section 4.1 can be directly generalized to h-BN [96]. Although there are two type single vacancies for h-BN, i.e. boron (V B ) and nitrogen (V N ), previous theoretical and experimental studies showed that V B is much easier to form than V N [97,98]. Similar to that in graphene, Fe (Ru, Os) or Co (Rh, Ir) group atom was chosen as the B-type atom to be embedded into V B and a 5d TM atom (i.e. Ta, W, Re, Os, Ir or Pt) as A-type atom was placed over the B-type atom to form the hypothetical A-B@Dh-BN structure. After relaxation, the B-type atoms typically move outwards in a range of 1.3-1.8Å as shown by h 2 in Figure 8(a). Meanwhile, their nearest nitrogen neighbors are pulled out of the h-BN plane (h 1 ) by 0.4-0.8 Å. The distances between A-type and B-type atoms, h 3 , are in a range of 2.0-2.3Å, which are relatively shorter than those for TMDs on benzene and grapheme [84,85,92]. The binding energies of the B-type atoms to Dh-BN are larger than 8.5 eV as shown in Table 3. This means that the B-type atoms are strongly anchored at V B , or in the other words, the B-type atoms can be quickly attracted to the defected sites on h-BN. Here, the binding energies of B-type atoms are defined as where E(Dh-BN) and E(B@Dh-BN) are the total energies of Dh-BN and B@ Dh-BN cases, respectively; E(B) is the total energy of an isolated B-type atom. Figure 9 shows the MAEs calculated by the torque method [Equation (5)]. It can be seen that some systems prefer in-plane magnetization, suggested by the negative MAEs. They are not suitable for applications because their azimuthal MAEs are typical very small [92]. Interestingly, there are twelve systems that have positive MAEs larger than 30 meV. In particular, Ta-Fe, W-Os and Ir-(Co, Rh, Ir)@Dh-BN have MAEs larger than 100 meV. Using the Arrhenius-Neel model and an attempt frequency of 10 10 Hz, a magnetic unit with MAE ~100 meV may hold its magnetization for a few seconds at the liquid nitrogen temperature or years at 10 K. Therefore, these A-B@Dh-BN structures are potential candidates as magnetic units in spintronics devices for information storage and operation.
The possibility to fabricate these hypothetical nanostructures and use them in spintronics devices relies on their structural stability, i.e. whether they can resist dissociation or alternation of the TMDs. While the B-type atoms are tightly pinned at the V B sites, the A-B@Dh-BN system may break down if the A-type atoms drift  Table 3 for cases with large positive MAEs, ΔE A-B are larger than 2.5 eV for all systems except W-Rh, indicating that B@Dh-BN may serve as attractive centers to the A-type atoms. The nitrogen atoms around V B are rather active so there is a competition between metal-metal bonds and metal-nitrogen bonds. Therefore, the A-type atom may drift away from B-type atom if the metal-nitrogen bonds are stronger than the metal-metal bonds. From the energy profiles of the different pathways in Figure 10, it can be seen that the vertical geometry is not favorable for W-Os@ Dh-BN, Ta-Fe@Dh-BN and Ir-Co@Dh-BN since their A-type atoms tend to drift toward adjacent N atom, resulting in lower total energies. On the contrary, Ir-Rh@Dh-BN and Ir-Ir@Dh-BN prefer the vertical geometry and the energy barrier for the A-type atom diffusing out of vicinity of the B-type atom is at least 0.75 eV. Nonetheless, Ir-Rh@Dh-BN may undergo A⇔B switch, and becomes Rh-Ir@Dh-BN with an energy gain of 0.57 eV and MAE of −30.3 meV. Therefore, only Ir-Ir@Dh-BN is a promising candidate for practical use with both large MAE (~126 meV) and high stability of the vertical geometry. Another factor to consider is the Coulomb correlation in TM atoms, which may affect the MAE. It was found that the MAE of Ir-Ir@Dh-BN retains large (>100 meV) in a reasonable range of Hubbard U (≤2.5 eV) for the Ir-5d orbitals.
The PDOS of Ir-5d orbitals and Fermi level dependent MAE are plotted in Figure 11. From Figure 11(a), one can see that Dh-BN has a strong influence on the electronic structure of the B-type Ir atoms, especially for the PDOS near to the  Fermi level. The PDOS peaks of the d xy and d x 2 −y 2 orbitals of the B-type Ir atom in the minority spin channel split and have similar amplitudes around −0.7 eV and slightly above the Fermi level. In contrast, these two orbitals are almost degenerate for the A-type Ir atom and stay near the Fermi level in the minority spin channel. According to Equation (7), the SOC interaction between them leads to large positive MAE(dd) because of the small denominator. Interestingly, it can be seen from Figure 11(b) that the MAE may rapidly reduce to zero when the Fermi level shifts upwards, as a consequence of the full occupation of the d xy and d x 2 −y 2 orbitals of the A-type Ir atom in the minority spin channel. Therefore, the MAE of Ir-Ir@Dh-BN can be conventionally tuned by applying a gate voltage or external electric field to adjust the position of Fermi level. This provides an additional degree of freedom to control its magnetic states required as required for information writing process.

Transition metal dimers on two-dimensional polyphthalocyanine framework
As discussed in Section 4.1, to produce the structure A-B@NDV in graphene, the nitrogenized divacancies should be generated first. Intriguingly, 2D polymeric Fe-Phthalocyanine was synthesized in experiment recently [99]. Here the NDVs already exist uniformly and each NDV contains one Fe atom. This material can be considered as 2D polyphthalocyanine framework (p-Pc) with Fe atoms embedded at the centers of phthalocyanines. Furthermore, the central Fe atoms may be replaced by other TM atoms conveniently [100]. Therefore, the structure A-B@NDV can be achieved by depositing A-type atoms (usually 5d TM atoms) on 2D p-Pc. In our recent paper [101], we studied the electronic and magnetic properties of the homogeneous dimers of all 5d TM atoms on 2D p-Pc (notated as TM 2 @p-Pc), as shown in Figure 12. Similar to the TMDs on graphene and h-BN, the vertical geometry was assumed at first. The B-type atoms are only slightly higher than the p-Pc plane. As listed in Table 4, the N-B bond lengths for all cases are very close, within 2.06 ± 0.07 Å; the A-B bond lengths range from 2.13 to 2.99 Å, due to the large difference of the atomic radii among these TM atoms. All A-type atoms except Au have large spin moments, while those of the B-type atoms are much smaller. Interestingly, the Ir 2 @p-Pc possesses the largest MAE of 47.2 meV MAEs which is about twice of the thermal excitation energy at room temperature (~25 meV). Therefore, it is promising for the applications in spintronics devices at room temperature.
The PDOS of the Ir-5d orbitals in Ir 2 @p-Pc are plotted in Figure 13. It can be seen that d x 2 −y 2 and d xy orbitals of the A-type Ir atom are nearly degenerate and they contribute most part of the spin magnetic moment. By calculating the nonzero angular momentum matrix elements using Equation (7), it was revealed that the dominant contribution to the large positive MAE is from the coupling between the d x 2 −y 2 and d xy orbitals of the A-type Ir atom through L z .
The structural stability was estimated by calculating the energy profiles of the A-type atom diffusing along the pathways ( Figure 12). As shown in Figure 14 for Ta 2 @p-Pc, Os 2 @p-Pc and Ir 2 @p-Pc that have relatively large MAEs, the vertical geometry of Ta 2 @p-Pc is unstable and the A-type Ta atom will tilt away from the B-type Ta atom to reach more energetically favorable geometry. The energy barrier for the A-type Ir atom to diffusion away from the B-type Ir atom is about 0.8 eV,  which implies good structural stability. The Os 2 @p-Pc is very stable and the energy barrier for Os diffusion is larger than 2.5 eV. Interestingly, the MAE of Os 2 @p-Pc is also large, 40.7 meV, which makes Os 2 @p-Pc another promising candidate for the applications in spintronics devices at room temperature.

Two-dimensional metalloporphyrin
Except for building TMDs to produce large MAEs in 2D magnetic materials, we predicted in our recent paper [103] that modifying the functional radicals that do not interact directly with the TM atoms in 2D metaloporphyrin framework (TM@Pp) can also enhance the MAE significantly. The atomic structure of TM@ Pp is shown in Figure 15(a). The central 5d TM atom (Ta, W, Re, Os, or Ir) which provides magnetic moment and magnetic anisotropy of a TM@Pp. Table 5 lists some key properties of TM@Pp. The bond lengths between the TM and N atoms lie in the range of 2.07 ± 0.06 Å, showing strong constraint from the polyporphyrin framework. For an embedded TM atom, its binding with the polyporphyrin framework can be characterized by the adsorption energy where E(TM@Pp), E(Pp) and E(TM) are the total energies of TM@Pp, the pure polyporphyrin framework and an isolated TM, respectively. As listed in Table 5, the binding between the TM atoms and polyporphyrin framework is very strong, with E ad of 10.5~12.4 eV. This is reasonable since four strong TM-nitrogen bonds form through the binding between the TM atom and polyporphyrin framework. In fact, large adsorption energies were also found in similar system -metal mesotetraphenylporphines where the adsorption energies vary from 6.3 to 10.8 eV for (11)  notes: The energy of initial vertical structure is set to zero. The corresponding diffusion path is shown in Figure 1. the late 3d TMs from Fe to Zn [104]. Generally, the 5d TM atoms are more active than their 3d counterparts, which results in larger adsorption energies.
It can be seen from Table 5 that all considered systems are spin polarized and the spin moments are mainly localized on the central TM atoms. The W and Re atoms possess largest local spin moments of 2.7 and 2.8 μ B , respectively, while smallest local spin moment is on the Ir atom, 0.6 μ B . Although the SOC effect    in all TM atoms are strong, only the W, Re and Os atoms induce large magnetic anisotropies when they are placed at the center of polyporphyrin framework, with MAEs of about 24 meV for W@Pp and Re@Pp and −46 meV for Os@Pp. A negative MAE corresponds to the in-plane preferential spin orientation, and the in-plane magnetic anisotropy is very small for Os@Pp (less than 1 meV), which implies that the spin orientation of Os@Pp is undetermined even at very low temperature. On the contrary, a positive MAE indicates the perpendicular preferential spin orientation without additional uncertainty for spin orientation. Therefore, large perpendicular magnetic anisotropy -corresponding to large positive MAE -is beneficial for applications in spintronics devices [26]. Clearly, the MAEs of W@Pp and Re@Pp (~24 meV) are quit large compared to most 3d magnetic nanostructures [105] and close to the thermal excitation energy at room temperature (~26 meV). However, the MAE is still not large enough for applications at room temperature. The magnetic anisotropy of a magnetic nanostructure can be affected by subtle change of the ligand field around the central TM atom, due to the modification of electron occupancy around the chemical potential [92,93]. On the other hand, the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) of organic materials can be manipulated conveniently by adjusting functional radicals with different electron donating capabilities. It has been revealed that the electron-donating capability of some common functional radicals is in the order of: hydrogen (-H) < methyl (-CH 3 ) < hydroxyl (-OH) < amino (-NH 2 ) [106]. Therefore, replacing the H atoms on the edges of the polyporphyrin framework, respectively, by methyl (-CH 3 ), hydroxyl (-OH), and amino (-NH 2 ) radicals, the magnetic properties of the modified TM@Pp are expected to vary with electron-donating capability of the corresponding functional radicals. The modified TM@Pp is denoted as TM@M-Pp, TM@H-Pp, and TM@A-Pp, as shown in Figure 15(b)-(d). After structural relaxation, the TM-N bond lengths and the adsorption energies change little as listed in Table 6.  In most cases, the total spin moments of modified TM@Pp remain nearly unchanged, but those of modified W@Pp change significantly because the electron occupancy increases more in minority spin states than in majority spin states. In addition, the local spin moments on the W and Re atoms decrease notably after hydroxyl and amino modifications. M S,W decreases to 2.5 and 2.3 μ B from 2.7 μ B , respectively; while M S,Re decreases to 2.5 μ B from 2.8 μ B . Interestingly, the MAE is altered substantially for all cases, as shown in Figure 16. For W@Pp and its derivatives, the MAE increases monotonously with the increasing electron donating capability of the functional radicals and reaches 36.7 meV in W@A-Pp, with largest increment by about 50% compared to that in the original W@Pp (24 meV). More strikingly, incorporation of hydroxyl and amino radicals in Re@Pp can enhance the MAEs from 24 to 52 meV and 61 meV, respectively. Such giant MAEs endow the possible applications of these materials at relatively high temperature. Figure 17 shows the PDOS of the Re-5d orbitals in Re@Pp and its derivatives. The states of most 5d components locate near the Fermi level, while the states of d xy orbital locate about 4~5 eV higher than the Fermi level (outside the energy range in Figure 17) due to the repulsion from N 2p orbitals. It can be seen that the majority spin states of d z 2, d x 2 −y 2, d xz/yz are all occupied, while the corresponding minority spin states are unoccupied, resulting in large local magnetic moment on the Re atom. Owing to the geometrical symmetry of the polyporphrin framework, d xz and d yz states are almost doubly degenerated and are much more delocalized. After modification with hydroxyl and amino, the electron occupancy on d xz/yz decreases from 0.80 e to 0.78 e and 0.77 e, manifested by the reduction of the PDOS of d xz/yz as shown in Figure 17. The energy levels of d x 2 −y 2 orbital move towards the Fermi level, leading to reduction of the exchange splitting of the majority and minority spin states. Most strikingly, the minority spin state of the d x 2 −y 2 orbital becomes partially occupied after modification, with electron occupations of about 0.08 e, 0.23 e and 0.19 e for methyl, hydroxyl and amino radicals, respectively. Therefore, the local spin moment on Re atom reduces compared to the unmodified system (see Table 6). Obviously, modifying functional radicals is an effective way to engineer the electron occupancy and hence tailor the MAE. Figure 18 shows the Fermi-level dependent total and decomposed MAEs of Re@Pp and its derivatives. It can be seen that, for both the original and modified Re@Pp, the coupling between cross spin states (ud + du term) dominates from ~−0.5 eV to the Fermi level and is the main origin of the positive MAEs. In contrast, the coupling between minority spin states (dd term) leads to huge negative contribution to the MAEs above the Fermi level, resulting in sharp drop of the MAE above the Fermi level. This feature implies that the MAE can be turned easily from large positive into large negative by injecting extra electron, which offers extra degree of freedom to engineer the MAE of these systems. In addition, the plateau right below the Fermi level of the MAE (ud + du) becomes narrower in width and higher in height as the electron donating capacity of the functional radicals increases, which leads to significant enhancement of the total MAE, especially for hydroxyl and amino. Based on Equation (7), it was found that the contribution from L z [i.e. the first term in Equation (7)] is negligible, while that from L x dominates. The coupling between d xz/yz in majority spin channel and d x 2 −y 2 and d z 2 in minority spin channel through L x [i.e. the sencond term in Equation (7)] contributes to large positive MAE(ud + du) and hence to the total MAE. Furthermore, the energy separation between the corresponding orbitals, i.e. the denominator E u,↓ − E o,↑ in Equation (7), reduces greatly [see Figure 17], resulting in significant increase of MAE.
For W@Pp and its derivatives, the enhancement of MAE by chemical modification of functional radicals is also remarkable. Figure 19 shows the PDOSs of the W-5d orbitals in W@Pp-based systems. By comparing the four systems, the energy level shifts of the W-5d orbitals are distinct. After modification, the energy levels of the d x 2 −y 2 and d z 2 orbitals move downwards, due to the indirect repulsion from the electron donating groups. The magnitudes of the energy level shifts increase with increasing electron donating capacity of the functional radicals. In addition, the substitutional radicals make the d xz/yz band become narrower and move toward the Fermi level. Consequently, the state of the d xz/yz orbital in majority spin channel becomes partially occupied for hydroxyl and amino modified W@Pp systems, leading to significant reduction of the local spin moment on W atom (see Table 6). The large positive MAEs of W@Pp-based systems are contributed from the coupling between d xz/yz in majority spin channel and d x 2 −y 2 and d z 2 in minority spin notes: here, 'uu' , 'dd' and 'ud + du' notate the contributions from the coupling between majority spin states, between minority spin states, and between cross spin states, respectively. E 0 F stands for the natural Fermi level. Figure reproduced from Ref. [103]. channel, the same as the situation in Re@Pp-based systems. However, the states of the d xz/yz orbital in W@Pp-based systems are more delocalized, which reduces the contribution to MAE. Therefore, the enhancement of MAE in W@Pp-based systems is less significant than that in Re@Pp-based systems.
Apart from the magnetic anisotropy, another crucial factor for applications of TM@Pp is the magnetic coupling behavior. To serve as a magnetic building block in spintronics devices, ferromagnetic (FM) ground state is favorable, because antiferromagnetic (AFM) state introduces extra degree of freedom to control [107]. The exchange energy can be calculated as Here E AFM and E FM denote the total energies of AFM and FM states, respectively. As listed in Table 7, W@Pp-based systems prefer AFM coupling, while Re@Pp-based systems prefer FM coupling. Interestingly, the exchange energies of W@Pp series are closely dependent on the electron donating capacity of the functional radicals. As the electron-donating capacity increases, the AFM coupling is gradually strengthened, as indicated by the increasing magnitudes of the exchange energies in Table 7. In contrast, the exchange energies of Re@Pp-based systems are over 180 meV and change little with different functional radicals.  Accordingly, the Curie temperature (T C ) can be estimated using the Ising model and the mean field theory (MFT) [101,108]. Considering the first nearest magnetic coupling, the Hamiltonian can be written as Here, J C is the first nearest magnetic coupling parameter, M i is the magnetic moment at site i. Then we can obtain J C for the TM@Pp as For Re@Pp-based systems, the magnetic moment M i on each site is 3 μ B , and E ex is ~185 meV, so the calculated J C is 1.28 meV. By solving the partition function, the Curie temperature can be expressed as T C = 5γJ C /k B , where γ is the first nearest coordination numbers of the Re atoms and k B is the Boltzmann constant. Hence, the estimated T C is 297 K. It is known that T C is usually overestimated by MFT. More precise Monte Carlo simulations yields T C of about 200 K, suggesting the magnetism can be retained at relatively high temperature.
The feasibility of chemical modifications in experiments was evaluated by calculating the reaction heat (ΔH) for possible reaction routes. Here, the reaction heat is defined as where E(Pro) and E(Rea) refer to the total energies of the products and reactants, respectively. For the proposed reactions in Table 8, the values of ΔH are negative, indicating exothermic reactions. Most importantly, the reaction heat is remarkably large for the substitution of hydrogen by hydroxyl or amino, implying good feasibility for the proposed modification under proper circumstances.

Organic single-molecule magnet on Cu(110) surface
Molecular spintronics has attracted considerable attention in recent years, driven by the accelerating miniaturization of electronic devices [49][50][51]. Metalphthalocyanine (commonly referred to as MPc) molecules are particularly intriguing due to the highly tunable electronic and magnetic properties [109]. With the unusual intermediate triplet spin state (S = 1), FePc molecule is unique among various MPc molecules [110,111]. Tsuhakara et al. reported that the magnetization of FePc molecule reorients from in-plane direction to the perpendicular direction relative to the molecule plane when it is adsorbed on the O-Cu (110) surface [denoted as O-Cu (110)] [112]. It is interesting to reveal the underlying mechanism for the SRT. Furthermore, it is fantastic if the spin orientation can be controlled by conventional approaches. The adsorption of FePc on O-Cu(110) was found neither physisorption nor chemisorption, but with mixed features of both chemisorption and physisorption [93]. The FePc molecule prefers the O-top site with its Fe-N axis along 30° (denoted as the α-type geometry) away from the [001] axis of the Cu lattice [112]. The magnitude of E b (0.36 eV) is small, yet the FePc molecule deforms remarkably and causes a significant surface reconstruction on the O-Cu(110) substrate. The Fe-O bond length is only 1.94 Å, indicating a strong attraction between the two atoms. In contrast, the carbon rings and substrate repel each other, implying that physisorption is more preferable. Therefore, these effects result in the mixed features of chemisorption and physisorption, i.e. a strong ionic bond but with a very small binding energy.
According to Bader's charge analysis scheme [113], the iron atom of FePc transfers 0.43 electrons to the oxygen atom underneath. As a consequence, the spin magnetic moment (M S ) of FePc/O-Cu(110) enhances to 2.40 μ B compared to 2.00 μ B in the freestanding case. This value agrees excellently with the experimental data, 2.30 ± 0.02 μ B [112]. Significant spin-polarization is induced around the O and Cu atoms adjacent to Fe, with M S of 0.15 μ B for O and 0.03 μ B for Cu, respectively. The PDOS for both the free and the O-Cu(110) supported FePc molecule were plotted in Figure 20. Since the freestanding FePc molecule has a D 4h symmetry, the Fe-3d orbitals split into four groups: b 1g (xy) and b 2g (x 2 −y 2 ) for the in-plane components, along with a 1g (z 2 ) and e g (xz and yz) for the out-of-plane components. However, the actual energy spectrum of the Fe-3d orbitals in Figure  20(a) is somewhat different from this simple assignment, because of the strong interaction with N-2p orbitals. For example, the lowest peak in the minority spin channel comprises b 2g , e g and a 1g features all together; the peak right above the Fermi level combines the b 2g and e g components. When the FePc molecule is placed on O-Cu (110), the Fe-a 1g orbital becomes delocalized as shown by the broad PDOS features in Figure 20(b). This manifests the strong hybridization between the Fe-a 1g and O-p z orbitals. Significantly, the two PDOS peaks across the Fermi level become more 'pure': with the b 2g state below E F and the e g state above E F . As a result, the contour plot of the charge density difference shows the intra-molecular charge transfer from the Fe-e g orbital to the Fe-b 2g orbital in Figure 20(c).
The origin of the MAE can be estimated from Equation (7). For convenience of discussion, we present E MCA (dd) in details only because it contributes most part of the total E MCA for both free FePc and FePc/O-Cu(110) as discussed later. The energy diagram of Fe-3d orbitals in the minority spin channel and all non-vanishing L z and L x elements across these states were shown in Figure 21 . As a result, two thin L x lines (i.e. with small L x elements) and two thick L z lines (i.e. with large L z elements) intercept the Fermi level and E MCA (dd) thereby becomes positive. Following this argument, the number of L x lines exceeds that of L z lines if its Fermi level shifts up to above the e g state as seen in Figure 21 Note that E MCA of FePc/O-Cu(110) changes rapidly near E F − E 0 F = 0, with a positive slope as shown in Figure 21(d). This offers an opportunity to tune the magnetic anisotropy of FePc/O-Cu(110) by applying an EEF (ε). Since the screening in the region between the molecule and the substrate is rather weak, the electronic potential around the FePc molecule may easily shift to lower (higher) value with respect to the substrate by a positive (negative) EEF, and so do the PDOS peaks of Fe-3d orbitals. The field dependence of the Fe-e g peak was used as example as shown in Figure 22(a). It appears that the magnitude of the energy shift with ε = −1.0 V/Å is larger than that with ε = 1.0 V/Å. The presence of the EEF alters the electron population of the FePc molecule and also its magnetic moment, as shown in Figure 22(b). At ε = +1.0 V/Å, the Bader charge of FePc molecule becomes −0.66 e, compared to +0.43 e for the zero-field case. When ε is −1.0 V/Å, the Bader charge of FePc molecule is +1.39 e. The inset in Figure  22(a) shows that a negative EEF causes charge depletion from the Fe-e g state to the substrate. The large range of M S in Figure 22  Note that the calculated E MCA (ε) curve in Figure 22(b) closely follows the trend of E MCA (E F − E 0 F ) in Figure 21(d), as predicted by the rigid band model analysis. On the positive side of ε, E MCA first increases to its summit at ε = 0.25 V/Å and then drops gradually afterward. For negative ε, E MCA decreases rapidly and changes its sign near ε = −0.5 V/Å. This is caused by electron depletion from the Fe-e g orbital as well as by the involvement of Cu-d states as shown in the inset of Figure 22(a). This external electric field-induced SRT is very important for magnetic recording and spintronics applications since one has a means to switch the easy axis of FePc/O-Cu (110) between the in-plane and perpendicular direction.
Although the MAEs of free and O-Cu (110) supported FePc molecule are small, the substrate induced SRT and electric field controlled magnetic anisotropy can be generalized to other organic single-molecule magnets which possess large MAEs, to search for potential candidates for the applications in spintronics devices.

Conclusions and outlook
As discussed through several examples in the preceding text, high-quality density-functional calculations can provide reliable results for magnetocrystalline anisotropy in transition-metal thin films and magnetic nanostructures. The insights revealed through analyses of electronic structures are imperative in guiding experimental procedures for the design of novel magnetic nanomaterials. To precisely determine MAEs of magnetic nanostructures, a general treatment including noncollinear ordering is necessary, as was done in several recent studies [114,115]. The correlation and orbital polarization effects on MAEs should also be examined, especially when the system has only a few atoms [116,117]. Nevertheless, the torque method used in our studies captured the key factor of the origin of MAEs and hence offers effective approaches to engineer and control the MAEs.
It should be pointed out that the fabrication of the magnetic nanostructures proposed above requires delicate growth condition, so it is still difficult to achieve these magnetic nanostructures currently. However, intense experimental research on the magnetic anisotropy of magnetic nanostructures has been inspired and great progresses have been made since the discovery of a giant MAE for Co atoms adsorbed on a Pt(111) surface [53]. Interestingly, artificial TM nanostructures can be created on semiconducting surfaces using scanning tunneling microscope (STM) tip, in order to obtain intentional magnetic properties [105,[118][119][120][121]. The Fe atoms of a linear chain on Cu 2 N surface couple antiferromagnetically, and the bistable Néel states can be 'read' and 'write' by a spin polarized STM tip, therefore it may be used for information recording [105]. The MAE of a single Co atom has been boosted to the theoretical limit (~60 meV), when it is placed on MgO(100) surface [119]. However, the magnetic remanence can be realized only at extremely low temperature (<1 K) and the magnetic lifetime is short (~10 2 μs). This is because the magnetic states of a quantum magnet need to be protected from quantum tunneling of the magnetization as well as from scattering with electrons and phonons of the substrate, in addition to the magnetic anisotropy. In a recent experiment, the blocking temperature and magnetic lifetime were significantly enhanced for the Ho atom adsorbed on MgO(100) surface (30 K and 1500 s) [120,121]. Nevertheless, the blocking temperature is still too low to hold the magnetic remanence at practical temperatures. Our proposals provide new choices to realize magnetic bistability at practical temperatures with long magnetic lifetime. It can be expected that the atomic structures may be achieved with delicate growth conditions.