Scattering of Dirac Electrons by Randomly Distributed Nitrogen Substitutional Impurities in Graphene

The propagation of wave packets in a monolayer graphene containing a random distribution of dopant atoms has been explored. The time-dependent, two-dimensional Weyl-Dirac equation was solved numerically to propagate an initial Gaussian-type wave front and to investigate how the set of impurities influences its motion. It has been observed that the charge transport in doped graphene differs from the pristine case. In particular, nitrogen substitutional doping reduces the charge mobility in graphene due to backscattering effects.


Introduction
Graphene-a single layer of carbon atoms arranged in a honeycomb pattern-is the basic material for the family of low-dimensional allotropes of carbon known as carbon nanomaterials, such as carbon nanotubes and fullerenes, which exhibit unusual physical properties [1].The high electron mobility, among other characteristics, makes graphene a subject of great interest in contemporary physics [2][3][4][5].The fact that graphene is a zero-gap semiconductor may be a disadvantage for some applications.A solution is to chemically dope the structure so as to shift the Fermi level or to open a sizable gap, or both [6].In particular, simple substitutional doping with B [7], N [8], and other elements can be realized by a variety of methods, either during the growth or after the graphene synthesis [9].In this respect, low-energy implantation of the ionized dopant atoms is a promising technique [10].
Nitrogen is an important element for producing n-doped nanocarbons.For that reason, doping carbon-based nanostructures with N has been investigated since a long time [11].In recent years, different techniques have been developed to incorporate nitrogen in graphene at concentration of a few atomic percent [12][13][14].Along with substitutional nitrogen impurities, pyridinic and pyrrolic defects are created in variable proportions.In heavily-doped graphene compounds, the latter defects confer new and remarkable properties [14] to the samples and play an important role on their transport properties [12].In particular, field-effect transistors with strong n-type characteristics, which remain stable in air for a long period, have been realized from post-synthesis N-doped graphene samples [13].
It is of course important to understand the effects of chemical doping on the electronic structure of graphene and on its quantum transport properties [15].The electronic properties of graphene doped with nitrogen atoms have been analyzed by several groups [16][17][18][19][20][21][22][23].Substitutional N dopants shift the Fermi level above the Dirac energy, modify the electronic density of states and open a small band gap [16,17,19].Also, theory predicts that charge transport in graphene becomes diffusive in the presence of N or B impurities at a concentration above ∼0.5 at.% [24].All these effects affect the conductivity of the samples in a non-trivial way.There are more carriers available for transport in these systems than in pure graphene but, on the other hand, scattering by the dopant atoms reduces the electron relaxation time [23] and, hence, the mobility.
All the theoretical works referred to above were performed in the framework of the Schrödinger equation.The two-dimensional Weyl-Dirac equation (Dirac equation for massless fermions) is a continuous-medium approximation of the Schrödinger equation of graphene that is valid at low energies, up to 1-2 eV excitations [39,40].In this framework, the dynamics of wave packets has been investigated recently in an ideal monolayer graphene [41,42] and, subsequently, in a random potential landscape that models charge puddles [43,44].The scattering of Schrödinger and Weyl fermions by potential and magnetic barriers has also been the subject of recent studies [45,46].For solving the time-dependent Schrödinger and Weyl-Dirac equations, a numerical method based on the split-operator technique has been used [47,48].
In the present work, we investigate the propagation of a wave packet in monolayer graphene containing a random distribution of nitrogen dopants for various concentrations up to 0.5 at.%.Within the continuum model, the dopants were modeled by attractive potential wells around randomly localized points in the graphene plane, and the Weyl-Dirac equation of motion for the wave packet was numerically solved.These calculations allowed us to investigate how the set of impurities influences the electron transport properties of the doped graphene sample.The transmission of the wave packet across the perturbed region was investigated for various concentrations of dopants.The obtained results suggest that the conductance at a fixed Fermi energy slightly decreases with increasing concentration of dopants.Being a continuous-medium approximation, the model ignores fine details linked to the interactions between impurities, responsible for sublattice segregation observed by STM [49] and predicted theoretically to take place below about 0.2 at.% concentration of N [21].Also, the model disregards any inter-valley scattering mechanism, because the two Dirac cones of graphene are decoupled.Therefore, the present calculations evaluate the sole contribution of intra-valley mechanism, which dominates the elastic collisions of electrons and holes in graphene on weak-scatter defects.
Although less perturbing than its pyridine-like and other forms [50,51], substitutional nitrogen is not exactly a weak scatterer.This dopant was chosen because of its intrinsic technological importance for some applications of graphene [8,10].The results are not expected to be quantitative, because of the approximation made.Nevertheless, the formalism developed here is not without interest, for it may serve as a "numerical laboratory" to check sophisticated theoretical formulation of the Drude conductivity of Dirac-Weyl fermions [52,53].In fact, the continuous-medium model worked out in this paper is for graphene what the jellium model is for normal metals.It is not necessary to emphasize how much important the physics of jellium is for a qualitative understanding of the metallic state.Let is suffice to mention that the calculations presented below correspond to a system that would represent 40 millions atoms if it were described at the atomic level.

Theoretical Model
We consider an infinite graphene sample, by setting periodic boundary conditions in both x and y directions, with a period 128 nm in the former and 8192 nm in the latter.The two ends of the ribbon in the y direction are then considered as leads, as a Gaussian shaped pulse of charge carriers is launched from one lead (at y = y 0 = −80 nm) and propagates through the system towards y → +∞.Around a given region of the graphene plane, scattering potential wells mimic the effects of substitutional N dopants, while the dopants concentration is related to the number of scattering centers per area unit.The perturbed region covers the full width of the ribbon and has a spatial extension of 40 nm in the y direction.As initial wave packet, we consider a Gaussian wave front as sketched in Figure 1, namely, a wave packet that is invariant along the x-direction, but which has a width a y in the y-direction, multiplied by a pseudo-spinor (1 i) T and by a plane wave where N is a normalization constant and k is the average wave vector, which, in graphene, is related to the wave-packet energy E by k = E hv F , where v F is the Fermi velocity.Since this pseudo-spinor is such that σ y = 1 and σ x = σ z = 0, where σ is the Pauli matrix vector, this choice of pseudo-spinor is convenient: according to Heisenberg equation, the velocity is given by where and the wave packet is therefore expected to propagate in the direction of its pseudo-spinor (in this case, the y-direction), with constant velocity v F , no matter how much energy it has [44,45].
| (x,y,0)|² We chose to work with a wave front, instead of a circularly symmetric wave packet, in order to avoid zitterbewegung (a trembling motion of the wave packet) [47,54] along the x-direction.In the latter case, one would not have been able to define an exact value of k x and, consequently, σ -and, ultimately, the velocity-would not have been a constant of motion.All the results in the remainder of this paper were obtained for a wave packet of width a y = 20 nm and energy E = 100 meV.In the following, the wave-packet energy will remains pinned at this 100 meV value for all N concentrations.Our primary objective is to investigate how the concentration of scattering centers affects the transport properties while keeping constant the Fermi energy.Kubo-Greenwood calculations indicate that the conductance of graphene containing 4 at.%random N impurities does not vary with energy above the Dirac point [26].This result makes it non-essential to adapt the Fermi energy to the doping level.In fact, the Fermi energy of the system depicted in Figure 1 is imposed by the graphene leads.The chosen value of 100 meV corresponds to graphene doped at a charge density n e = 10 12 cm −2 , which could be the result of defects or be obtained by other means such as an external gate potential.
The propagation of the wave packet is obtained by applying the time-evolution operator where the Hamiltonian H, assumed to be time-independent, is the one for low energy electrons in graphene in the presence of a potential landscape V(x, y), [39,40] with I as the 2 × 2 identity matrix and wave functions written as pseudo-spinors is the probability of finding the electron in the sub-lattice A (B) of graphene.We separate the potential and kinetic energy terms of the time-evolution operator through the split-operator technique [55,56], where terms of order higher than O(∆t 3 ) are neglected as an approximation, by taking a small time step ∆t.The advantage of this factorization is that it allows us to perform multiplications in real and reciprocal spaces separately, thereby avoiding writing the momentum operator as a derivative, just by making a Fourier transform of the functions and writing p = h k.Furthermore, the exponential of Pauli matrices can be re-written exactly as matrices, which simplifies the calculations even more [46,48].We performed wave-packet propagation with a time step as small as ∆t = 0.1 fs and kept track of the probabilities of finding the electron before (P 1 ), inside (P 2 ), and after (P 3 ) the scattering region.
The large-t limit of the probability P 3 was then identified as the transmission probability through the perturbed region.
Notice that performing a Fourier transform while solving Equation ( 5) is similar to applying periodic boundary conditions to the system.Thus, in order to avoid spurious reduction (enhancement) of the transmission (reflection) probability P 3 (P 1 ) when the wave packet reaches the upper edge of the numerical grid, we considered a very long grid in the propagation direction.Specifically, the graphene ribbon contained 1024 × 65,536 grid points regularly spaced by 0.125 nm in both x and y directions.Fourier components of the wave packet were computed on discrete wave-vector nodes regularly distributed in reciprocal space with 0.049 nm −1 and 0.00077 nm −1 step sizes in the transverse and longitudinal directions, respectively.
Nitrogen doping of graphene usually leads to randomly distributed substitutional sites [49,[57][58][59].Accordingly, doping was simulated by a random distribution of scattering centers within an area of A = 128 nm × 40 nm around the y = 0 axis of the ribbon, as sketched in Figure 1.The N atoms, which transfer 0.35-0.45electron to the graphene host [10,60], are positively ionized.As revealed by ab-initio calculations [30,61], the perturbation produced on the π electrons extends somewhat around each impurity, due to Coulomb attraction.Each dopant atom is thus represented by a potential well of depth |U| and width σ [61]: where ρ is the distance to the atom.The parameters to be used for well-separated nitrogen atoms are b = 0.15 nm and U = −4.0eV [19,61].In order to avoid any overlapping of the potential wells generated by the dopants, the selection of the scattering centers was constrained by the requirement that the distance between two nitrogen centers remains greater than 2 nm (for dopant concentrations 0.1 at.%-0.4 at.%) or 1.8 nm (for dopant concentration 0.5 at.%).

Results and Discussion
As previously mentioned, the choice of the pseudo-spinor components of the initial wave packet (Equation ( 1)) is such that in pure graphene, the wave packet conserves its spatial shape and travels at the constant velocity v F set to 8.7 × 10 5 m/s = 0.87 nm/fs.This soliton-like propagation is clearly seen by the P 2 probability displayed in Figure 2a in dashed line as a function of time, together with the probabilities of finding the electron in region 1 (P 1 ) and 3 (P 3 ).These curves are shown for comparison with the doped case.After 96 fs, the center of the wave packet has reached the coordinate y = 0 at the middle of region 2, P 2 reaches its maximum value.The wave packet has almost completely left the region 2 after 175 fs, when the probability P 3 saturates to one.Figures 2b, 3 and 4 represent the probabilities of finding the electron before (P 1 ), inside (P 2 ), and after (P 3 ) the scattering region, as a function of time, for wave packets propagating through the nitrogen-doped region, with concentrations 0.1 at.% and 0.4 at.%, respectively.The case 0.1 at.% is also illustrated in Figure 2a on a smaller time interval.The effects of doping show up after 100 fs.P 1 drops to a very small value after 120 fs and then starts increasing slowly, as a consequence of reflection of the wave packet by the scatterers.Just like in pristine graphene, P 2 reaches its maximum after 100 fs, decreases rapidly immediately after, but then much slower after 150 fs.This is because part of the wave packet remains trapped in the scattering region.Due to the finite extension of the perturbation region, the trapped fermions eventually escape from region 2 on both sides, giving rise to a non-zero reflection probability and to a transmission probability smaller than one.The wave function escapes as small packets, giving rise to amplitude oscillations on the probability density functions P i (i = 1, 2 and 3).However, since these oscillations come from the wave packet being bounced back and forth by the scatterers in the doped region, they are configuration dependent, so that averaging between several impurity distributions for the same dopant concentration leads to smooth curves, as illustrated in Figure 4.This figure reveals that the transmission probability fluctuates by less than 2% from one configuration to the other.The observed small configuration dependence of the results is partly due to the wave packet extending uniformly across the width of the graphene strip (see Figure 1), which produces an averaging process over the lateral positions of the dopant atoms.Thus, due to a prohibitively long time of computation, just one configuration was generated for the other dopant concentrations.It can be clearly observed that P 3 saturates faster for the smallest values of N concentration, which suggests that the dwell time of the wave packet inside the scattering region becomes longer as the concentration of nitrogen increases.The dwell time, whose definition follows from Equation (7), is listed in Table 1 for different dopant concentrations.The data indicate that the dwell time increases by about a factor of two when the N concentration ranges from 0.1 to 0.5 at.%.The overall increase of time with increasing impurity concentration is the result of the wave packet bouncing back and forth several times inside the perturbed region.By comparison, the transmission time of the wave packet in pure graphene is 150 fs, here defined as the time interval between the moment P 2 starts to increase and the moment when P 2 falls back to zero (Figure 2a).In the presence of scattering centers, the probability of finding the electron after the perturbed region decreases slightly, but significantly, as the concentration of dopant increases, as shown by the blue curve in Figures 2b, 3 and 4. The asymptotic value of P 3 was extrapolated by an exponential law assumed to govern its large-t behavior: The parameters T (transmission coefficient), α and γ were fitted to the numerical calculations of P 3 (t) for t > 200 fs.The obtained transmission coefficient T together with the parameter γ −1 , which is considered as an indicator of the dwell time, are listed in Table 1 vs. N concentration c.From Landauer theory of quantum transport [62], T for a given dopant concentration is a measure of the electrical conductance G(c) of the 40 nm doped area at the energy of 100 meV.As mentioned above, the conductivity of N-doped graphene varies weakly with energy above the Dirac point [26].This means that the arbitrary choice of 100 meV does not entail the trends of the results.The trends are clearly that the relative conductance G(c)/G(0), estimated by the ration T(c)/T(0), decreases with increasing concentration with an average slope 0.19 when c is in at.%.This behavior is in agreement with the results of DFT calculations of graphene conductivity performed for higher N concentrations with the Kubo formula [23].These calculations, in addition, pinpoint the dominant role played by the π electrons in the decrease of conductance with increasing doping.This finding validates the use of π-electron Hamiltonian considered here.To finish this discussion, it is worth mentioning that a significant reduction of electrical conductance of graphene upon B or N doping has been reported experimentally [10].
The observed decrease of conductance G(c) reflects the decrease of mobility µ = G(c)/n e e with increasing doping [63].In principle, the mobility should decrease with increasing the length of the doped region.However, Landauer-Büttiker calculations reported in ref. [26] show that the conductivity of N-doped graphene does not depend significantly on the channel length for electron energies above the Dirac point, as considered here.For that reason, it is expected that the data of Table 1 are not much sensitive to the 40 nm width of the perturbed region.
Further comparisons of the data of Table 1 with the existing literature is not easy for several reasons.Landauer-Büttiker DFT calculations for doped graphene nanoribbons depend strongly on the edges (zigzag or armchair) and on the width of the ribbon [33][34][35].Furthermore, the conductance is spin dependent in the case of narrow zigzag nanoribbons [34,37], because of magnetic ordering of the edge spins.The Kubo-Greenwood formalism is better adapted to two-dimensional graphene then the Landauer-Büttiker approach is, but there is no systematic calculations of the conductance versus dopant concentration, except in Ref. [23] as already mentioned above.
In the present work, there is no free edges due to the application of periodic boundary conditions.The system under investigation behaves similar to an infinite graphene sheet due to the large width of the graphene strip used.This argument is further supported by the experimental demonstration that transverse quantization in graphene nanoribbons has negligible effects on its room temperature conductivity, already when it is 30 nm wide [64].
Lastly, it is worth emphasizing that the 100-meV energy of the wave packet is the center of a Gaussian distribution whose standard deviation δE = hv F /a y is related to the spatial width a y of the wave packet.The wave packet therefore probes the transport properties of the doped region in an energy interval of half width δE.With the used value a y = 20 nm, δE = 29 meV, namely room temperature.In other words, smearing effects of the Fermi-Dirac distribution at room temperature was implicitly included in the calculations.

Conclusions
In summary, we have investigated the transmission of a wave packet through nitrogen dopant atoms randomly distributed in a graphene strip by solving the time-dependent Weyl-Dirac equation with periodic boundary conditions.Our results demonstrate that the dwell time of the wave packet inside the scattering region is significantly increased by the presence of the nitrogens.The observed backscattering probability of charge carriers by the dopant atoms is responsible for a decrease of the transmission of the wave packet and a concomitant reduction of the electron mobility.

Figure 1 .
Figure 1.(Color online) Sketch of the graphene ribbon, illustrating the geometry of the problem addressed in this work.The colors represent the potential energy V(x, y) (in eV) generated by a given configuration of N atoms located in the perturbed region.The graphene parts below and above the doped region serve as leads.The shaded area represents the initial Gaussian wave packet (Equation (1)) in the bottom lead and propagating upwards.

Figure 2 .
Figure 2. (Color online) (a) Probabilities of finding the electron before (P 1 ), inside (P 2 ), and after (P 3 ) the scattering region, as a function of time, for a E = 100 meV wave packet propagating in pristine graphene (dashed line) and through a 40 nm wide doped region with 0.1 at.%N (full line); (b) Same as in (a) over a longer time interval for the 0.1 at.% doping concentration.The shortest distance between impurities is 2 nm.The number of scattering centers for this concentration is 195.

Figure 3 .
Figure 3. (Color online) Probabilities of finding the electron before (P 1 -black curve), inside (P 2 -red curve), and after (P 3 -blue curve) the scattering region, as a function of time, for a E = 100 meV wave packet propagating through the 40-nm wide nitrogen doped area with 0.4 at.% concentration.

Figure 4 .
Figure 4. (Color online) Effect of configurational averaging on the probabilities P 1 (black curve), P 2 (red curve), and P 3 (blue curve) for the N concentration 0.3 at.%.The dashed lines represent the probability evolution obtained in a small time interval for one configuration of dopants, the full lines represent the same curves averaged over nine independent configurations.

Table 1 .
Transmission coefficient T and dwell time γ −1 (Equation (7)) of a 40 nm wide strip of graphene doped with a given concentration c of nitrogen atoms.