RKKY Interactions in Graphene: Dependence on Disorder and Gate Voltage

RKKY Interactions in Graphene: Dependence on Disorder and Gate Voltage

Hyunyong Lee Division of Advanced Materials Science, Pohang University of Science and Technology (POSTECH), Pohang 790-784, South Korea    E. R. Mucciolo Department of Physics, University of Central Florida, Orlando, Florida 32816, USA    Georges Bouzerar Institut Néel, CNRS, départment MCBT, 25 avenue des Martyrs, BP 166, 38042 Grenoble Cedex 09, France School of Engineering and Science, Jacobs University Bremen, Bremen 28759, Germany    S. Kettemann Division of Advanced Materials Science, Pohang University of Science and Technology (POSTECH), Pohang 790-784, South Korea School of Engineering and Science, Jacobs University Bremen, Bremen 28759, Germany s.kettemann@jacobs-university.de
July 11, 2019
Abstract

We report the dependence of Ruderman-Kittel-Kasuya-Yoshida (RKKY) interaction on nonmagmetic disorder and gate voltage in graphene. First the semiclassical method is employed to rederive the expression for RKKY interaction in clean graphene. Due to the pseudogap at Dirac point, the RKKY coupling in undoped graphene is found to be proportional to . Next, we investigate how the RKKY interaction depends on nonmagnetic disorder strength and gate voltage by studying numerically the Anderson tight-binding model on a honeycomb lattice. We observe that the RKKY interaction along the armchair direction is more robust to nonmagnetic disorder than in other directions. This effect can be explained semiclassically: The presence of multiple shortest paths between two lattice sites in the armchair directions is found to be responsible for the reduced disorder sensitivity. We also present the distribution of the RKKY interaction for the zigzag and armchair directions. We identify three different shapes of the distributions which are repeated periodically along the zigzag direction, while only one kind, and more narrow distribution, is observed along the armchair direction. Moreover, we find that the distribution of amplitudes of the RKKY interaction crosses over from a non-Gaussian shape with very long tails to a completely log-normal distribution when increasing the nonmagnetic disorder strength. The width of the log-normal distribution is found to linearly increase with the strength of disorder, in agreement with analytical predictions. At finite gate voltage near the Dirac point, Friedel oscillation appears in addition to the oscillation from the interference between two Dirac points. This results in a beating pattern. We study how these beating patterns are effected by the nonmagnetic disorder in doped graphene.

I Introduction

A large effort has been devoted to the study the electronic transport properties of graphene due to the unusual nature of the quasiparticles in this material, which are massless chiral Dirac fermions. A recent experiment indicating that vacancies in graphene may induce local magnetic moments Chen et al. (2007) renewed the interest in investigating magnetic properties as well. It is belived that the carrier-mediated effective exchange coupling between local moments, or Ruderman-Kittel-Kasuya-Yoshida (RKKY) interaction may play a crucial role in establishing how magnetism develops in graphene. Probing these properties locally is not completely out of reach: Using spin-polarized scanning tunneling spectroscopy, a direct measurement of the RKKY interaction in a conventional metal has been done by measuring the magnetization curves of individual magnetic atoms adsorbed on a platinum surface.Meier et al. (2008) Similar experiments may soon be possible on graphene.

Analytical and numerical studies of the RKKY interaction in clean graphene at the neutrality point have been reported.Saremi (2007); Black-Schaffer (2010); Sherafati and Satpathy (2011a, b) In this context, a dominant feature of graphene is the bipartite nature of its honeycomb lattice. Due to particle-hole symmetry of bipartite lattices, the RKKY interaction always induces ferromagnetic correlation between the magnetic impurities on the same sublattice, but antiferromagnetic correlation between the ones on different sublattices.Saremi (2007) At the neutrality point, the dependence of the RKKY interaction on the distance between two local magnetic moments in graphene is found to be , in contrast to the standard behavior in conventional two-dimensional metallic systems where is expected. In doped graphene, but not too far from the neutrality point, the behavior changes to and two different length scales control the RKKY interaction: The wavelength corresponding to the inverse of the distance between the two Dirac points in momentum space, , and the Fermi wavelength .Sherafati and Satpathy (2011b)

Since the RKKY interaction is mediated by itinerant electrons in the host metal, nonmagnetic defects influence directly these interactions. On-site potential fluctuations scatter the phase of the electron’s wave function as well as cause random changes in its amplitude, altering any interference phenomenon observed in the clean system. The effects of weak nonmagnetic disorder on the exchange interactions in conventional metals have been carefully studied.de Chatel (1981); Jagannathan et al. (1988); Zyuzin and Spivak (1986); Lerner (1993); Bulaevskii and Panyukov (1986); Lee et al. (2012) Numerical work in the strongly localized regime has been done in one-dimensional disordered chains.Sobota et al. (2007) These studies found that in the diffusive regime the RKKY interaction, when averaged over disorder configurations, decays exponentially over distances exceeding the mean free path scale . However, its geometrical average value has the same amplitude as in the clean limit. It is only in the localized regime, on length scales , i.e., exceeding the localization length , that its geometrical averaged values are exponentially suppressed so that the RKKY interaction is cutoff over distances exceeding . In our previous work we have studied the effect of nonmagnetic disorder on the RKKY interaction in undoped graphene using the kernel polynomial method (KPM).Lee et al. (2012) The unexpected and interesting result found was that, at weak disorder, the RKKY interaction along the armchair direction is not influenced by nonmagnetic disorder as much as in the other directions.

Motivated by this unexplained behavior, in this paper we start by employing the semiclassical method to evaluate the RKKY interaction in graphene. This calculation helped us understand the -dependence and the origin of the direction dependence of the sensitivity to disorder seen in Ref. Lee et al., 2012. We then used the numerical KPM to calculate the RKKY interaction in graphene in order to study the effect of disorder and of doping in more detail.

In order to study the disordered conduction electrons in graphene numerically, we employ the single-band Anderson tight-binding model on a honeycomb lattice,

(1)

where is the hopping integral, is an annihilation operator of an electron at site , is the corresponding creation operator, is the uncorrelated on-site disorder potential distributed uniformly between , and indicates nearest-neighbor sites. Periodic boundary conditions are used and the lattice constant and are set to be unity in all calculations.

Below, we start with the derivation of the RKKY interaction expression in clean graphene using the semiclassical method.

Ii Semiclassical approach to the RKKY interaction

(a)                               
(b)                               
Figure 1: (Color online) Schematic diagrams of (a) the propagating paths of an electron in clean system and (b) the graphene lattice. The two sublattices are denoted as A and B and the two representative directions (zigzag and armchair) are indicated as dashed gray lines. is the angle between the displacement vector of the magnetic impurity and the unit vector .

Bergmann interpreted the RKKY oscillation as an interference of the conduction electron’s wave function scattered by the magnetic impurity and calculated the interference using the Huygens’ principle in a three-dimensional metal.Bergmann (1987) According to the Huygens’ principle in two dimensions,Cirone et al. (2002) the amplitude of a wave which propagates from a source at position decays with distance and gains a phase factor at a position given by

(2)

where is the amplitude of the wave at the source and the extra factor comes from the asymptotic form of the Bessel function in two dimensions (). If an electron propagates from to the origin in graphene (Fig. 1a), the amplitude gets a phase factor

(3)

where the wave vector is expanded around the two neighboring Dirac points and with a relative wave vector . During its propagation, the electron gets another phase factor, , where is the kinetic energy of the electron near the Dirac point, is the Fermi velocity and is the propagation time from to the origin. After the electron is scattered by a magnetic impurity at the origin, it travels back to the position . The amplitude then gains an additional modulation of the amplitude which is given by

(4)

where depends on the properties of the magnetic impurity at the origin. When the electron goes back to , its direction is opposite compared to the first propagation and this is the reason why the signs are opposite in the phase factors in and . In the time domain, however, it propagates in the same direction so that the phase factor related with the time () has the same sign. This is consistent with the diagrammatic expression for the RKKY interaction which has two retarded Green’s functions. For a closed path, a loop of the opposite direction makes a contribution with the same weight, so that the modulation of the charge density of an electron which has the energy is given by

(5)

such that the total charge modulation is given by

(6)

Here, for free electron, is the volume of the system, is the Fermi-Dirac distribution function, is the density of state at the Fermi level, is the pseudogap exponent, is the spatial dimension, and is the total time it takes to return to .

The total charge density in (Eq. (6)) can be easily related to the RKKY interaction.Bergmann (1987) Since in graphene and , the resulting distance dependence is , which is consistent with previous works.Saremi (2007); Black-Schaffer (2010); Sherafati and Satpathy (2011a) We conclude therefore that the existence of the pseudogap at the Dirac point causes the faster spatial decay of the RKKY interactions in graphene compared to conventional two-dimensional system (, ). Even though Eq. (6) does not give detailed information of about the RKKY interactions, such as the sign dependence on the sublattice and an extra phase factor depending on the direction between point and the origin, it agrees with previous results obtained by the Green’s function method,Sherafati and Satpathy (2011a); Kogan (2011) namely,

(7)
(8)

where is an angle between the position vector and the zigzag-AA direction ( as described in Fig. 1b.

Iii Numerical method

We start with a general expression for the RKKY exchange coupling constant in terms of the off-diagonal matrix elements of the local density of states ,Roche and Mayou (1999); Lee et al. (2012)

(9)

where , is the local coupling constant between the localized magnetic impurity and the itinerant electrons in the host metal, is the magnitude of the impurity spin, are the site indices of magnetic impurities located at position  (), and is the Fermi energy. Zero temperature () is assumed. Using the KPM, one may calculate the matrix elements of the local density of state efficiently,Weiße et al. (2006)

(10)

where is the th Chebyshev polynomial of the first kind, , is an electronic Hamiltonian, and are the Jackson kernels coefficients. The sum is taken up to a cutoff number . One may obtain the expansion coefficients using the recurrence relation of Chebyshev polynomials, namely, with and . Implicit in Eq. (10) is the normalization of the energy spectrum to a band of unity width.

Iv Disordered graphene at half filling (, )

We have investigated the effect of on-site nonmagnetic disorder at the charge neutrality point (). For each value of , 1600 different disorder configurations were generated and the corresponding matrix elements of the density of state (Eq. 10) evaluated through the KPM with a Chebyshev polynomial cutoff number on a lattice with sites. We have previously reported that in the diffusive regime (), which was characterized by determining the mean free path and the localization length , the armchair direction is the only direction in which the averaged value of the RKKY interaction over disorder configurations does not alter its sign.Lee et al. (2012) To illustrate this effect, we calculated the RKKY interactions along the directions () which pass through only the same sublattice and the results are shown in Fig. 2 together with the ones for the zigzag-AA and armchair-AA directions. Notice that the amplitudes of the RKKY interaction are multiplied by the cube of the distance in order to make the effect of nonmagnetic disorder more transparent. As expected from previous studies,de Chatel (1981); Jagannathan et al. (1988); Zyuzin and Spivak (1986); Lerner (1993); Bulaevskii and Panyukov (1986) the amplitude of the averaged interactions decays exponentially over length scales exceeding the mean free path ,

(11)

where indicates averaging over disorder configurations and is the amplitude of the RKKY interaction in a clean system. On-site disorder breaks the sublattice symmetry. Consequently, the sign of the RKKY interaction between the moments which are localized at the same sublattice fluctuates, allowing both ferromagnetic and antiferromagnetic correlation (Fig. 2a, b, c). The importance of the sublattice symmetry was highlighted in our previous work by considering the hopping disorder.Lee et al. (2012) When randomness is added to the hopping integral without on-site disorder, it does not break the sublattice symmetry and the RKKY interaction never changes its sign for magnetic moments sitting in either the same or different sublattices, even for a fixed disorder configuration. Interestingly, the averaged RKKY interaction amplitude along the armchair direction with on-site disorder does not change sign (Fig. 2d), while for a particular disorder configuration it randomly changes both sign and amplitude.

(a)     , zigzag-AA    
(b)                   
(c)                    
(d) , armchair-AA
Figure 2: (Color online) Plots of the RKKY interaction strength multiplied by the cube of the distance, , along the (a) (zigzag-AA) (b) (c) and (d) (armchair-AA) direction in the diffusive regime, as averaged over 1600 different disorder configurations. A lattice with sites and a polynomial degree cutoff of are used in these numerical calculations.

In a dirty metal the charge density is modified by disorder due to the random phase factors. When these random phase factors are simply averaged over disorder configurations, they gain an exponentially decaying factor. Due to the geometrical anisotropy in graphene, the influence of these random phase factors depends on the path direction. According to Feynman’s path integral representation, the coupling is dominated by the shortest path between two magnetic impurities. There always exists an even number of shortest paths along the armchair direction (Fig. 3b), but only one along the zigzag direction (Fig. 3a). Therefore, for the zigzag direction, the electron wave is scattered by the same disorder twice with the same phase (e.g. in Fig. 3a) when it returns to the origin. However, along the armchair direction, there are closed paths which include different impurities and thereby different scattering phases (e.g. and in Fig. 3b). If we average the modulation of the charge density over the phases and we find

(12)

and

(13)

where is the average over disorder configurations and is the modulation of the charge density in the clean system [Eq. (6)]. The modulation along the armchair direction [Eq. (13)] is dominated by the term coming from the closed loop connected by two unrelated paths, so that the average value is exponentially closer to the clean value, while the zigzag direction is dominated by [Eq. (12)]. In this sense, one may expect the RKKY interaction in the armchair direction to be less sensitive to the nonmagnetic weak disorder than in other directions. This may explain our previous results Lee et al. (2012) and the results shown in Fig. 2, where the averaged value in the armchair direction remains susbstantially larger than in the other directions and does not change its sign.

(a) zigzag     
(b) armchair      
Figure 3: (Color online) Schematic diagrams of the shortest path from the origin to the first lattice point along (a) the zigzag-AA direction and (b) the armchair-AA direction. denotes the phase shift due to nonmagnetic disorder and denote the modulation of the charge density amplitude.

In order to investigate in more detail how much nonmagnetic disorder affects the RKKY interaction and how the latter depends on direction in the graphene lattice and on particular lattice points, we have calculated the distribution of the interaction amplitude for a disorder strength and at lattice points which are on the zigzag or armchair direction. The results are shown in Fig. 4 a-e. In these calculations, different disorder configurations are used with for a lattice of sites, while the number of Chebyshev polynomials is fixed. When the distance between the magnetic impurities is smaller than the mean free path, the RKKY interaction multiplied by the cube of the distance has a shape which does not depend on the distance. Interestingly, the distribution has three different shapes along the zigzag-AA direction (Fig. 4a-c), which are repeated periodically every third site in that direction. However, there is only one type of distribution along the armchair direction, as can be seen in  (Fig. 4d). The oscillating factor in the semiclassical expressions for the charge density, in Eq. (6), takes only the values , along sites in the zigzag-AA direction. The RKKY interaction on sites which give the same numerical factor, either or , are found to have the same distribution. In order to directly compare them with each other, we plot together the distributions of the RKKY interactions obtained for these sites in Fig. 4e. The RKKY interaction at sites where the oscillating factor takes the value (-1/2) along the zigzag-AA direction has the broadest (green) distribution, while the RKKY interaction at site where the oscillating factor takes the value (1) has the narrowest (red) distribution.

The unaveraged RKKY interaction, i.e., that obtained for a particular disorder representation, is found to remain long ranged as in the clean system. However, we find that the oscillations acquire a random phase that also changes with the distance between magnetic impurities. Therefore, the average value does not characterize the RKKY interaction at distances exceeding the mean free path. To find the typical value, we evaluated the geometrical average of the interaction amplitude, which is defined by

(a)           zigzag1            
(b)           zigzag2            
(c)           zigzag3            
(d)          armchair            
(e)                                 
(f)                               
Figure 4: (Color online) Plots of the distribution of the RKKY interaction amplitude () multiplied by the cube of distance, , for the (a) first (b) second (c) third of the triplet of points along zigzag-AA direction (see Fig. 4 f for the definition of these points), and (d) for the armchair direction. The disorder strength is fixed to . For comparison, these distributions are plotted together in (e). A lattice with sites and a polynomial degree cutoff of are used. The lattice constant is set to unity.
(14)
(a)                                                     
(b)                                                     
Figure 5: (Color online) (a) The averaged density of states in a double logarithmic plot (b) and the geometrical averaged RKKY interactions along the zigzag-AA direction. For the density of state calculation, a lattice with sites and a polynomial degree cutoff of are used. For the RKKY interaction calculation, the same lattice size and polynomial cutoff of Fig. 2 are used. The black lines in (a) represent fittings to the relation with for , respectively.

The typical value is found to have the same power-law decaying behavior with distance as the amplitude of the RKKY interaction in the clean limit.Zyuzin and Spivak (1986); Bergmann (1987); Lerner (1993); Sobota et al. (2007); Bulaevskii and Panyukov (1986) Before a direct evaluation of the geometrical average, we have calculated the density of state for two disorder strengths   to observe how weak disorder affects the pseudogap at the neutrality point (), since, as we have seen above, the power of the pseudogap is directly related to the unconventional power-law decay of the RKKY interaction. For the density-of-states calculation, we used the KPM Weiße et al. (2006); Amini et al. (2009) method with sites and a polynomial degree cutoff of . As one may see from the plot of the density of states in  (Fig. 5a), the pseudogap is still not filled for weak disorder, has the same power law as in the clean limit, and only the slope around the neutrality point is changed,Amini et al. (2009)

(15)

where is the density of states and is the energy measured in units of the hopping amplitude . The slope depends on the disorder strength and is obtained by fitting the data in Fig. 6.

Using the Born and T-matrix approximations, an analytical study has reported that there is a logarithmic correction to the density of states around the neutrality point, yielding in the presence of disorder.Hu et al. (2008) This logarithm correction does not change the power of the distance dependence of the RKKY interaction [Eq.(6)]. Consequently, the geometrical average of the RKKY interaction is expected to have the same power-law exponent as the clean system (). The direct numerical calculation shown in Fig. 5b strongly supports this conclusion.

(a)                                                     
(b)                                                     
Figure 6: (Color online) Plots of (a) the absolute value of the RKKY interaction in a disordered sample as a function of distance and (b) of the distribution of the logarithm of the RKKY interaction amplitude for and different disorder strengths. The inset of (b) shows the distribution of the amplitude itself for . A lattice with sites and a polynomial degree cutoff of are used. The black dashed lines represent the resulting distribution functions with and for , respectively (see text).

Figure 6a shows how the RKKY interaction in a strongly disordered sample gets suppressed by several orders of magnitude when the strength of the nonmagnetic random potential is increased. In order to investigate the broadening of the distribution, we employed realizations of the disorder potential and calculated the RKKY interaction amplitude for a fixed distance . The results are shown in Fig. 6b. The inset indicates that the interaction strength follows a distribution with very long assymetric tails. The squared amplitude () has a distribution similar to log-normal with a width that increases with disorder strength. Using a field-theoretical approach valid in the metallic regime, Lerner found Lerner (1993) that the increase in the strength of the nonmagnetic disorder leads to a crossover in the shape of the distribution function from a broad non-Gaussian with the very long tails in the weak disorder regime to a completely log-normal distribution in the region of strong disorder regime. This may explain the crossover of the distribution which we observe from weak disorder to strong disorder  (). In order to directly compare we have fitted the results with a log-normal functional form, which is given by

(16)

where is the number of realizations and . The data are shown together with the fitting curves in Fig. 6b (black dashed lines).

The width  of the distribution in Eq. (16) has been analyzed as functions of the disorder strength and is shown in Fig. 7. The red line in Fig. 7 is a linear fitting curve () and agrees with the analytical prediction by Lerner, who has used the renormalization group method to obtain Lerner (1993)

(17)

where is the mean free path of electrons, which we studied in our previous work.Lee et al. (2012)

Figure 7: (Color online) Plot of the width of the distribution of RKKY interaction amplitudes as a function of the disorder strength  (in unit of ). The red line represents a linear fitting curve  which yields .

V Doped graphene ()

v.1 Clean system ()

By controlling in the Hamiltonian of Eq. (1), we investigate how the RKKY interaction evolves with the Fermi level. Results are shown in Fig. 8, where the interaction amplitude is multiplied by in order to emphasize its oscillatory behavior. In these calculations, we used a lattice with sites and a polynomial degree cutoff . Near the neutrality point a beating pattern appears as shown in Fig. 8a, b. It consists of a superposition of waves with the wave vectors and , where is the Fermi wave vector originating from the Friedel oscillations at finite Fermi energy. Recently, the following analytical expressions for the beating pattern were derived using lattice Green’s functions:Sherafati and Satpathy (2011a)

(18)

and

(19)

where is the RKKY coupling function at the neutrality point [see Eq. (8)] and is the Meijer-G function. Note that term within brackets describes the isotropic dependence of the oscillations on the Fermi momentum . The external prefactor, , on the contrary, is strongly anisotropic, depending on the vector given by the momentum difference between the two neighbored Dirac points . To make a comparison with our calculations, the function represented by Eq. (19) is also presented in Fig. 8 by a red dashed line. Excellent agreement is found. One can estimate the wavelength of the long oscillation, which appears at finite Fermi level, using the dispersion relation near the neutrality point, which is given by

(20)

where is the momentum relative to the Dirac point and is the Fermi velocity.Castro Neto et al. (2009) The Fermi wavelength is found to be about for . This coincides with the period seen in  (Fig. 8a). As expected from Eq. (19), the oscillations with large period seen in (Fig. 8b) are isotropic.

(a)                                              
(b)                            
(c)                            
Figure 8: (Color online) (a) The RKKY interaction amplitude multiplied by in doped graphene at the chemical potential along the zigzag-AB direction. Density plots of the RKKY interaction amplitude for (b) and (c) for all directions. In (a), the result of culations used the kernel polynomial method and the lattice Green’s function method are represented as solid blue and dashed red line, respectively. A lattice with sites and a polynomial degree cutoff of were used. The lattice constant is set to unity.

We have also calculated the RKKY interaction amplitude for highly doped graphene (). The results are multiplied by the square of the distance, , and are shown in the density plots of  (Fig. 8c). The behavior cannot be described by Eqs. (18) and (19) which are valid only close to the neutrality point. When the Fermi level is exactly at the van Hove singularity at  (), the ordering pattern of the RKKY interactions along the zigzag direction is reversed. In other words, the correlation between impurities on zigzag-AA or BB is always antiferromagnetic and ferromagnetic for zigzag-AB or BA pairs. At the same time, the interactions are strongly suppressed for the other directions. This is in accordance with a result of a previous study.Uchoa et al. (2011)

v.2 Disordered system ()

We have also calculated the RKKY interaction amplitude in disordered doped graphene and the results are shown in Fig. 9. A lattice with sites and a polynomial degree cutoff of were used in these numerical calculations. In our previous work,Lee et al. (2012) we reported that the amplitude of the averaged RKKY interaction along the zigzag-AA direction in the ballistic regime () increases with weak disorder (Fig. 4f). Even though the unusual oscillation coming from the interference of two Dirac points still exists (Fig. 9a), the amplitude of the envelope of the averaged RKKY interaction decreases with disorder strength  and the period () of the oscillation coming from the finite Fermi energy is modified by the disorder as in a conventional metal.

(a)                                 
(b)                                 
Figure 9: (Color online) The average RKKY interaction strength multiplied by the square of the distance, , along the (a) zigzag-AA and (b) armchair-AA directions at the Fermi energy is plotted in the diffusive regime. 1600 different disorder configurations are used in the averaging. The same lattice size and polynomial cutoff as in Fig. 2 were used.

Vi Discussion and Conclusions

In conclusion, we have employed a semiclassical description of the RKKY interaction in terms of the modulation of the charge density by the presence of magnetic impurities in order to rederive an analytical expression applicable to pure graphene. This semiclassical approach provides not only a simpler derivation, but also an intuitive interpretation of the origin of the oscillating behavior of the RKKY interaction at zero doping as an interference of the two degenerate Dirac points. Moreover, the origin of the unusual power-law decay of the RKKY interaction in pure graphene at the neutrality point is clearly related to the pseudogap in the density of states. Using a Feynman’s path integral scheme, we could also trace the origin of the anisotropic sensitivity of the RKKY interaction to disorder to the presence of multiple shortest paths between two magnetic impurities in the armchair direction, showing that this direction is more robust to disorder than the zigzag one.

As an extension of our previous study,Lee et al. (2012) we have calculated and studied the RKKY interaction in doped and disordered graphene in detail using the kernel polynomial method. Fnite gate voltage breaks particle-hole symmetry and the resulting finite Fermi surface yields Friedel oscillations, so that the sign of the RKKY interaction between the impurities localized on the same sublattice now oscillates with distance. When the Fermi level is exactly at the van Hove singularity (), the ordering pattern of the RKKY interactions along the zigzag direction is reversed. At the same time, the interactions are strongly suppressed for all other directions.

In order to study the anisotropic influence of nonmagnetic disorder, we evaluated the RKKY interactions along two different directions between the zigzag and armchair directions. As reported in our previous study and expected from the semiclassical approach, in the diffusive regime the armchair direction is not affected by the nonmagnetic disorder as much as the other directions. We have also found that, in the ballistic regime (), the distribution of the RKKY interactions along the zigzag direction is not universal but depends on the lattice sites at which the pair of magnetic impurities sit. We identified three different representative shapes, which repeat themselves periodically. By an accurate evaluation of the density of states around the neutrality point in weakly disordered regime (), we confirmed that the linear dispersion relation is still valid and the pseudogap is not filled. This is in full agreement with the fact that the geometrical average of the RKKY interaction in the diffusive regime decays as in the clean system, namely, as , and not as , which is the usual behavior for two-dimensional metals. In the localized regime (), the geometrical average is also exponentially suppressed at distances exceeding the localization length and the distribution of the strength of the RKKY interaction shows a crossover from the non-Gaussian shape with very long tails to the completely log-normal form when increasing the disorder strength. We have analyzed the width of the log-normal distribution and confirmed that it increases linearly with the amplitude of the disorder potential .

In this work we showed that the KPM method is efficient and accurate for studying interactions in disordered systems. In order to minimize the computation time while keeping the highest accuracy, the convergence of the calculations with respect to the Chebyshev polynomial cutoff degree and the system size have been investigated in detail (see Appendix A). The proper cutoff to reach convergence is found to increase linearly with the distance between two magnetic impurities, as seen in Fig. 10 b. For a given distance , a system size about five times larger than is found to give convergent results. In comparison with the exact numerical diagonalization, the KPM is found to be much faster. It can be implemented for very large system sizes, even in disordered ones, where thousands of realizations are needed in order to yield meaningfull statistics.

(a)                                                         
(b)                                             
Figure 10: (Color online) Plots of (a) the RKKY interaction as function of the polynomial degree in a lattice with sites. (b) The smallest cutoff number that yields a convergent result as function of the distance between the magnetic impurities .

Appendix A Convergence of the Kernel Polynomial Method Calculations

When using the KPM, the calculation of the Chebyshev polynomials using recurrence relations consumes most of the computation time. Therefore, we investigated first the relation between the cutoff number  and the convergence of the results in order to be able to minimize and optimize the calculations. For clean graphene, a lattice with sites was used in these calculations. When a cutoff number is not sufficient, the amplitude of the RKKY interaction deviates from the expected power-law behavior, as indicated by the blue and red lines in Fig. 10a). When we determined the smallest cutoff number such that the variance of the amplitude of the RKKY interaction is less than 5%, we found that it increases linearly with the distance between two magnetic impurities as seen in Fig. 10b). This linear relation between the distance and the cutoff number allows the rapid calculation of the RKKY interaction amplitude between the magnetic moments even at large distances .

In order to minimize the KPM calculation time further, we have also studied the smallest system size which yields convergent results, as shown in Fig. 11. denotes the longest distance used in the calculation, and the cutoff number is . One can observe a good convergent behavior (see green and black line in Fig. 11) when the system size is larger than . The exact diagonalization method also yields a proper result when .Black-Schaffer (2010) However, the KPM does not require matrix diagonalization and therefore is a much faster computational tool.

(a)
(b)
Figure 11: (Color online) Plots of the RKKY interaction in terms of a system size . is used as the longest distance . A cutoff number is used in these calculations.
Acknowledgements.
This research was supported by WCU (World Class University) program through the National Research Foundation of Korea funded by the Ministry of Education, Science and Technology(Grant No. R31-2008-000-10059-0), Division of Advanced Materials Science. E.R.M. acknowledges partial support through the NSF DMR (Grant No. 1006230). E.R.M. and G.B. thank the WCU AMS for its hospitality. H. Y. Lee thanks Jacobs University Bremen for its hospitality.

References

  • Chen et al. (2007) J.-H. Chen, L. Li, W. G. Cullen, E. D. Williams, and M. S. Fuhrer, Nat. Phys. 7, 1745 (2007).
  • Meier et al. (2008) F. Meier, L. Zhou, J. Wiebe, and R. Wiesendanger, Science 320, 82 (2008).
  • Saremi (2007) S. Saremi, Phys. Rev. B 76, 184430 (2007).
  • Black-Schaffer (2010) A. M. Black-Schaffer, Phys. Rev. B 81, 205416 (2010).
  • Sherafati and Satpathy (2011a) M. Sherafati and S. Satpathy, Phys. Rev. B 83, 165425 (2011a).
  • Sherafati and Satpathy (2011b) M. Sherafati and S. Satpathy, Phys. Rev. B 84, 125416 (2011b).
  • de Chatel (1981) P. F. de Chatel, J. Magn. Mater. 23, 28 (1981).
  • Jagannathan et al. (1988) A. Jagannathan, E. Abrahams, and M. J. Stephen, Phys. Rev. B 37, 436 (1988).
  • Zyuzin and Spivak (1986) A. I. Zyuzin and B. Z. Spivak, JEPT Lett. 243, 43 (1986).
  • Lerner (1993) I. V. Lerner, Phys. Rev. B 48, 9462 (1993).
  • Bulaevskii and Panyukov (1986) L. N. Bulaevskii and S. V. Panyukov, JETP Lett. 43, 240 (1986).
  • Lee et al. (2012) H. Lee, J. Kim, E. R. Mucciolo, G. Bouzerar, and S. Kettemann, Phys. Rev. B 85, 075420 (2012).
  • Sobota et al. (2007) J. A. Sobota, D. Tanasković, and V. Dobrosavljević, Phys. Rev. B 76, 245106 (2007).
  • Bergmann (1987) G. Bergmann, Phys. Rev. B 36, 2469 (1987).
  • Cirone et al. (2002) M. A. Cirone, J. P. Dahl, M. Fedorov, D. Greenberger, and W. P. Schleich, Journal of Physics B: Atomic, Molecular and Optical Physics 35, 191 (2002).
  • Kogan (2011) E. Kogan, Phys. Rev. B 84, 115119 (2011).
  • Roche and Mayou (1999) S. Roche and D. Mayou, Phys. Rev. B 60, 322 (1999).
  • Weiße et al. (2006) A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, Rev. Mod. Phys. 78, 275 (2006).
  • Amini et al. (2009) M. Amini, S. A. Jafari, and F. Shahbazi, EPL (Europhysics Letters) 87, 37002 (2009).
  • Hu et al. (2008) B. Y.-K. Hu, E. H. Hwang, and S. Das Sarma, Phys. Rev. B 78, 165411 (2008).
  • Castro Neto et al. (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  • Uchoa et al. (2011) B. Uchoa, T. G. Rappoport, and A. H. Castro Neto, Phys. Rev. Lett. 106, 016801 (2011).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
145105
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description