Stable local moments of vacancies, substitutional and hollow site impurities in graphene

Since its isolation in the laboratory, graphene has stimulated extensive research activity among condensed matter and material physicists [1–3]. The nature of low-energy excitations, which are described by the Dirac equation, and its two-dimensional structure, are distinct features of graphene [4], which are responsible for many of its novel properties. For example, two-dimensionality allows for tuning the carrier type and density by a gate voltage. On the one hand, the flat structure of graphene allows us to add/substitute various atoms on graphene [5]. On the other hand, the two-dimensional structure further allows us to conveniently remove some of the carbon atoms and create vacancies.

Several ab initio calculations have been done on graphene to explore point defects and impurities [6–18]. As an example, even before the isolation of a single layer of carbon atoms from graphite by the Manchester group, first-principle calculations, by using spin-polarized density functional theory and implementing the generalized gradient approximation, revealed that the presence of carbon adatoms on graphene leads to magnetic local moment ∼0.5 μ_B [10]. Furthermore, the formation of magnetic states from different adatoms and vacancies have been investigated in graphene and carbon nanotubes [11–16]. Ab initio studies showed that the properties of local magnetic moments depend on the gate voltage, position of the impurity and curvature of the carbon layer [6–9]. Interestingly, experimental reports revealed that vacancies and adatoms were responsible for spin-half magnetism in graphene [19]. The magnitude of local moments, which originate from vacancy (p-wave) and top site (s-wave) impurity, can be tuned experimentally by electron (hole) doping [20] in agreement with first-principle studies [8, 17, 18]. This experiment emphasizes the difference between local moment formation for the vacancy and top site impurity. The potential distortion resulting from the creation of a vacancy leads to a shallow impurity state which hybridizes with orbitals of three neighboring carbon atoms. Substitutional and hollow site adatoms have similar symmetry and non-local hybridization with neighboring carbon atoms. The possibility of such a non-local hybridization as opposed to top site adatoms, after simple Fourier transformation from localized Wannier states to corresponding Bloch wave-functions results in a momentum dependence in the hybridization matrix element between the state localized on vacancy and the π-bands of graphene. After linearization around the Dirac points at low energies, the mixing function acquires a p-wave dependence on momentum [21, 22],

$$\begin{equation} V({\vec k})=(k_x-{\rm i}k_y) \tilde V. \end{equation} \tag{ 1 }$$

The above form is consistent with the fact that the pseudo-spin structure of the Bloch wave functions in graphene allows them to hybridize with external states in both l = 0 and l = 1 angular momentum channels [22, 23]. To be precise, the same form holds for the hybridization between the localized states of a substitutional adatom as illustrated in figure 1. Formally, the same result applies to the hybridization of Bloch states of graphene with a localized state of a hollow site adatom. The only small difference in the latter case is that the impurity state in the hollow site case hybridizes with both sublattices through $V({\vec k})$ and $V^*({\vec k})$ functional forms [24]. But the contribution of the hybridization with the two sublattices to the Hartree self-energy is identical. Therefore, the problem of local moment formation in hollow-site adatoms is not formally different from the vacancy or substitutional adatom. Once we work out the case of a vacancy or a substitutional adatom, to obtain a result corresponding to the hollow-site adatom it is sufficient to make the replacement $V\to \sqrt 2 V$. In this paper, 'vacancy' will represent all the above family of p-wave mixing cases while the top site stands for s-wave mixing. This new possibility of a p-wave mixing does not exist in normal metallic hosts where the wave function is described by a scalar function. The combination of the spinorial (two-component) nature of the Dirac wave functions in graphene with the C₃ symmetry is responsible for the new opportunity of mixing in the l = 1 angular momentum channel between the localized state and the Bloch states.

Zoom In Zoom Out Reset image size

In this work, we will focus on the local moment formation arising from the p-wave hybridization and will compare the local moment properties with the s-wave hybridization case. The phase diagram of magnetic states in graphene has been previously obtained by Uchoa and co-workers [25]. They studied the channel l = 0 corresponding to the momentum-independent hybridization function. Within the Hartree mean-field theory, they obtained a phase diagram for the magnetic states that was significantly different from the corresponding one for ordinary metals [26]. In this work, having the material's application in mind, we focus on the momentum-dependent hybridization function of the p-wave form relevant to vacancies, substitutional and hollow site adatoms and investigate the local spectral properties of the single impurity Anderson model (SIAM) both within and beyond the Hartree mean-field approximation. Within the mean-field theory, in contrast to the momentum-independent case [25], we find that with the p-wave mixing, the magnetic states are formed in a much larger region of the parameter space. This could be interpreted as the parametric robustness of the local moments arising from p-wave hybridization. Then, within the equation of motion approach, we proceed one step beyond the Hartree decoupling scheme and find that the local moments arising from the p-wave hybridization display stability against quantum fluctuations. These findings shed light on the experimentally observed formation of spin-half states which are attributed to vacancies in graphene [19] and the observed difference in the behavior of local moments between vacancy and top site impurities [20].

The Anderson model for graphene with a vacancy can be written as,

$$\begin{equation} H = {H_{\rm g}} + {H_{\rm hyb}} + {H_{\rm imp}}. \end{equation} \tag{ 2 }$$

The first term of the above Hamiltonian describes π electrons of graphene. The second term describes hybridization (mixing) of the impurity state with Dirac fermions in graphene. To create a vacancy we assume that a carbon atom from A-sublattice has been removed. In principle, the removal of a vacancy breaks the translational symmetry of the underlying graphene host. But neglecting the loss of translational symmetry introduces a small effect of the order of the relative concentration of impurities. On the other hand, a vacancy can be thought of as a strong local potential U₀ the effective result of which is the formation of a resonance state around the vacant site. Therefore, an effective way to model a vacancy is to model it with a very shallow localized state around the vacant site. These terms after a standard Hartree mean-field factorization will be written as

$$\begin{eqnarray} &&\begin{array}{l} {H_{\rm g}} = - t\sum\limits_{{\vec k}\sigma} {(\phi( {\vec k})a_{{\vec k}\sigma}^\dagger {b_{{\vec k}\sigma}} + {\phi ^*}( {\vec k} )b_{{\vec k}\sigma}^\dagger {a_{{\vec k}\sigma}})} \\ {H_{\rm hyb}} = \frac{1}{\sqrt N}\sum\limits_{{\vec k}\sigma} V({\vec k}) b_{{\vec k}\sigma}^\dagger {d_\sigma} + {\rm h.c.}\\ {H_{\rm imp}} = \sum\limits_\sigma {{\varepsilon _\sigma}d_\sigma^\dagger {d_\sigma}} ~~~~;{\varepsilon _\sigma} \equiv {\varepsilon_0} + U\langle {{n_{\bar \sigma}}} \rangle \end{array} \end{eqnarray} \tag{ 3 }$$

where $\phi({\vec k})=\sum_{n=1}^3\exp(-{\rm i}\vec \delta_n.{\vec k})$ is the form factor associated with three neighbors connected by vectors δ_n, with n = 1, 2, 3. The creation operators $a^\dagger_{{\vec k}\sigma}$ and $b^\dagger_{{\vec k}\sigma}$ create Bloch electrons on A and B sublattices with spin σ, respectively and $d^\dagger_\sigma$ creates an electron with spin σ in the impurity state associated with the vacancy. Here ε₀ is the impurity level, which after incorporating the Hubbard repulsion U splits into two sub-levels denoted by ε_σ. The hybridization function $V({\vec k})$ as emphasized in equation (1) is proportional to the form factor $\phi({\vec k})$ which in the continuum limit becomes proportional k_x − ik_y. After the Hartree mean-field factorization implicit in equation (3), the impurity Green's function can be written as,

$$\begin{equation} G_{d\sigma}^{-1} (\omega)= {\omega - {\varepsilon _\sigma} - \Sigma^{\Delta}_d(\omega)} \end{equation} \tag{ 4 }$$

where,

$$\begin{equation} \Sigma^{\Delta}_d(\omega) = \sum_{\vec k} |V({\vec k})|^2 \frac{\omega}{\omega^2-t^2|\phi({\vec k})|^2} \end{equation} \tag{ 5 }$$

The above integral can be easily evaluated to give,

$$\begin{eqnarray} &&\fl - \frac{{2\omega\Delta}}{{\pi {t^2}}}\left[ {{D^2} + {\omega ^2}\ln \bigg( {\frac{{\vert {{\omega ^2} - {D^2}} \vert}}{{{\omega ^2}}}} \bigg) + {\rm i}\pi \omega \vert \omega \vert\theta ( {D - \vert \omega \vert} )} \right] \nonumber\\ \end{eqnarray} \tag{ 6 }$$

where $\Delta \equiv \frac{{\pi {V^2}}}{{{D^2}}}$, and D ∼ 7 eV is the high-energy band cutoff chosen such that the total number of electrons in the linearized band are the same as the π-band. Note that in contrast to the momentum-independent case—see equation (23)—where the imaginary part of the $\Sigma^\Delta_d$ is proportional |ω|, here the p-wave momentum-dependence of the hybridization gives rise to anomalous |ω|³ dependence in the imaginary part of $\Sigma^\Delta_d$. Inserting the above equation in equation (4) the impurity Green's function becomes,

$$\begin{equation} G_{d\sigma}(\omega) = \frac{1}{{{Z^{-1}}( \omega)\omega - {\varepsilon_\sigma} +{\rm i}|\omega|^3\Delta \theta(D-|\omega|)/t^2}}, \end{equation} \tag{ 7 }$$

where Z⁻¹ (ω) is given by,

$$\begin{equation} {Z^{- 1}}( \omega ) = 1 + \frac{{{2V^2}}}{{{t^2}{D^2}}}\left[ {{D^2} + {\omega ^2}\ln ( {{{\vert {{\omega ^2} - {D^2}} \vert}/ {{\omega ^2}}}} )} \right]. \end{equation} \tag{ 8 }$$

The imaginary part of ${\Sigma^{\Delta}_d}(\omega)$ shows the broadening of the localized level due to hybridization. Here, in contrast to normal metals, the broadening displays strong ω-dependence. The ω-dependence in the p-wave case is even stronger than the s-wave case. Then the local spectral function becomes,

$$\begin{equation} A_{d\sigma} (\omega) = \frac{\Delta}{{\pi {t^2}}} \frac{{{{\vert \omega \vert}^3}\theta ( {D - \vert \omega \vert} )}}{{{{( {{Z^{- 1}}( \omega )\omega - {\varepsilon _\sigma}} )}^2} + {{{\Delta ^2}{\omega ^6}} / {{t^4}}}}} \end{equation} \tag{ 9 }$$

which in turn gives the occupation of the σ sub-band as,

$$\begin{equation} {n_\sigma} = \frac{\Delta}{{\pi {t^2}}}\int_{- \infty}^\mu {{\rm d}\omega \frac{{{{\vert \omega \vert}^3}\theta ( {D - \vert \omega \vert} )}}{{{{( {{Z^{- 1}}( \omega )\omega - {\varepsilon _\sigma}} )}^2} + {{\Delta ^2}{\omega ^6}} / t^4}}} \end{equation} \tag{ 10 }$$

In this equation, the n_↑ is given as an integral involving ε_↑ which itself depends on n_↓ as given by equation (3). These equations for n_↑ and n_↓ must be solved self-consistently. Equation (9) must be compared to equation (9) of [25] corresponding to s-wave hybridization of top site impurities. In order to make the comparison, we have introduced function Z⁻¹(ω) along with the notations of this reference. The function Z⁻¹(ω) in our case differs from the ${\vec k}$-independent hybridization in two respects: (i) The logarithmic part of the ω-dependence has acquired an anomalous ω² factor compared to the ${\vec k}$-independent case. (ii) Comparing $\Sigma^\Delta_d$ with that of the s-wave case,—equation (23)—in the p-wave case an additive D² term appears. To complete the self-consistency cycle, one needs to accurately evaluate the integrals connecting n_↑ and n_↓. The authors of [25] argue that neglecting the ω-dependence of Z⁻¹(ω) allows for analytic evaluation of the integral at the cost of introducing a few percent error. But as we will show in the following, since such errors are repeated several times throughout the self-consistency cycle, the final results are sensitive to the ω-dependence of the Z⁻¹(ω) function.

In this section, we present numerical results for the self-consistent solution at the level of Hartree mean-field equations. To begin with, we revisit the top site hybridization case discussed in [25] and show how the approximation of neglecting the ω-dependence of Z⁻¹(ω) leads to a different phase diagram. Then we focus on the case of the p-wave hybridization function pertinent to vacancy and hollow site impurities.

3.1. Local moments from top site impurities

In this case, the hybridization $V_{\vec k}$ does not depend on ${\vec k}$, and the occupation of spin σ impurity sub-band is given by [25],

$$\begin{equation} n_{\sigma}=-\frac{1}{\pi}\int^{\mu}_{-D}{\rm d}\omega \frac{\Delta\vert\omega\vert\theta(D-\vert\omega\vert)}{[Z^{-1}(\omega)\omega-\varepsilon_{\sigma}]^2+\Delta^2\omega^2}, \end{equation} \tag{ 11 }$$

where

$$\begin{equation} Z^{-1}(\omega)= 1 + \frac{V^2}{D^2} \ln\bigg(\frac{\vert D^2 - \omega^2 \vert}{\omega^2}\bigg). \end{equation} \tag{ 12 }$$

Here we do not apply the approximation Z⁻¹(ω) ≈ Z⁻¹(ε_σ) and retain the full energy dependence of Z⁻¹. The price one has to pay is to do the integrations necessary to get n_σ numerically. We define two standard variables X = DΔ/U and Y = (μ − ε₀)/U in terms of which the phase diagram is traditionally constructed [26].

The result is shown in figure 2. In this plot the colored region represents the intensity plot of the magnetic region 〈n_↑〉 ≠ 〈n_↓〉. In this plot we retain the full energy dependence of Z⁻¹(ω) and the integrals are evaluated numerically. To compare our numerically derived phase diagram with approximate results where the energy dependence of Z⁻¹ (ω) has been neglected, we demonstrate the difference between these two cases in figure 3. As can be seen, there are differences between the two schemes. The numeric result compared to the approximate one shows that the above-mentioned approximation underestimates the magnetic moment. In addition, with small values of parameter X and close to regions Y = 1, this change in local moment strength becomes as large as 80%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although neglecting the energy dependence in Z⁻¹(ω) leads only to small error in each evaluation of the integral, repeating this procedure in the self-consistency process, propagates and enhances the errors. With this point in mind, in the following we avoid using such an approximation and focus on the self-consistent solution of the mean-field equations for the p-wave hybridization function.

3.2. Local moments from vacancies or hollow site impurities

Now we proceed by self-consistently solving the mean-field equation (9) for vacancy-induced local moments. In figure 4 we compare the self-consistently determined values of occupation numbers n_↑ and n_↓ for a fixed value of X = 5.0 in the two cases corresponding to top site impurity and vacancy. For very small and very large values of Y the n_↑ and n_↓ curves coincide which means the magnetic moment is zero. For intermediate values of Y the curves corresponding to n_↑ and n_↓ split. The amount of splitting determines the magnetic moment. As can be seen, the local moment formation region in vacancy sets in smaller values of Y on the negative values. Moreover, in the case of vacancy the magnetic moment inside the local moment region stays more or less constant and abruptly vanishes by approaching the boundary, while in the top site case, the local moment smoothly vanishes as the boundary is reached.

Zoom In Zoom Out Reset image size

Repeating the above analysis for other values of X, we obtain the phase diagram depicted in figure 5. This plot has been generated for V/D = 0.14 and two different values of ε₀/D = 0.029 and ε₀/D = 0.043 as indicated in the figure. The plots corresponding to ε₀/D = 0.043 in this figure can be compared to the one in figure 2 of the top site impurities as they correspond to the same set of parameters. As can be noted by comparing figures 2 and 5, in the case of p-wave hybridization, the local magnetic moments are found in a much larger region of parameter space. This can be interpreted as the parametric robustness of local moments arising from vacancies or hollow site impurity states compared to those due to top site impurity states. Another important aspect can be noticed by comparison of the p-wave phase diagram to the numerically exact plot in figure 2 of the s-wave case is that in the latter case the local moment region dominantly lies in the positive values of Y, a feature shared by normal metals, while in the former case, the local moment region is substantially extended to negative values of Y and has been elongated as a stripe-shaped region. Such elongation has been reported in a similar study for bilayer graphene [27] and has been shown to arise from a large contribution in the real part of the self-energy in the vacant case which arises from a constant shift in the scale of bandwidth D. Such self-energy shifts can be absorbed into ε₀ defining a largely shifted effective $\tilde \varepsilon_0$. This is consistent with the fact that small changes in ε₀ do not change the elongated nature of the phase diagram. Moreover, up to quite large values of X ∼ 12—which are already much larger than physically conceivable values for physical graphene—the width of the range of Y values for which the local moment is formed does not vary much. The larger magnetic moment formation area in the phase space can be interpreted as the robustness of the local moments formed by vacancies compared to top site impurities, and is consistent with the observed spin-half magnetic moment in vacant graphene [19].

Zoom In Zoom Out Reset image size

3.3. Internal structure of the impurity orbital in hollow site adatoms

Large adatoms prefer the hollow site position when added to a graphene sheet [7] while small atoms such as hydrogen prefer the top site position [5]. Therefore, implicit assumption of s-wave orbital structure for the impurity atom itself is only feasible for top site positions. With larger adatoms in the hollow site positions, higher multiplicities and more complex symmetry patterns may appear. In this section, we address the question of local moment formation for impurity orbitals that preserve the time reversal symmetry within the Hartree mean-field theory. For the case of hollow site adatoms, the s-wave symmetry is formally identical to what we considered in the vacancy case and that is why we treated these two cases on the same basis. The only difference between them is that the $\Sigma^\Delta_d$ in the s-wave hollow site case is twice that of the vacancy case. The reason is simply related to the fact that now in addition to three nearest neighbors from a given sublattice, three nearest neighbors from the other sublattice also emerge. Hence the phase diagram of local moments in the case of s-wave adatoms in the hollow site position is identical to that of vacancy after a re-scaling of hybridization strength V, i.e.: $V^{\rm hollow~site}=\sqrt 2 V^{\rm vacancy}$.

The p orbitals alone can not form a linear combination compatible with the six-fold rotational symmetry of graphene. The linear combination of d orbitals compatible with the symmetries of honeycomb hexagon requires complex coefficients which break the time reversal symmetry. Therefore, the simplest non-trivial higher angular momentum pattern compatible with the rotational symmetry of the honeycomb lattice is the f-wave pattern. For this purpose we consider in the hollow site configuration a hybridization Hamiltonian of the following form,

$$\begin{eqnarray} &&{H_{\rm hyb}} = \frac{V}{\sqrt N}\sum_{{\vec k}\sigma} {( \phi({\vec k}) b_{{\vec k}\sigma}^\dagger + {\rm e}^{{\rm i}\alpha}{\phi ^*}({\vec k}) a_{{\vec k}\sigma}^\dagger ) d_\sigma + {\rm h.c.}} \end{eqnarray} \tag{ 13 }$$

In the above Hamiltonian, when the parameter α takes the value π it corresponds to the f-wave pattern of the atomic orbital. Obviously α = 0 reduces to the s-wave atomic orbital where the impurity is hybridized with both A and B sublattices equivalently, while in the f-wave case the hybridization pattern changes the sign of V after each 2π/6 rotation.

The equation of motion for the Green's function of the above Hamiltonian will now depend on the phase α,

$$\begin{equation} G_{d\sigma}^{-1}(\omega) = \omega - \varepsilon_{0} - \Sigma^{\Delta}_{d}(\omega) \end{equation} \tag{ 14 }$$

where the dependence on α is through the self energy:

$$\begin{eqnarray} &&\fl \Sigma^{\Delta}_d(\omega) = {V^2}\sum\limits_{\vec k} \frac{\omega}{\omega^2-t^2|\phi({\vec k})|^2}(g({\vec k},\omega)+{\rm e}^{-2{\rm i}\alpha}g^*({\vec k},\omega)),\nonumber\\ &&\fl g^*({\vec k},\omega) \equiv \frac{{v_{\rm F}^2{k^2}}}{{{t^2}}}( {1 - {{{v_{\rm F}}k{{\rm e}^{{\rm i}(3\theta + \alpha )}}} / \omega}} ). \end{eqnarray} \tag{ 15 }$$

For both α values of zero and π, we have e^−2iα = 1. Furthermore, the θ integration of e^i(3θ+α) produces zero. Therefore, the self energy at the Hartree level will not depend on α, and the difference between the s-wave and f-wave internal structures of the impurity orbital does not appear at the mean-field level.

So far we have investigated within the mean-field level the formation of local moments for two situations corresponding to the s-wave and p-wave hybridization function $V({\vec k})$. We have further checked that in the hollow site case the internal structure of the impurity orbital itself does not matter within the mean-field theory. When the local moments are formed in the Hartree mean-field level, the question will be, what happens to the local magnetic moments when quantum fluctuations beyond the mean field are taken into account? A proper treatment of the effect of fluctuations leads to the dynamical screening of the magnetic moment (Kondo screening). In this section, we address this question by going one step beyond the Hartree mean field within the equation of motion approach. To do this we use the standard notation for the Fermionic correlation functions:

$$\begin{equation} \langle{\langle{d_\sigma(t)|d_\sigma^\dagger(t^\prime)}\rangle}\rangle = -{\rm i}\theta(t-t^\prime) \langle{\lbrace d_\sigma(t)|d_\sigma^\dagger(t^\prime)\rbrace}\rangle. \end{equation} \tag{ 16 }$$

Writing the equation of motion for the above Green's function in frequency domain gives,

$$\begin{equation} (\omega - \varepsilon _0 - \Sigma^{\Delta}_d(\omega))\langle{\langle{{d_\sigma}|d_\sigma ^\dagger}\rangle}\rangle = 1 + U\langle{\langle{{d_\sigma}{n_{\bar\sigma}}|d_\sigma ^\dagger}\rangle}\rangle. \end{equation} \tag{ 17 }$$

In the Hartree approximation the correlation function on the right of equation (17) is approximated by $\langle{n_{\bar\sigma}}\rangle\langle{\langle{d_\sigma}|d_\sigma ^\dagger}\rangle\rangle$, which closes the equations, and the effect of the Hubbard U will be to replace $\varepsilon_0\to\varepsilon_\sigma=\varepsilon_0+U \langle n_{\bar\sigma}\rangle$. But here we do not close the equations at this level and proceed one step further by writing another equation of motion for the correlation function $\langle{\langle{d_\sigma}{n_{\bar\sigma}}|d_\sigma ^\dagger}\rangle\rangle$. This gives,

$$\begin{equation} \omega \langle{\langle{d_\sigma}{n_{\bar \sigma}|d_\sigma ^\dagger}\rangle}\rangle = \langle \lbrace{d_\sigma} {n_{\bar \sigma}},d_\sigma ^\dagger\rbrace\rangle + \langle{\langle[{d_\sigma}{n_{\bar \sigma},H]|d_\sigma ^\dagger}\rangle}\rangle. \end{equation} \tag{ 18 }$$

Let us first focus on the simplest case which is the top site hybridization. In this case, calculating the (anti)commutation relations needed in the above equation we find:

$$\begin{eqnarray} &&({\omega - {\varepsilon _0} - U})\langle {\langle {{d_\sigma}{n_{\bar \sigma}}|d_\sigma ^\dagger}\rangle} \rangle = \langle {{n_{\bar \sigma}}}\rangle \nonumber\\ &&\qquad + V\sum\limits_{\vec k}\langle\langle {{b_{{\vec k}\sigma}}{n_{\bar \sigma}} + {d_\sigma}({d_{\bar \sigma}^\dagger {b_{p\bar \sigma}} - b_{p\bar \sigma}^\dagger {d_{\bar \sigma}}})|d_\sigma ^\dagger} \rangle\rangle. \end{eqnarray} \tag{ 19 }$$

At this stage, the decoupling scheme can be applied to the correlation function on the right hand side of equation (19) as follows:

$$\begin{eqnarray} &&\fl \langle\langle {{b_{{\vec k}\sigma}}{n_{\bar \sigma}} + {d_\sigma}({d_{\bar \sigma}^\dagger {b_{p\bar \sigma}} - b_{p\bar \sigma}^\dagger {d_{\bar \sigma}}})|d_\sigma ^\dagger} \rangle\rangle \nonumber\\ &&\approx \langle{n_{\bar \sigma}}\rangle \langle\langle {b_{{\vec k}\sigma}}|{d_\sigma ^\dagger}\rangle\rangle + \langle{d_{\bar \sigma}^\dagger {b_{p\bar \sigma}} - b_{p\bar \sigma}^\dagger {d_{\bar \sigma}}}\rangle\langle\langle {d_\sigma}|{d_\sigma ^\dagger} \rangle\rangle. \end{eqnarray} \tag{ 20 }$$

The above decoupling has been performed in such a way to avoid off-diagonal correlation functions in the spin-indices. Hence the operators carrying $\bar{\sigma}$ spin index are taken out of the correlations in the form of expectation values [28]. Finally, to close the set of equations of motion, we write the equation of motion for $\langle\langle {b_{{\vec k}\sigma}}|{d_\sigma ^\dagger}\rangle\rangle$, which gives,

$$\begin{equation} \langle {\langle {{b_{p\sigma}}|d_\sigma ^\dagger} \rangle} \rangle = V\frac{\omega}{\omega^2-t^2|\phi({\vec k})|^2}\langle {\langle {{d_\sigma}|d_\sigma ^\dagger} \rangle} \rangle. \end{equation} \tag{ 21 }$$

Therefore, at the present approximation, the local Green's function becomes,

$$\begin{equation} \langle {\langle {{d_\sigma}|d_\sigma ^\dagger} \rangle} \rangle = {\left[ {\omega - {\varepsilon _0} - \Sigma^{\Delta}_d ( \omega ) - \Sigma ^{\prime}( \omega )} \right]^{- 1}}, \end{equation} \tag{ 22 }$$

where $\Sigma^\Delta_d$ describes the hybridization of local electrons with the band continuum and for the top site situation is given by [25]

$$\begin{eqnarray}&&\fl \Sigma^\Delta_d(\omega)=-V^2\frac{\omega}{D^2} \ln\left(\frac{|\omega^2-D^2|}{\omega^2} \right) -{\rm i}V^2\frac{\pi|\omega|}{D^2}\theta(D-|\omega|).\nonumber\\ \end{eqnarray} \tag{ 23 }$$

In the p-wave case, as is shown in the appendix, within the present approximation, the form of the local Green's function remains the same as for equation (22). The only difference is that now the hybridization function $\Sigma^\Delta_d$ for the hollow site case is given by equation (6). Since the Hubbard U acts only locally in the impurity orbital, at this level of approximation, the self-energy correction coming from the hybridization $\Sigma^\Delta_d$, and those coming from the Hubbard term Σ', are separable and hence the latter turns out to be insensitive to the nature of hybridization of the impurity orbital with neighboring orbitals. The interaction effects arising from the Hubbard U are given by the self-energy Σ'(ω) as follows:

$$\begin{equation} \Sigma ^{\prime}(\omega) \equiv \frac{{U({\omega - {\varepsilon _0}})\langle {{n_{\bar \sigma}}} \rangle}}{{\omega - {\varepsilon _0} - U( {1 - \langle {{n_{\bar \sigma}}} \rangle})}}. \end{equation} \tag{ 24 }$$

Note that in this case, unlike the Hartree approximation where the mere effect of the Hubbard U is to shift ε₀ by $U\langle n_{\bar\sigma}\rangle$, in the present approximation, the self-energy not only is not of a simple shift form, but has acquired a non-trivial ω-dependence. Note that truncation of equation of motion is not a well-controlled approximation as there is no small parameter to control the accuracy. But it has been shown by Lacroix that by keeping the hierarchy of equations of motion up to the third level it is capable of capturing the Kondo effect [29]. Moreover, we compare the s- and p-wave hybridizations within this same approximation, irrespective of the accuracy of the present truncation.

In figure 6 we have compared the self-consistent solutions for the occupation of the ↑ and ↓ spin sub-bands within the Hartree and our post-Hartree approximation for the top site impurity. The plot has been generated for a fixed value of X = 5.0 and ε₀ = 0.029 as function of Y. Within the Hartree approximation, plots corresponding to n_↑ and n_↓ are split in magnetic region. This indicates the formation of local moments at the Hartree level. However, as soon as we employ the self-energy (24) to go beyond the Hartree, the n_↑ and n_↓ plots collapse on each other. Therefore, the local moments formed at the Hartree level are destroyed by quantum fluctuations beyond the mean-field level.

Zoom In Zoom Out Reset image size

In figure 7 we plot the same set of data as in figure 6 but for the p-wave hybridization function. The interesting feature to note in this case is that, even by including the self-energy effect equation (24) which goes a step beyond the Hartree approximation, there will still be a region in the parameter space where the local moment survives the fluctuations. Therefore, the local moments formed by hollow site impurities (or vacancies) when compared to those due to top site impurities, are not only more robust in terms of the extension of the region in the parameter space where local moment is formed, but they are robust against the quantum fluctuation effects captured by equation (24). This makes vacancies and hollow site impurity states ideal sources of local magnetic moments in graphene which are hard to destroy. This can account for the observed magnetic moments in vacant graphene [19].

Zoom In Zoom Out Reset image size

In this paper, we investigated the formation of local moments in a graphene sheet, in two situations corresponding to s-wave and p-wave functional dependence in the hybridization function $V({\vec k})$ within the single impurity Anderson model. The former corresponds to the top site impurities, while the latter may correspond to hollow site and substitutional impurities or vacancies. First, we noted that the output of the self-consistency cycle is sensitive to the precision of the integrations. Any approximation in the integrations propagates the errors through the self-consistency cycle and may give rise to different results. Then, we compared the Hartree mean-field level phase diagram of the above two cases. We found that vacancies and hollow site impurities give rise to a much larger region in the phase space where the magnetic moments are formed. In this sense the local moments arising from impurity's orbitals on the hollow sites or from vacancies are robust against change in the parameters. Furthermore, we showed that within the Hartree mean field the internal structure of the impurity orbital itself does not affect the local moment properties. Then, within the equation of motion method, we proceeded one step beyond Hartree approximation and derived a self-energy due to the Hubbard U term which takes into account quantum fluctuations beyond the Hartree mean field. Within such a post-Hartree mean field we obtained the self-consistently determined occupation numbers of the ↑ and ↓ spin impurity sub-bands. We found that in the top site case, the quantum fluctuations in the post-Hartree approximation can destroy the local moments formed at Hartree level, while in the case of vacancies and hollow site impurities, for the same set of parameters, the local moments survive the quantum fluctuations within our approximation. Therefore, we conclude that the p-wave hybridization functions give rise to local moments which are not only robust in parameter space, but are immune to quantum fluctuations beyond the Hartree mean field. The present stability analysis is consistent with the observed spin-half magnetic states in vacant graphene [19] and sheds light on the differences in the local moment properties of the top site and more probable cases of p-wave mixing [20].

The authors thank T Tohyama for useful discussions. This research was completed while SAJ was visiting the Yukawa Institute for Theoretical Physics through the fellowship S13135 from the Japan Society for the Promotion of Science. MM thanks the Institute for Research in Fundamental Sciences, Tehran, Iran for hospitality during his visit to the School of Physics. MM also thanks the Abdus Salam International Centre for Theoretical Physics, Trieste, Italy for financial support to visit the center.

In this appendix, we calculate the impurity Green's function one step beyond the Hartree mean-field approximation for hollow site impurity. The equation of motion for the local Green's function is given by equation (17). Avoiding the decoupling of $\langle\langle n_{\bar{\sigma}}d_{\sigma}|d^\dagger_{\sigma}\rangle\rangle$ and writing an equation of motion for it we obtain,

$$\begin{equation} \omega \langle{\langle{d_\sigma}{n_{\bar \sigma}|d_\sigma ^\dagger}\rangle}\rangle = \langle {n_{\bar \sigma}}\rangle + \langle{\langle[{d_\sigma}{n_{\bar \sigma},H]|d_\sigma ^\dagger}\rangle}\rangle, \end{equation} \tag{ A1 }$$

which after evaluation of the necessary commutation relations becomes,

$$\begin{eqnarray} &&({\omega - {\varepsilon _0} - U})\langle {\langle {{d_\sigma}{n_{\bar \sigma}}|d_\sigma ^\dagger}\rangle} \rangle = \langle {{n_{\bar \sigma}}}\rangle \nonumber\\&&\qquad + V\sum\limits_{\vec k}\langle\langle {{c_{{\vec k}\sigma}}{n_{\bar \sigma}} + {d_\sigma}({d_{\bar \sigma}^\dagger {c_{p\bar \sigma}} - c_{p\bar \sigma}^\dagger {d_{\bar \sigma}}})|d_\sigma ^\dagger} \rangle\rangle. \end{eqnarray} \tag{ A2 }$$

In the above equation the operator $c^{\dagger}_{{\vec k}\sigma}$ is defined by

$$\begin{equation} c^{\dagger}_{{\vec k}\sigma} \equiv \phi({\vec k})b^{\dagger}_{{\vec k}\sigma} + \phi^*({\vec k})a^{\dagger}_{{\vec k}\sigma}. \end{equation} \tag{ A3 }$$

In the hollow site configuration the impurity hybridizes with carbon atoms from both sublattices. Therefore, the operator $c_{{\vec k}\sigma}^\dagger$ in equation (A2) plays the same role as $b_{{\vec k}\sigma}^\dagger$ in equation (19). We apply the same decoupling scheme which gives,

$$\begin{eqnarray} &&\fl \langle\langle {{c_{{\vec k}\sigma}}{n_{\bar \sigma}} + {d_\sigma}({d_{\bar \sigma}^\dagger {c_{{\vec k}\bar \sigma}} - c_{{\vec k}\bar \sigma}^\dagger {d_{\bar \sigma}}})|d_\sigma ^\dagger} \rangle\rangle \nonumber\\ &&\approx \langle{n_{\bar \sigma}}\rangle \langle\langle {c_{{\vec k}\sigma}}|{d_\sigma ^\dagger}\rangle\rangle + \langle{d_{\bar \sigma}^\dagger {c_{{\vec k}\bar \sigma}} - c_{{\vec k}\bar \sigma}^\dagger {d_{\bar \sigma}}}\rangle\langle\langle {d_\sigma}|{d_\sigma ^\dagger} \rangle\rangle. \end{eqnarray} \tag{ A4 }$$

Hence the Green's function at this approximation for the hollow site configuration will be,

$$\begin{equation} \left\langle {\left\langle {{d_\sigma}|d_\sigma ^\dagger} \right\rangle} \right\rangle = {\left[ {\omega - {\varepsilon _0} - \Sigma^{\Delta}_d \left( \omega \right) - \Sigma ^{\prime}\left( \omega \right)} \right]^{- 1}}. \end{equation} \tag{ A5 }$$

In the above equation the correlation effect has been encoded in self energy Σ' which in both p-wave and s-wave cases is given by equation (24). The only difference between the two possible hybridization functions appears in $\Sigma^\Delta_d$ which for the s-wave case is given by equation (23), while in the p-wave case it is given by equation (6).