Ferromagnetism in nitrogen doped graphene

TL;DR

This paper uses first-principles calculations to explore the mechanisms behind ferromagnetism in nitrogen-doped graphene, identifying key defect complexes and proposing a direct exchange mechanism, with co-doping enhancing magnetic properties.

Contribution

It provides a detailed microscopic understanding of ferromagnetism in nitrogen-doped graphene and suggests strategies to enhance magnetic properties through defect engineering and co-doping.

Findings

Nitrogen complexes contribute to local and itinerant magnetic moments.

Magnetic interactions depend on defect type, concentration, and geometry.

B co-doping enhances ferromagnetic behavior.

Abstract

Inducing a robust long-range magnetic order in diamagnetic graphene remains a challenge. While nitrogen-doped graphene is reported to be a promising candidate, the corresponding exchange mechanism endures unclear and is essential to tune further and manipulate magnetism. Within the first-principles calculations, we systematically investigate the local moment formation and the concurrent interaction between various defect complexes. The importance of adatom diffusion on the differential defect abundance is discussed. The individual nitrogen complexes that contribute toward itinerant and a local magnetic moment are identified. The magnetic interaction between the complexes is found to depend on the concentration, complex type, sublattice, distance, and orientation. We propose that the direct exchange mechanism between the delocalized magnetic moment originating from the itinerant $\pi$-electron at the prevalent graphitic complexes to be responsible for the observed ferromagnetism. We show that B co-doping further improves ferromagnetism. Present results will assist in the microscopic understanding of the current experimental results and motivate experiments to produce robust magnetism following the proposed synthesis strategy.

Tables

Table 1. Table 1: Magnetic moment of the defect complexes in Figure  1 severely depends on the treatment of exchange-correlation energy. We compare the moments calculated with conventional PBE, hybrid HSE06 ( a HF subscript 𝑎 HF a_{{}_{\rm HF}} =0.25) and SCAN meta-GGA functionals, where the later better describes the moment. While the HSE06 calculations were carried out with 4 × \times 3 ( 3 3 \sqrt{3} , 3) a 0 subscript 𝑎 0 a_{0} supercell, the 10 × \times 5 ( 3 3 \sqrt{3} , 3) a 0 subscript 𝑎 0 a_{0} supercell was used for the calculations with PBE and SCAN density functionals.

Defect	Magnetic moment (${\mu }_B$)
	PBE	HSE	SCAN
(a) Monovacancy	1.5	1.6	2.0
(b) N adatom	0.9	1.0	1.0
(c) Graphitic N	0.0	0.1	1.0
(d) Meta-graphitic N pair (NAA)	0.0	0.4	2.0
(e) Para-graphitic N pair (N${}_{{\mathrm{AB}}^{\prime }}$)	0.0	0.0	0.0
(f) Graphitic 3N (G3N)	0.9	1.0	1.0
(g) Triazinic N (T)	0.9	1.0	1.0
(h) Trimerized pyridinic N (TPY)	0.2	1.0	1.0
(i) Trimerized pyrolic N (TPI)	1.0	1.0	1.0

Full text

I Introduction

The absence of intrinsic magnetism in graphene limits its application in spintronics and is being actively persuaded since it was first isolated. (Novoselov et al., 2004; Tucek et al., 2018) However, the magnetism can be induced in graphene due to the $p_{z}$ electron imbalance in bipartite hexagonal sublattice, (Boukhvalov et al., 2008; McCreary et al., 2012; González-Herrero et al., 2016) electron localization at the zigzag edge, (Nakada et al., 1996; Fernández-Rossier and Palacios, 2007; Yazyev and Katsnelson, 2008; Kabir and Saha-Dasgupta, 2014; Magda et al., 2014; Ganguly et al., 2017; Slota et al., 2018) interaction with a spin-orbit coupled ferromagnetic substrate, (Swartz et al., 2012; Wang et al., 2015; Leutenantsmeyer et al., 2017) structural defect creation, (Esquinazi et al., 2003; Yazyev and Helm, 2007; Červenka et al., 2009; Ugeda et al., 2010; Nair et al., 2012, 2013; Zhang et al., 2016) and heteroatom doping. (Krasheninnikov et al., 2009; Blonski et al., 2017) However, a long-range magnetic order in graphene can only emerge through the exchange coupling between these induced local moments.

Hydrogen absorption on graphene creates an imbalance in $p_{z}$ electrons between the two sublattices, which in turn generates a net magnetic moment. (González-Herrero et al., 2016) Further, the spin-polarized state is observed to extend over several nanometers through the scanning tunneling microscopy (STM) experiment. This spatially extended magnetic state provides a channel for direct coupling between the induced magnetic moments at large distances. The half-filled flat band at the Fermi level produces localized magnetic moments at the zigzag edge of a graphene flake, (Nakada et al., 1996) and a long-range order has been predicted. (Fernández-Rossier and Palacios, 2007; Yazyev and Katsnelson, 2008; Kabir and Saha-Dasgupta, 2014; Magda et al., 2014; Ganguly et al., 2017; Slota et al., 2018) However, the edge magnetism is predicted to be very intricate as it is strongly affected by the shape and size of the flakes, charge doping, chemical functionalization at the edge, and on-site Coulomb interaction. (Kabir and Saha-Dasgupta, 2014; Ganguly et al., 2017) Proximity-induced exchange-coupled ferromagnetic order sets in graphene sheet at room temperature while places on an atomically flat ferromagnetic thin film insulator. (Swartz et al., 2012; Wang et al., 2015; Leutenantsmeyer et al., 2017) The smallest structural defect, a lattice vacancy in graphene generates a semi-localized magnetic moment due to the $\pi$ electron imbalance in addition to a $\sigma$ contribution from the dangling bond. (Yazyev and Helm, 2007; Ugeda et al., 2010; Nair et al., 2012, 2013; Zhang et al., 2016) Thus, proton irradiated graphene and graphene with a network of point defects show magnetic ordering. (Esquinazi et al., 2003; Červenka et al., 2009) Similarly, the larger lattice defects such as grain boundaries with dangling bonds are theoretically predicted to have magnetic moment. (Akhukov et al., 2012) However, the experimental evidence of a magnetic ordering remains scarce and controversial. (Červenka et al., 2009; Nair et al., 2012; Wang et al., 2009; Matte et al., 2009; Sepioni et al., 2010; Ney et al., 2011) While a room temperature ferromagnetism is predicted, (Červenka et al., 2009; Wang et al., 2009; Matte et al., 2009) defected graphene with vacancies and other $sp^{3}$ defects lead to paramagnetism and devoid of any magnetic ordering down to low temperature. (Nair et al., 2012; Sepioni et al., 2010; Ney et al., 2011) Furthermore, in the cases where ferromagnetic signals have been detected, the origin of such long-range ordering remained unexplained. The presence of $d$-electrons in transition-metal impurities induces a local moment in pristine and defected graphene. (Krasheninnikov et al., 2009; Nafday et al., 2016) However, only a paramagnetic ordering has been observed in this transition-metal doped graphene samples that were investigated through the X-ray magnetic circular dichroism measurement. (Eelbo et al., 2013)

In contrast, the substitutional and $sp^{3}$ impurities have attracted much attention with intriguing but elusive magnetic properties. The itinerant $sp$-electron Stoner ferromagnetism in a narrow impurity band is fundamentally different from the magnetism induced by the $d$-band of transition metals. (Edwards and Katsnelson, 2006) Such $sp$-electron ferromagnetism may lead to a higher Curie temperature than the usual dilute magnetic semiconductors. The adatom impurities like hydrogen and fluorine exhibit strong absorption and induce magnetic moments (Boukhvalov et al., 2008; McCreary et al., 2012; Hong et al., 2011; Kim and Cho, 2013) due to sublattice imbalance in $\pi$-electrons following Lieb’s theorem. (Lieb, 1989) While the hydrogenated graphene show magnetic order at 5 K, (González-Herrero et al., 2016) the fluorine doped graphene reveals paramagnetic behaviour. (Červenka et al., 2009; Nair et al., 2013) In this context, the substitutional N-doped graphene attracted much attention, which exhibits robust ferromagnetism. (Ito et al., 2015; Blonski et al., 2017; Liu et al., 2016; Li et al., 2015; Miao et al., 2016; Khurana et al., 2013) However, the corresponding Curie temperature and saturation magnetization are reported to show a strong dependence on the growth, annealing strategies, and nitrogen concentration, which hinder a careful microscopic analysis. A variety of nitrogen complexes are generated during the non-equilibrium synthesis, while N-concentration and the type of defect complex determines the electronic and transport properties. (Blonski et al., 2017; Ito et al., 2015; Blonski et al., 2017; Liu et al., 2016; Li et al., 2015; Miao et al., 2016; Khurana et al., 2013; Zhao et al., 2011; Usachov et al., 2011; Wang et al., 2012; Lv et al., 2012; Zabet-Khosousi et al., 2014; Tison et al., 2015a; Fujimoto and Saito, 2011; Lherbier et al., 2013; Deng et al., 2011) Recently, the room-temperature ferromagnetism was reported in boron and nitrogen-doped turbostratic carbon films grown via pulsed laser deposition technique, while the saturation magnetization is observed to increase with increasing N-doping. (Mondal et al., 2018) However, even with such efforts, the nature and microscopic mechanism of magnetic interaction are not well understood in N-doped graphene, which is necessary to tune and manipulate magnetism toward its successful applications.

Within the first-principles calculations, here we systematically investigate the electronic and magnetic properties of various N-defect complexes in graphene. Further, how these defects interact with each other and lattice vacancies are studied to understand the mechanism responsible for the observed ferromagnetism. We report that the long-range magnetic interactions show a strong dependence on the N-concentration, complex type and on which sublattice they reside, distance, and orientation. We observe that the magnetism in the graphitic and vacancy containing N-complexes are fundamentally very different, and propose that the observed ferromagnetic order originates from the direct exchange between the extended magnetic states of the graphitic defects. Further, we investigate the crucial role of boron co-doping in the ferromagnetism.

II Computational details

The spin-polarized density functional theory based calculations (Kresse and Hafner, 1993; Kresse and Furthmüller, 1996) were performed within the projector augmented wave formalism, (Blöchl, 1994), and the wave functions were expanded in a plane-wave basis with 600 eV kinetic energy cut-off. The exchange-correlation energy was described using the conventional Perdew-Burke-Ernzerhof (PBE) form for the generalized gradient approximation (GGA). (Perdew et al., 1996) The structures were fully relaxed until all the forces became smaller than 0.015 eV/Å  threshold using the PBE functional. The Brillouin zone was sampled using a $\Gamma$-centered 8$\times$6$\times$1 $k$-grid for the structural optimization, and a much dense $\Gamma$-centered 24$\times$18$\times$1 $k$-grid was used for the evaluation of density of states (DOS) and magnetization density. The subsequent N-doping was carried with the 4$\times$3$(\sqrt{3},3)a_{0}$ and 10$\times$5$(\sqrt{3},3)a_{0}$ rectangular supercells, where $a_{0}$ is the C–C bond length in graphene. The periodic images of the supercell along the direction perpendicular to the graphene sheet are separated by at least 12 Å distance. It is known that the self-interaction error inherent to the conventional GGA functional results in spurious electron delocalization and incorrect magnetic solutions for adatom doping and vacancy defects in graphene. (Kim and Cho, 2013; Ronchi et al., 2017; Valencia and Caldas, 2017) In this regard, we recalculated the electronic and magnetic properties using the Heyd-Scuseria-Ernzerhof (HSE06) (Heyd et al., 2003) hybrid functional as well as using the strongly constrained and appropriately normed (SCAN) meta-GGA functional. (Sun et al., 2015) While the inclusion of fractional exact exchange for HSE06 hybrid functional improves the description of the electronic structure and the magnetization density, it is computationally expensive. In contrast, since the semilocal SCAN meta-GGA density functional satisfies all the known exact constraints for the exchange-correlation energy and the self-interaction error is reduced, it should better describe the thermodynamic, electronic and magnetic properties of materials. Indeed the SCAN functional predicts these properties correctly for diverse materials with a similar computational cost to the conventional GGA. (Sun et al., 2016; Isaacs and Wolverton, 2018) In the present context of defect magnetism and the concurrent interaction between them, we further validate the applicability of the SCAN functional. The results are obtained within the SCAN functional if it is not otherwise stated.

III Results and Discussion

Before we investigate the magnetic interactions between the moment generating defect complexes in N-doped graphene, we first address the local magnetism produced by these individual defects. In addition to the lattice vacancies, the common defects in N-doped graphene are chemisorbed, graphitic and vacancy containing nitrogen complexes (Figure 1). (Blonski et al., 2017; Ito et al., 2015; Blonski et al., 2017; Liu et al., 2016; Li et al., 2015; Miao et al., 2016; Khurana et al., 2013; Zhao et al., 2011; Fujimoto and Saito, 2011; Usachov et al., 2011; Wang et al., 2012; Lv et al., 2012; Lherbier et al., 2013; Zabet-Khosousi et al., 2014; Tison et al., 2015a) As the treatment of exchange-correlation may affect the electronic and magnetic description of these defects, here we compare the PBE, hybrid HSE06, and SCAN meta-GGA functionals.

III.1 Magnetism of individual defects

Single vacancy. Magnetic properties of single vacancy defect in graphene are still studied and debated. Yazyev and Helm (2007); Nair et al. (2012, 2013); Ronchi et al. (2017); Valencia and Caldas (2017) The Jahn-Teller distorted V${}_{1}(5|9)$ vacancy [Figure 1(a)] generates defect induced localized $\sigma$ and $\pi$ states, Vσ and Vπ (Figure 2). Thus, the magnetic moment of V${}_{1}(5|9)$ vacancy has two components. (Zhang et al., 2016; Nair et al., 2013) The $\pi$-electron imbalance in the bipartite graphene lattice should generate 1 $\mu_{B}$ moment according to the Lieb’s theorem. Additionally, the singly occupied dangling $\sigma$-bond should provide with another 1 $\mu_{B}$ moment, leaving the total vacancy moment to 2 $\mu_{B}$. However, the $\pi$-magnetism is suppressed within the conventional GGA description of exchange-correlation due to improper delocalization of electrons in the Vπ orbital. Thus, the magnetic moment of V${}_{1}(5|9)$ is underestimated to 1.5 $\mu_{B}$ within the PBE exchange-correlation functional (Figure 2, Table 1), similar to that predicted earlier. Yazyev and Helm (2007) Further, the moment depends on the supercell size, and thus on the vacancy concentration. (sup, ) The picture is slightly improved while a fraction $a_{{}_{\rm HF}}$ of Hartree-Fock exact exchange is considered in the HSE06 functional and 1.6 $\mu_{B}$ moment is predicted at the vacancy for $a_{{}_{\rm HF}}$=0.25. A further increase in $a_{{}_{\rm HF}}$ to 0.4 ensures a 2 $\mu_{B}$ moment at the vacancy. In contrast, treating the exchange-correlation with SCAN meta-GGA functional describe the magnetism correctly with 2 $\mu_{B}$ moment for the vacancy (Figure 2, Table 1).

Similarly, the corresponding electronic structure depends on the exchange-correlation functional. While a metallic solution emerges within the PBE and HSE06 functionals, the SCAN functional correctly predicts a gapped band structure. This result is consistent with the diverging resistivity observed in irradiation induced defected graphene with about 1012 cm*-2* vacancy concentration, which indicate insulating behaviour. (Chen et al., 2009) This was latter theoretically corroborated in graphene with high vacancy concentrations and an insulating solution with a large bandgap was proposed. (Lherbier et al., 2012; Ulybyshev and Katsnelson, 2015) Further, while the defect induced Vσ state shows a large spin-splitting above 2 eV originating from the JT reconstruction, the splitting for the Vπ is small, about 250 meV (Figure 2). This prediction is consistent with the recent STM experiment, where the spin-polarized $\pi$ density of states (DOS) peaks at the defect site are observed to be separated by 20-60 meV. (Zhang et al., 2016) It should be noted here that in such experiments the Vπ spin-splitting may be affected by a metastable vacancy geometry or by its interaction with the substrate.

We find the vacancy magnetism to be robust against carrier doping with a realistic concentration, (Das et al., 2008) as the position of the Fermi level is tuned (Figure 2). The graphene with a neutral V${}_{1}(5|9)$ defect is intrinsically hole doped, and one extra electron on a large 10$\times$5 ($\sqrt{3}$, 3)$a_{0}$ supercell corresponds to 1.89$\times$1013 cm*-2* carrier density. The extra electron occupies the unfilled V${}_{\pi}^{\downarrow}$ state destroying the $\pi$-electron imbalance in the two bipartite sublattices. Thus, the $\pi$-component of the vacancy moment is completely quenched (Figure 2). In contrast, the $\sigma$-component of the magnetic moment is protected against carrier doping. Increasing the electron doping beyond 1.89$\times$1013 cm*-2* does not further affect the magnetic moment (Figure 2) At higher electron density the unoccupied bulk $\pi$-states get filled, while the unoccupied V${}_{\sigma}^{\downarrow}$ state lies much higher in energy and is impossible to populate with any realistic carrier density. The carrier concentration dependent differential charge density corroborates this picture. (Babar and Kabir, 2018) The scenario with hole doping is different only qualitatively. The $\sigma$-component is unaffected, and hole doping does not quench the $\pi$-moment completely indicating the particle-hole asymmetry in defected graphene (Figure 2). The semi-localized V${}_{\pi}^{\uparrow}$ state is only partially depopulated and a fractional $\pi$-moment survives. (Nair et al., 2013; Yndurain, 2014) Therefore, while the $\pi$ magnetism in graphene vacancy could be completely quenched with sufficient carrier doping, the spin-$\frac{1}{2}$ $\sigma$-magnetism is robust.

In the context of differential defect abundance in N-doped graphene, it is relevant to investigate the C and N adatom diffusion on graphene lattice. The C and N adatoms are absorbed at the bridge site with 1.51 and 0.9 eV binding energy. We calculate the activation barrier for adatom diffusion within climbing image nudged elastic band method. (Henkelman et al., 2000) In agreement with the previous calculations, (Lehtinen et al., 2003; Hou and Terakura, 2015) the calculated activation barriers are found to be 0.43 and 0.86 eV, respectively, for C and N diffusion (Supplemental Material). (sup, ) Such low activation barrier results into a very high diffusivity $D_{\rm C/N}$ at the experimental growth and annealing temperatures (Figure 3), and the adatoms will quickly recombine with the graphene vacancies. Therefore, due to high adatom mobility, it is expected that the graphene vacancies will be annihilated to decrease in concentration and from graphitic N-complexes, while annealed at high-temperature or non-equilibrium growth conditions. (Blonski et al., 2017; Ito et al., 2015; Blonski et al., 2017; Liu et al., 2016; Li et al., 2015; Miao et al., 2016; Khurana et al., 2013; Zhao et al., 2011; Usachov et al., 2011; Wang et al., 2012; Lv et al., 2012; Zabet-Khosousi et al., 2014; Tison et al., 2015a) This prediction is consistent with the experimental observation of pyridinic to graphitic defect conversion during postannealing at high-temperature. (Usachov et al., 2011) Further, the diffusion mediated self-healing has been observed in scanning transmission electron microscopy. (Zan et al., 2012; Robertson et al., 2013)

N-adatom. The chemisorbed N-adatom [Figure 1(b)] is a p-type defect with the impurity state above the Fermi level. A localized magnetic moment of 1 $\mu_{B}$ is generated from the in-plane N-orbitals (Table 1). However, due to the high N diffusion on graphene (Figure 3), such adatom defect may not be present and effectively contribute to magnetism in graphene samples that are prepared and annealed at high-temperature. In contrast, high N diffusion will further promote the formation of graphitic N-complexes.

Graphitic nitrogen defects. The appearance of Dirac point (DP, $E_{\rm D}$) in graphene band structure is very sensitive to external perturbations, and the breaking of sublattice symmetry induces a substantial gap at the DP. (Ohta et al., 2006; Giovannetti et al., 2007; Zhou et al., 2008) Thus, all the graphitic N-defects are expected to generate a gap at $E_{\rm D}$. Single graphitic nitrogen [Figure 1(c)] introduces an impurity state near the Fermi level and induces a 0.60 eV gap at $E_{\rm D}$. An overall semiconducting electronic structure emerges and the results are consistent with the STM measurement. (Tison et al., 2015b) Further, an exchange splitting of 18 meV is observed between the majority and minority spin channels, which results in itinerant $\pi$-magnetism along with a small 0.06 $\mu_{B}$ moment localized at the N-site. However, the total moment generated by this defect strongly depends on the treatment of exchange-correlation energy (Table 1) and defect concentration shown in Supplemental Material. (sup, ) While the SCAN functional correctly generates 1 $\mu_{B}$ moment, the moment strongly depends on the fractional exact exchange $a_{{}_{\rm HF}}$ in HSE06 functional. While the conventional HSE06 functional ($a_{{}_{\rm HF}}$=0.25) underestimates the moment (Table 1), it gradually increases with increasing $a_{{}_{\rm HF}}$ and generates 1 $\mu_{B}$ moment for $a_{{}_{\rm HF}}$=1. In contrast, in agreement with an earlier report, (Park et al., 2013) the conventional PBE functional fails to produce any spin-splitting and a non-magnetic solution emerges. Therefore, the SCAN meta-GGA functional provides the correct magnetic picture over other exchange-correlational functionals.

The graphitic complexes with two nitrogens NAA and N${}_{\rm AB^{\prime}}$ have been experimentally identified [Figure 1(d) and (e)]. (Deng et al., 2011; Lv et al., 2012; Zabet-Khosousi et al., 2014; Tison et al., 2015a) Due to the sublattice asymmetry in these n-type complexes, the DP is destroyed, and within the SCAN functional a 0.31 and 0.50 eV gap emerges at the $E_{\rm D}$ for NAA and N${}_{\rm AB^{\prime}}$ complexes, respectively. The overall electronic structure is found to be defect specific as a half-metallic (semiconducting) solution is observed for the NAA (N${}_{\rm AB^{\prime}}$) defect in the large 10$\times$5 ($\sqrt{3}$, 3)$a_{0}$ supercell. The overall results here are in agreement with the experimental results as the character of electronic DOS near the Fermi level is consistent with the STS measurements. (Lv et al., 2012; Tison et al., 2015a) The exchange-split bands generate itinerant magnetism and a small moment of 0.06–0.08 $\mu_{B}$ is generated at the N-sites, which are coupled ferromagnetically (antiferromagnetically) while the nitrogen atoms occupy the same (different) sublattice. This fact results into 2 and 0 $\mu_{B}$ moments for the NAA and N${}_{\rm AB^{\prime}}$ pairs, respectively (Table 1). While the results from the different exchange-correlation functionals are compared, a similar picture emerges as discussed for the single graphitic N-defect.

While the triazine (T) complex [Figure 1(g)] is thermodynamically slightly more favorable over the graphitic-3N (G3N) complex [Figure 1(g)], both of these defects are formed at the non-equilibrium growth conditions.(sup, ) The graphitic-3N complex could also be produced when a highly mobile carbon adatom diffuses into the lattice vacancy in the trimerized pyridinic and pyrrolic complexes in the bulk regions of graphene. In both the defects [Figure 1(f) and Figure 1(g)], all three nitrogen atoms occupy the same sublattice leading to a sublattice asymmetry, which in turn destroys the DP and consequently open up a gap at the $E_{\rm D}$ (Figure 4). The overall electronic structure is half-metallic for the T-complex, and the same is found to be semiconducting for the G3N-complex (Figure 4). It is important to note that such half-metallicity induced asymmetric electron transport for different spin components could have significant implication in spin-based devices. To gain further insight into the magnetism, we investigate the differential charge density and the corresponding magnetization density (Figure 5). The C$\rightarrow$N $\pi$-electron transfer creates Lieb’s imbalance and an exchange split $\pi$-band emerges, which is responsible for the itinerant magnetism that is confirmed form the magnetization density [Figure 5(a) and (b)]. In addition, a hybridized nitrogen impurity state is observed near the Fermi level, which is responsible for a tiny N-moment of $\sim$ 0.05 $\mu_{B}$. Note the C-atom at the centre of G3N-defect acquire a comparatively larger localized moment of 0.25 $\mu_{B}$ [Figure 5(b)]. These together result into a total moment of 1 $\mu_{B}$ for these graphitic defect complexes (Table 1).

Thus, in general, the electronic and magnetic structure of the different graphitic defect complexes are remarkably similar. First, all of them open up a gap at the $E_{\rm D}$ due to sublattice asymmetry, which is consistent with the the angle-resolved photoemission spectroscopy results.  (Usachov et al., 2011) Second, the exchange-split $\pi$-band generates itinerant magnetism, while the localized moment at the N-site is small about 0.03–0.09 $\mu_{B}$. While the overall electronic structure of the graphitic complexes are defect specific, all of them are found to be n-type, which are in overall agreement with the STS results. (Lv et al., 2012; Tison et al., 2015a, b)

Vacancy containing nitrogen defects. Now we explore the defect complexes wherein the N atoms are incorporated in a vacancy. The trimerized pyridine (TPY), where three N atoms occupy the same graphene sublattice surrounding a vacancy [Figure 1(h)] is energetically favoured over the mono and dimerized pyridines. (Fujimoto and Saito, 2011) The 1.33 Å C$-$N bonds are consistent with experimental values for pyridinic complexes. The N atoms form an equilateral triangle with 2.61 Å sides, which is larger than the triangular unrelaxed vacancy with 2.46 Å. Further, this complex is found to be p-type and the Bader analysis indicates a 0.9-1.12 $e$/N charge transfer to nitrogen from carbon.

The magnetism in such vacancy containing nitrogen defects is rather intricate, and a detailed investigation into the electronic DOS (Figure 4) and the magnetization density [Figure 5(c)] unravel the origin of 1 $\mu_{B}$ moment for this complex. Unlike the bare vacancy, the $\sigma$ component to the moment is completely ruled out due to the presence of occupied $\sigma$-state that originates from the undercoordinated N atoms (Figure 4). In contrast, the moment is of entirely $\pi$ character and originates from the localized N-$p_{z}$ and itinerant C-$p_{z}$ states with an exchange-splitting of 11 meV. As a result, a sizeable localized moment of 0.15 $\mu_{B}$ is observed at the N-sites for this complex and a fraction of the total moment originates form the itinerant magnetism [Figure 5(c)]. In contrast to the SCAN functional, the PBE results provide a very contrasting picture. First, the moment is substantially underestimated (Table 1), and further, the moment originates from the localized $\sigma$-state at the N-sites without any $\pi$ contributions. However, the STM and STS investigations of the monomerized and trimerized pyridine reveal an electron density distribution that is consistent with the $p_{z}$ states of the vacancy defect and thus corroborate the results obtained with the SCAN functional. (Tison et al., 2015b; Kondo et al., 2012; Zhang et al., 2016)

The trimerized pyrrolic (TPI) complex is formed when a single N atom is absorbed in a divacancy to form a pentagon and rest of the undercoordinated C atoms are also substituted by N in Figure 1(i). The pentagonal C$-$N bonds (1.41 Å) are larger in comparison with the other C$-$N bonds (1.33 Å). Similar to the TPY-complex and bare vacancy, this complex has p-type character. The electronic DOS calculated within the PBE and SCAN functionals show qualitatively similar features, and the unoccupied $\sigma$-states originating from the N-impurities appear just above the Fermi level (Figure 4). The calculated magnetization density indicates that a large fraction $\sim$ 80% of the total 1 $\mu_{B}$ moment is localized on the nitrogen atoms and only a small fraction originates from the delocalized $\pi$-band [Figure 5(d)]. Further, due to the proximity of the unoccupied $\sigma$-states to the Fermi level, a slight perturbation through the charge doping will adversely affect this localized $\sigma$ component of the moment.

Therefore, the magnetism in the vacancy containing N-complexes is fundamentally different from the graphitic defects. While the magnetism in graphitic defects is entirely itinerant, a very significant fraction of the moment is localized at the N-sites in vacancy containing N-complexes. Thus, in the absence of a strong exchange path, one would not expect that such localized magnetic moments in these vacancy containing defect complexes to contribute significantly to the high-temperature ferromagnetism. Indeed, the role of this particular defect on the observed ferromagnetism has remained controversial.  (Ito et al., 2015; Blonski et al., 2017; Li et al., 2015)

III.2 Magnetic interactions

Although nitrogen-doped graphene is reported to show ferromagnetism, the microscopic mechanism is not yet understood and resulted in conflicting arguments. (Ito et al., 2015; Blonski et al., 2017; Liu et al., 2016; Li et al., 2015; Miao et al., 2016; Khurana et al., 2013) The appearance of a diverse Curie temperature with varied growth and annealing conditions makes the microscopic analysis very complicated as the differential defect abundance may fluctuate drastically in the samples. A long-range magnetic order in nitrogen doped graphene can only emerge when the moment generating defects are exchange coupled. Thus, we investigate the magnetic interactions between the abundant defect complexes.

The single vacancy in graphene lattice generates a robust semi-localized moment. However, the interaction between the vacancy moments is sublattice dependent. Vacancies in the same sublattice are coupled ferromagnetically, whereas they prefer antiferromagnetic alignment when placed on different sublattices. Further, the magnetic interaction sharply falls to zero below 1 nm distance due to semi-localized nature of the moment. (Scopel et al., 2016) Thus, a randomly distributed vacancy in graphene will not lead to any magnetic ordering, which is consistent with the results predicting paramagnetism. (Nair et al., 2012; Sepioni et al., 2010; Ney et al., 2011) Moreover, it should be noticed in this context that, the existence of a vacancy in the N-doped graphene is rare due to the very high C and N adatom mobility at the high-temperature growth conditions (Figure 3), which leads to self-healing of vacancy or the formation of graphitic N-defect.

Single graphitic nitrogen defect is similar to the pristine vacancy. For two substitutional atoms on the (different) same sublattice, the spin-configuration is (anti-) ferromagnetic. Moreover, the interaction strength decays fast with increasing separation and follows an inverse power law, and the long-range interaction may develop through Ruderman-Kittel-Kasuya-Yosida (RKKY) mechanism. (Park et al., 2013) Thus, a statistically random distribution of single graphitic nitrogen will not result in high-temperature ferromagnetism that has been reported so far. (Ito et al., 2015; Blonski et al., 2017; Liu et al., 2016; Li et al., 2015; Miao et al., 2016; Khurana et al., 2013)

A very complex picture emerges when we explore the interactions between the larger defects. Here, we address two distinct regimes of N-doping that is motivated by several experimental observations. (Blonski et al., 2017; Liu et al., 2016; Miao et al., 2016; Mondal et al., 2018) It has been unanimously pointed out in the literature that ferromagnetism emerges at a critical N-concentration around 5 at. % and above which the saturation magnetization and coercive field increase further with increasing the concentration. Thus, we consider N-concentration below 5 at. % and in the 8$-$13 at. % range, and investigate how the concentration influences the magnetic interaction. In this regard, we calculate the interaction energy $E^{\mathcal{DD}^{\prime}}_{\rm int}$ between two magnetic defects $\mathcal{D}$ and $\mathcal{D}^{\prime}$ as, $E^{\mathcal{DD}^{\prime}}_{\rm int}=(E^{\mathcal{DD}^{\prime}}_{\rm m}-E^{\mathcal{DD}^{\prime}}_{\rm nm})-(E^{\mathcal{D}}_{\rm m}-E^{\mathcal{D}}_{\rm nm})-(E^{\mathcal{D}^{\prime}}_{\rm m}-E^{\mathcal{D}^{\prime}}_{\rm nm})$ and shown in the Figure 6. The energies corresponding to the magnetic and non-magnetic solutions for the interacting defects are $E^{\mathcal{DD}^{\prime}}_{\rm m}$ and $E^{\mathcal{DD}^{\prime}}_{\rm nm}$, while $E^{\mathcal{D/D}^{\prime}}_{\rm m}$ and $E^{\mathcal{D/D}^{\prime}}_{\rm nm}$ represent magnetic and non-magnetic energies of the isolated defects.

First, we discuss the long-range magnetic interaction between the graphitic defects in the low N-concentration limit and investigate how it is affected by the spatial separation and direction between the defects (Figure 6). Two interacting T$-$T, G3N$-$G3N, and mixed T$-$G3N defects in a (10$\times$5)-($\sqrt{3}$,3)$a_{0}$ supercell represents 3 at. % N-concentration. At such a low concentration and for all these types of defect interactions, we notice some general trends [Figure 6(a) and 6(b)]. (i) These defects mostly prefer to occupy the same graphene sublattice, (ii) owing to the itinerant nature of the individual moments in these graphitic defects, the interaction between them is long-range, and (iii) the nature of magnetic interaction between the defects is entirely dictated by their respective sublattice [Figure 7(a)–(d)]. For the defects on the equivalent (opposite) graphene sublattice, the magnetic interaction is found to be ferromagnetic FM (antiferromagnetic AFM). The picture remains unaltered regardless of the spatial separation and relative orientation between the defects along the armchair and zigzag directions (Figure 6). Although the nature of the magnetic moment in these defects is completely different, the sublattice dependent magnetic interaction is in complete agreement with results predicted by the RKKY interaction between the localized moments on the bipartite graphene lattice. (Saremi, 2007; Brey et al., 2007) While these defects mostly prefer to occupy the same sublattice $\mathcal{D}_{\rm A}\mathcal{D}^{\prime}_{\rm A}$ over the $\mathcal{D}_{\rm A}\mathcal{D}^{\prime}_{\rm B}$ configuration, during the non-equilibrium growth process the defects are equivalently distributed over the two graphene sublattice. Therefore, the sublattice dependent FM and AFM interaction will restrict the observation of FM order at such a low concentration of nitrogen. This picture is consistent with the experimental observation that FM ordering is not found below 5 at. % N-concentration. (Blonski et al., 2017)

Now we turn our attention onto high N-concentration, and two interacting T$-$T, G3N$-$G3N, and mixed T$-$G3N defects in a (4$\times$3)-($\sqrt{3}$,3)$a_{0}$ supercell correspond to 12.5 at. % N-doping. The sublattice dependent magnetism for such high N-concentration changes dramatically. In contrast to the AFM interaction, a FM interaction with spontaneous magnetization emerges while the defects that are placed on the different sublattice [Figure 7(e)–(g)]. However, such FM structure is found when the defects are interacting along a particular direction, which is along the zigzag direction for T$-$T interaction [Figure 7(e)] and along the armchair direction for the G3N$-$G3N [Figure 7(f)] , and mixed T$-$G3N [Figure 7(g)] defect interactions. Thus, at high N-concentration, the graphitic T and G3N defects are responsible for the observed ferromagnetism. Further, the concentration of the graphitic T-complex may be increased by using the triazine-based precursor during the growth using the chemical vapour deposition technique. (Lu et al., 2013) Moreover, a suitable postannealing strategy may convert the vacancy containing pyridinic and pyrrolic defects to graphitic G3N-complex. (Usachov et al., 2011) Therefore, we conclude that such sample preparation strategies may give rise to robust ferromagnetism with improved Curie temperature. This discussion is consistent with the recent observation of room-temperature ferromagnetism in turbostratic N-doped graphene films, where the T-defect is found to be the most abundant complex. (Mondal et al., 2018)

Further, in the high concentration limit of 8.3 at. % N-doping, we investigate the interactions between the two NXX defects, and between the T-complex and a single graphitic N-defect. Both these interactions follow the usual sublattice dependent magnetic couplings. While the NAA$-$ NAA and TAAA$-$NA interactions are found to be ferromagnetic, the NAA$-$ NBB and TAAA$-$NB are observed to be antiferromagnetic. Therefore, such defect-defect interactions are not responsible for the observed ferromagnetism in the N-doped graphene.

In contrast to the graphitic N-defects, the vacancy containing complexes do not contribute to ferromagnetism. Moreover, in accordance with the experimental observation, (Usachov et al., 2011) we have discussed earlier that the C-adatom is highly mobile on the graphene lattice (Figure 3), and thus the vacancy containing TPY-defect complex will be converted into the graphitic G3N-complex during high-temperature growth and annealing. Therefore, the concentration of this defect will be negligibly small. However, for completeness, we have investigated the magnetic interaction between two TPY-complexes in both the low and high concentration limits. Unlike the graphitic T- and G3N-complexes and independent of N-concentration, the TPY$-$TPY magnetic coupling is AFM (FM) for defects on the opposite (equivalent) sublattice. Further, the picture is independent of defect orientation and separation [Figure 6(c)]. Thus, this vacancy containing defect will not play any role in the ferromagnetism in N-doped graphene.

III.3 Effects of boron co-doping

Motivated by the recent observation of room-temperature ferromagnetism in boron and nitrogen doped turbostratic carbon films, Mondal et al. (2018) we further investigate the effect of B co-doping on the ferromagnetism. Notably, we consider a small B co-doping in the limit of high N-concentration that is responsible for the observed ferromagnetism. Due to very-low B-concentration, we have only considered graphitic single B-defect and ignored the possibility of the larger B-complexes. We also ignored the chemisorbed B-defect, as in this configuration, the B-adatom is expected to be very mobile similar to the C and N adatoms (Figure 3), and consequently form single graphitic B-defect.

In contrast to the single N-defect, the graphitic B-defect leads to p-type behaviour due to electron deficiency. The C–B bonds are found to elongated (1.49 Å) compared to the C–N bonds (1.41 Å) in the graphitic N-defect and a B $\rightarrow$ C charge transfer takes place. However, similar to the single graphitic N-defect, the B-defect also results in an itinerant magnetism with 1 $\mu_{B}$ moment, where a very small 0.03 $\mu_{B}$ moment is localized at the B-site. In contrast to the repulsive and sublattice dependent N$-$N magnetic interactions, the N$-$B interaction is found to be attractive, and a completely compensated AFM structure emerges, which is independent of the sublattice. Further, we investigate the interaction between this graphitic B-defect and the prevalent graphitic N-defects, T and G3N complexes with 1 $\mu_{B}$ moment each. Surprisingly, the BA/B affects the overall picture significantly. The interaction between these graphitic defects and the BA/B is found to be attractive, which is opposite to their interaction with single graphitic N-defect. Independent of their respective sublattice, the interacting TAAA$-$BA/B [Figure 7(h)–(i)] and G3NAAA$-$BA/B configurations generate a total moment of 2 $\mu_{B}$, while for all cases, the compensated AFM magnetic structure is found to be energetically much unfavourable. A closer look at the magnetization density reveals that the resultant 2 $\mu_{B}$ moment arises due to the sublattice imbalance of the itinerant C-$p_{z}$ electrons. This finding is consistent with the recent experimental observation of increased saturation magnetization in N-doped graphene due to B co-doping. Mondal et al. (2018) Therefore, we argue that small B co-doping will enhance ferromagnetism in nitrogen-doped graphene.

IV Summary

The microscopic origin of ferromagnetism that is observed in N-doped graphene has remained elusive. (Ito et al., 2015; Blonski et al., 2017; Liu et al., 2016; Li et al., 2015; Miao et al., 2016; Khurana et al., 2013) The difficulty to comprehend arises due to the complex dependence of magnetism on N-content and on the corresponding differential defect distribution originating from the varied growth and annealing details. In this context, we systematically investigate the moment formation and the intricate magnetic interaction between various defect complexes within the first-principles calculations. The C and N adatom diffusion on graphene plays a very important role in healing vacancy containing defects and generating graphitic N-complexes. The magnetism in the graphitic N-complexes are found to be fundamentally different from the vacancy containing complexes. While an itinerant magnetism emerges in the graphitic N-complexes, a very significant fraction of the moment is localized at the N-sites for the later. By comparing the results with PBE and HSE06 exchange-correlation functional and available experimental data, we conclude that the SCAN semilocal meta-GGA density functional provides a better description of the electronic structure and magnetism in defected graphene. We further demonstrate that the vacancy magnetism is robust against carrier doping.

The magnetic interaction between the defects are found to be intricate and complex, which depends on the defect type, N-concentration, defect orientation and the corresponding distance between them. Regardless of the defect type and at low N-concentration, the magnetic interaction shows a generic sublattice dependence independent of defect orientation and distance. In agreement with the RKKY interaction leading to sublattice dependent magnetic interaction in bipartite graphene lattice, (Saremi, 2007; Brey et al., 2007) the long-range coupling is found to be FM (AFM) for the defects on the equivalent (opposite) sublattice. Thus, at low N-content, the statistical distribution of defects on both graphene sublattice will not produce FM ordering. This result is consistent with experimental observation. (Ito et al., 2015) Interestingly, the picture drastically changes at the high N-concentration owing to a direction dependent and sublattice independent FM interaction with spontaneous magnetization for the interaction between the graphitic T and G3N defects. While we reject the other type of defects through a comprehensive investigation, we conclude that graphitic T and G3N defects are predominantly responsible for the observed ferromagnetism. Further, we predict that an experimental strategy to increase the concentration of these graphitic defects will produce robust ferromagnetism with improved Curie temperature. Finally, we demonstrate that a small amount of boron co-doping will further enhance the ferromagnetism in N-doped graphene. Thus, while the present results help us to understand the scattered and often controversial experimental results, an overall microscopic understanding and a concurrent experimental synthesis strategy provided here will motivate new experiments to control magnetism in nitrogen and boron doped graphene.

Acknowledgements.

We acknowledge stimulating discussions with S. B. Ogale and R. Mandal. M.K. acknowledges funding from the Science and Engineering Research Board through Ramanujan Fellowship and the Department of Science and Technology through Nano Mission project SR/NM/TP-13/2016. The supercomputing facilities at the Centre for Development of Advanced Computing, Pune; Inter University Accelerator Centre, Delhi; and at the Center for Computational Materials Science, Institute of Materials Research, Tohoku University are greatly acknowledged.