Giant magnetic anisotropy of heavy p-elements on high-symmetry substrates: a new paradigm for supported nanostructures

Nanostructures with giant magnetic anisotropy energies (MAEs) are desired in designing miniaturized magnetic storage and quantum computing devices. Previous works focused mainly on materials or elements with d electrons. Here, by taking Bi–X(X = In, Tl, Ge, Sn, Pb) adsorbed on nitrogenized divacancy of graphene and Bi atoms adsorbed on MgO(100) as examples, through ab initio and model calculations, we propose that special p-element dimers and single-adatoms on symmetry-matched substrates possess giant atomic MAEs of 72–200 meV, and has room temperature structural stability. The huge MAEs originate from the p-orbital degeneracy around the Fermi level in a symmetry-matched surface ligand field and the lifting of this degeneracy when spin–orbit interaction (SOI) is taken into account. Especially, we developed a simplified quantum mechanical model for the design principles of giant MAEs of supported magnetic adatoms and dimers. Thus, our discoveries and mechanisms provide a new paradigm to design giant atomic MAE of p electrons in supported nanostructures.


Introduction
The physics of few-atom structures on surfaces has initiated tremendous interest in recent years [1][2][3][4][5][6]. In particular, adatoms and dimers supported on substrates are the minimal functional units available, and in the meantime, very rich physics has been revealed, which predicts great potential for their applications in electronics, spintronics, catalysis and quantum computation [5,[7][8][9][10][11][12][13][14]. Nowadays, the design, fabrication, and manipulation of functional supported few-atom structures are attracting more and more attention in the field of nanoscience and nanotechnology.
An important topic of few-atom systems is their magnetic anisotropic energy (MAE), which is a key property of magnetic materials [15,16], and defined as the maximum energy difference of the system when the magnetic moment points to different directions. The MAE creates the energy barrier that stabilizes the magnetic moment. There are two reasons why the MAE of few-atom systems is important. First, the MAE of few-atom structures are hopefully several orders of amplitude larger than the MAE in bulk because of the reduced symmetry constrain [17], which is important for tremendous technological fields such as high-density-recording, quantum computation, and molecular spintronics [15,[18][19][20][21]. Second, the MAE depends on the occupation of specific orbitals, which provides the possibility to manipulate the magnetic properties with external fields [22][23][24]. The MAE can be linked to several experimental techniques such as x-ray absorption spectra, x-ray magnetic circular dichroism and inelastic electron tunneling spectroscopy [21,25,26]. In order to fit the requirement of application, systems with large MAE are desired. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
As the spin-orbital coupling constant λ increases with Z, where Z is the atomic number [21], dimers containing 5d-elements were predicated to give the largest MAE in the d electron systems [43]. The main group late pelements, which lie on the right hand side of d-elements in the periodic table, typically have larger λ than the transition metals in the same period [54]. Normally, these p electron systems are not prototype systems related to magnetism so that MAE has been studied for only very few organic radicals [54]. However, the isolated p atoms, such as Bi and Tl, are magnetic. By designing a suitable surface ligand filed, the p-element magnetization can be preserved on surfaces. More interestingly, there are only six p orbitals while the number of d orbitals is ten, hence the coupling between orbitals associated with spin-orbital interaction (SOI) is simpler to handle in the p electron systems, which makes the design of p-element adsorption structures much easier for the purpose of optimizing MAE. To our knowledge, no attempt has been made to construct supported magnetic structures with these pelements. Figure 1 shows our design principle. The degeneracy of p x/y orbitals in an isolated atom could be preserved in a ligand field with certain symmetry, for example, C v 4 in a square lattice. In the single electron occupation case, this degeneracy can be lifted under the perturbation of SOI when the spin of the electron is out-of-plane (detailed derivation will be given later). This so-produced splitting will generate an energy difference between the spins pointing to in-plane and out-of-plane directions. Thus, the desired MAE can be produced in such systems with heavy p-elements.
To substantiate the above idea, in this paper, we construct a series of room-temperature-stable atomic magnets with p-elements by placing a Bi adatom or Bi-X dimers (X = Ga, Ge, In, Sn, Tl and Pb) on symmetrymatched surfaces, i.e., p-element dimers on nitrogenized divacancies (NDV) of graphene (figures 2(a), (b)), and p-element adatoms on MgO(100) (figures 2(c), (d)). Note that the NDV and similar structures have been proposed and synthesized [55,56]. We find that most of these structures have extremely large MAE (up to 200 meV) with vertical easy-magnetization-axis, suggesting that they can be applied in magnetic storage and quantum spin processing. Based on these discoveries, we develop a quantum mechanical model and propose that in order to get giant MAE with heavy p-elements, one can try to create p x/y degenerate orbitals around the Fermi level and move the p z orbital away from the Fermi level as far as possible. This provides us a paradigm for p atomic magnetic storage-device design. To elaborate the principles to generate giant p electron MAE, first we focus on systems of Bi-X on NDV (Bi-X@NDV) of graphene, and then we demonstrate that this principle can be applied to other p electron magnetic systems.

Methods
Calculations were performed using the plane-wave-basis-set Vienna ab initio simulation package (VASP) [57]. A 4 3 8 supercell with a K-point sampling of 6×6 was used for dimers on NDV of graphene and a 3×3 supercell with a K-point sampling of 10×10 was used for Bi adatom on MgO(100). The vacuum layer was set to be 25 Å. Projector augmented wave potentials were employed with a kinetic energy cutoff of 500 eV [58]. For the exchange-correlation functional, the Perdew-Burke-Ernzerh of generalized gradient approximation was utilized [59]. Gaussian smearing with 0.05 eV was used to smooth the occupancy in all calculations and the PDOS result is compared with the tetrahedron method with Blöchl corrections [60]. The structural relaxation was end till the forces on each atom is less than 0.01 eV/Å. The relativistic Hamiltonian containing all relativistic corrections up to order α 2 (where α is the fine-structure constant) is recast in the form of 2×2 matrices in spinspace by re-expressing the eigenstates of the total angular momentum in terms of a tensor product of regular angular momentum eigenstates l m , ñ | and the eigenstates of the z-component of the Pauli spin matrices. The relativistic effective potential consists of a term diagonal in spin-space which contains the mass-velocity and Darwin corrections. The nondiagonal contributions in spin-space arise from the spin-orbit coupling but also from the exchange-correlation potential when the system under consideration displays a noncollinear magnetization density [61]. The MAEs were calculated using the magnetic force theorem [62] with the convergence criteria of 10 −7 eV and were further checked by selfconsistent computations in Bi-Pb and Bi-Tl in a 4×4 supercell. The magnetic force theorem suggests that for a magnetic system, the variation on the directions of magnetic moments will not cause too much variation on the charge density so that one can calculate the energy differences between different magnetic configurations with non-selfconsistent computations based on a frozen charge density. We use a scalar-relativistic calculation to get a charge density, and use this charge density to perform two non-selfconsistent relativistic calculations with the magnetic moments out-of-plane and inplane. The differences between the energies in these two calculations is used as MAE. Because in this way one still diagonalizes the full Hamiltonian so that it includes more effects than a normal perturbative treatment. Note that for a full selfconsistent calculation, the configuration of the magnetic moment pointing to a specific direction is not guaranteed so that the MAE cannot be obtained. Our selfconsistent calculation suggests the easy axis is in the out-of-plane direction which is in good agreement with the MAE calculated by magnetic force theorem and hence excludes the possibility of magnetic moment quenching in non-selfconsistent calculations.

Geometrical and electronic structures
For Bi-X@NDV, by examining the total energies of a set of possible structures (see supplementary material, available at stacks.iop.org/NJP/20/043056/mmedia), we find that for X=Ga, Ge, In , Sn and Tl, the vertical structure (shown in figures 2(a) and (b)) is stable at room temperature (see supplementary material). For X=Pb, the vertical structure is unstable (see supplementary material). Nevertheless, the properties of all the six vertical Bi-X structures will be discussed on the same footing in the following to gain more insight into the principle of p-electron magnetism for giant MAE. By examining the adsorption energies of X at NDV and Bi on X@NDV, we find it is possible to synthesize the vertical structures by first depositing atom X on the NDV@graphene and then adsorbing atom Bi on atom X (see supplementary material).
We present the spin moments (M S ) of Bi and X alongside with the MAE of Bi-X dimers in table 1. The MAE is defined as MAE=E x −E z , where E z or E x is the total energy of the system when the magnetization-axis is along axis z (surface normal) or axis x as shown in figure 2. The magnetic moment is mainly distributed on the Bi atom while atom X has negligible M S in all these systems. The M S of Bi with third group elements are 0.8 μ B for Ga, In and Tl, while those with the fourth group elements are 0.6 μ B for Ge, Sn and Pb, respectively. More importantly, we note that all these systems have giant values of MAE together with vertical easy axis except for Bi-Ga. The vertical easy axis is crucial for applications because the energy difference between magnetization along in-plane directions is typically small (no more than 10% of the corresponding MAE in our cases) so that only vertical easy axis can uniquely lock the spin. The MAEs of Bi-Ge, Bi-In and Bi-Sn are comparable to the largest MAE reported for d-element dimers [51]. To our knowledge, the MAE of the Bi-Tl dimer of 203 meV is the largest value ever reported for the supported dimers, and this MAE value is very close to the record-breaking MAE value of an Os adatom on MgO [29]. Note also from table 1 that the MAE value increases along with the increasing atomic number of element X for X in the same group. However, as X is not spin-polarized, the MAE is from Bi.
In order to understand the origins of these giant MAEs as well as their dependence on the atomic number of X, next we study the spin-polarized projected density of states (PDOS) of the Bi atom in Bi-X@NDV both without and with the consideration of SOI. As the s orbital has no contribution to the SOI, only the contributions of the p orbitals are presented. Figure 3 shows the spin-polarized PDOS of Bi without SOI. Just as expected from our proposed principle in figure 1, a very clear feature of the PDOS of all the cases without SOI is that the p x and p y orbitals are quasidegenerate, which can be seen from the very small or even undiscernible energy difference Δ of these two orbitals. The quasi-degeneracy is associated with the approximately square ligand field of the NDV substrate in which the nearest neighbor N-N distances are very close as 2.71 and 2.81 Å. The difference between the N-N bonds of these two directions is the reason of generating Δ. And the Coulomb repulsions Un n   between the two electrons with different spins on the fully occupied orbital may be the reason that Δ is larger in X of group VI than in X of group III (Because in the former case one minority spin-orbital is fully occupied on Bi, which causes a Coulomb repulsion and distinguishes the double occupied and single occupied orbitals , while in the later case, there is only partially occupancy on such orbital which leads to a weak Coulomb repulsion and a smaller Δ). In order to be simplified, we will approximate the symmetry of the system by C v 4 rather than C v 2 in the following text. These degenerate orbitals are fully spin-polarized, and the minority spin p x/y orbitals just lie around the Fermi level in all systems except for Bi-Ga. The p z orbitals splits into two levels, with one far below the Fermi level and the other above it.
After taking SOI into account, the spin-projected PDOS of each p-orbital is shown in figure 4 for all the six cases. Here the spin-projected PDOS is obtained by the projected total DOS onto the two eigenstates of S z where z  is the direction of the magnetization direction of the total spin. For each case, two magnetization directions along the z axis (vertical direction) and x axis (parallel direction) are considered. When the magnetization is in the parallel direction, the PDOS changes quite neglectable as the one without SOI. However, when the magnetization is in the vertical direction, the PDOS peak of the p x/y around the Fermi level mix and split into two peaks labeled with p1 and p2. Since these p x and p y orbitals just lie around the Fermi level, the energysplitting which causes the distinct PDOS change between vertical-and parallel-magnetization has a direct relation to MAE (see below). In contrast, the change of p z orbitals is negligible in all cases due to m l = 0.

A simplified model
Such behaviors of the p orbitals with SOI in figure 4 and the origin of the giant MAE can be qualitatively understood by a simplified quantum mechanical model, i.e., a combination of a two-level Hamiltonian with a single electron occupation (for quasi-degenerate orbitals) and a second-order nondegenerate perturbation theory (for nondegenerate orbitals) [63]. These two theories will be used to deal with the couplings of the quasidegenerate p x −p y orbitals and the nondegenerate p x/y −p z orbitals, respectively. For most systems studied here, a simple two-level Hamiltonian is enough to give a clear physical picture to understand the energy lowering due to SOI, while for other systems, a two-level model is not enough and the nondegenerate perturbation theory is required to treat some higher order contributions, which will be seen below.
The SOI Hamiltonian can be approximately written as V L S l =ˆ·ˆ, where λ is an orbital and element dependent parameter containing the radial part of the orbital wave function. The angular and spin parts of p wavefunctions can be expressed as [64] p Y Y 2 1 x 1 1 1 where±labels the majority and minority spin, χ is the wave function in the spin-space, Y l m is the spherical harmonics. Since in our systems, the p z orbital is far from p x and p y in energy, we can work in a subspace spanned by the degenerate p x and p y orbitals first. The Hamiltonian H 0 expanded in such a subspace with the minority spin (in figure 3) can be expressed as a diagonalized matrix with diagonal elements being 0 and Δ, where Δ is the energy difference between the quasi-degenerate minority spin p x and p y orbitals around the Fermi level calculated by DFT without SOI (as shown in figure 3). When the magnetic moment is in the vertical direction, the SOI Hamiltonian is L S where L  and S  are angular moment ladder operator and spin ladder operator. By calculating the matrix elements between p x -( )and p y (−), we find that all these values are zero. So including SOI will not have much influence in energy for parallel magnetization. Thus, the MAE caused by degenerate orbitals can be expressed as where MAE D is the contribution of quasi-degenerate orbitals to MAE. We should note that the total energy lowering in the vertical magnetization direction requires the degenerate orbitals located very close to the Fermi level. Therefore, even the SOI can also cause the mixing of p x/y and splitting of spin-up orbitals, in all the systems we investigated, the spin-up orbitals are far away from the Fermi level thus have no contribution to MAE. The different sign of MAE of Bi-Ga can also be understood by the same reason. For Bi-Ga (see figure 3(a)), although it has the degenerate p x/y orbitals, they are far away from the Fermi level. The lowest unoccupied molecular orbital (LUMO) and highest occupied molecular orbital (HOMO) are at 0.65 and −0.8 eV, respectively. Since the λ of Bi is approximately 0.8 eV [65], the splitting by including SOI is λ/2 = 0.4 eV which is insufficient to change the occupation. This is verified in figure 4(a). Bi-Ga and its MAE of −179 meV cannot be understood by the above two-level model in the degenerate p x -p y subspace. Instead, the coupling between the nondegenerate orbitals may give the dominating contribution to the MAE in this case. In the framework of a nondegenerate perturbation, the first-order contribution is zero [66,67]. Up to the second-order, the correction from SOI perturbation is where u and o denote unoccupied and occupied, respectively, E u and E o are energies of corresponding orbitals, and E SOI ND is the SOI energy contributed by the coupling between the nondegenerate orbitals. So up to the secondorder perturbation, E 0 SOI ND  . Because the SOI is inversely proportional to the energy separation between E u and E o , from figure 3(a) the most dominating coupling may arise from the occupied p z around −2 eV and the occupied p x/y orbitals at −0.8 eV, coupling to the unoccupied p x/y around 0.6 eV. When the magnetization is in the z direction, all the above mentioned matrix elements are zero; but when the magnetization is along the x axis, . Therefore, the system is more stable when being magnetized in-plane. We should mention that such effects also exist in other Bi-X cases, with which it is easy to interpret the variation of MAE over the element number of atom X as shown in table 1. As mentioned above, the nondegenerate part is contributed from the coupling between p z and p y . In such cases, from equation (8) one can see that E SOI ND is proportional to the occupation of p z . By integrating the PDOS of p z shown in figure 3, we find that for those X in the same group, the occupation of Bi p z orbital becomes smaller with the increasing of element number (see supplementary material). The reduction of p z may be attributed to the reducing of covalency in the Bi-X bond which results from the increasing of metallicity of X. Namely, the decrease of p z electrons leads to the negative energy decrease due to the SOI in the parallel cases, and thus the MAE increases with atomic number in the same group. Finally, we tired DFT+U with U = 2 eV on Bi in Bi-Tl but found negligible influence on the PDOS of Bi, which means the onsite correlation is not important in our systems.
The above two mechanisms are illustrated in figure 5. In the presence of degenerate p x/y orbitals, the SOI between the p x and p y cannot lift their degeneracy when M s is in-plane but lift the degeneracy when M s is out-ofplane (M ⊥ ) and hence lowers the energy of the system. For the nondegenerate orbitals, the SOI between p z and p x/y has negligible effect when M s is out-of-plane; but it has significant effect on the energy lowering of the system when M s is in-plane (M P ). The final results depend on the competition between these two effects from both degenerate and nondegenerate orbitals. As λ<E u −E o typically, the coupling between nondegenerate orbitals is a higher-order smaller contribution compared to the coupling between the degenerate orbitals. So usually the degenerate effect dominates. We noticed that in some of the previous investigations, the coupling between the degenerate orbitals was overlooked in the interpretation of MAE, in spite of the presence of the degenerate orbitals [43,51,53]. Our work points out the necessity to consider the effects of degenerate orbitals and makes the picture on the physical origin of giant MAE complete.
The above two mechanisms suggest a general picture to realize giant MAE. In order to generate vertical magnetization, degenerate orbitals around Fermi level are recommended. This can be obtained through manipulating an adatom or dimer on a symmetry-matched substrate by electric fields or by doping the substrate. In order to reduce the in-plane MAE (associated with nondegenerate orbital coupling), the orbital splitting of m l =0 needs to be enlarged. As atom X lies above the graphene surface, the Bi atom which is responsible for the giant MAE has a larger distance from the substrate than an adatom has in adatom/MgO [27]. This implies that the influence of the surface spin-flip scattering will be smaller in our cases so the system we proposed is hopeful to have a longer spin relaxation time. Note that in a certain surface ligand field, the splitting of degenerate orbitals under a SOI could also take place for adatom/dimer with d or f orbitals. And this effect in d systems have been reported early [68][69][70][71]. However, the fresh material in p electrons is that for p electron systems, p z is the only orbital which is out-of-plane, so it is possible to tune this orbital separately without bringing a strong influence on the other orbitals. For example, if an enhanced interaction in the normal direction is given by the p- Figure 5. Illustration of the competition of two mechanisms of MAE from degenerate (upper part) and the nondegenerate (lower part) orbitals. The red lines indicate initially degenerate p x/y energy levels without SOI, the orange arrow indicates spin, and the blue letters indicate the effective contributions to MAE (M ⊥ for the degenerate case and M P for the nondegenerate case). element X, the p z splitting will be enlarged, but the p x/y orbitals may not shift too much. But it is impossible to tune d z 2 in such a way without bringing influence on d xz/yz . This suggests the p systems may have more convenience in an engineering view.
The orbital of p z always contribute to in-plane MAE. Thus, these conclusions should also hold for a broad range of magnetic p-electron systems. To further demonstrate that the mechanisms shown in figure 5 are generally applicable, we investigate a single Bi adatom on MgO(100). The most stable geometric structure is the atop adsorption on oxygen (figures 2(c) and (d)), and accordingly the quasi-degenerate orbitals appeared around the Fermi level with an MAE being 75 meV, which again suggests the proposed mechanisms hold in other analogical systems (see supplementary material).

Conclusion
In summary, we discovered that most of the vertical structures of Bi-X@NDV of graphene can be stable at room temperature, and Bi is spin-polarized in these systems. Most of these dimers show giant MAE with easy axis in the out-of-plane direction. In particular, the MAE of Bi-Tl@NDV is comparable to the largest MAE ever reported. The origin of giant MAE can be attributed to the SOI correction to the degenerate orbitals near the Fermi level which are preserved by the surface ligand field. Such degeneracy is relevant to an approximated local C v 4 symmetry of the system. These findings provide a general method to design giant MAE with p electron magnetism: one could find a substrate with a proper symmetry (approximated C v 4 for p in this paper) to create quasi-degenerate orbitals and then tune these orbitals close to Fermi level; and simultaneously move the orbitals with p z far from Fermi level. These discoveries and insights pave a pathway for further development of robust nanomagnetic units.