First-principles multiscale modelling of charged adsorbates on doped graphene

Adsorbed atoms and molecules play an important role in controlling and tuning the functional properties of two-dimensional (2D) materials. Understanding and predicting this process from theory is challenging because of the need to capture the complex interplay between the local chemistry and the long-range screening response. To address this problem, we present a first-principles multiscale approach that combines linear-scaling density-functional theory, continuum screening theory and large-scale tight-binding simulations into a seamless parameter-free theory of adsorbates on 2D materials. We apply this method to investigate the electronic structure of doped graphene with a single calcium (Ca) adatom and find that the Ca atom acts as a Coulomb impurity which modifies the graphene local density of states (LDOS) within a distance of several nanometres in its vicinity. We also observe an important doping dependence of the graphene LDOS near the Ca atom, which gives insights into electronic screening in graphene. Our multiscale framework opens up the possibility of investigating complex mesoscale adsorbate configurations on 2D materials relevant to real devices.


Introduction
Adsorbate engineering is a promising route towards modifying the electronic structure of graphene and other 2D materials for emerging functional device technologies [1][2][3][4]. For example, adsorbed atoms and molecules can change the carrier mobility [5,6], effective dielectric constant [7], chemical potential [8] and local magnetic moment [9] of graphene.
Charge transfer from adsorbates to the graphene substrate results in doping and the creation of Coulomb impurities. From a theoretical perspective, the description of such impurities and their interaction with the graphene substrate constitutes a challenging multiscale problem: while the charge Z transferred between the adsorbate and substrate is principally determined by the local chemistry (i.e. the coupling between adsorbate and substrate orbitals), the screening response of graphene is long-ranged with an extent of several nanometres [10]. Experimentally, the structural properties and screening response of adsorbed species on graphene as well as other point defects can be investigated with scanning probe techniques, such as scanning tunneling microscopy and spectroscopy or atomic force microscopy [8,[10][11][12][13][14].
Previous theoretical studies of adsorbates on graphene have focussed either on the description of the short-range charge transfer or the long-range graphene response. First-principles density-functional theory (DFT) has been used to determine the charge transfer for different adatom species; however, the supercells in these calculations were much too small to capture the long-range screening response of graphene [15][16][17][18]. Conversely, theoretical calculations based on tightbinding (TB) or continuum Dirac models [10,[19][20][21], while able to access large length scales, require external parameters (most importantly, the impurity charge Z) and so do not possess the unbiased predictive power of first-principles methods. Furthermore, such models typically employ uncontrolled approximations for the description of electron-electron interactions, which are either completely neglected [10,12,19,22] or treated within linear-response [23,24] or Thomas-Fermi (TF) theory [25][26][27]. Thus, there is a need for a parameter-free theory capable of bridging the length scales relevant to graphene-adsorbate interactions.

First-principles multiscale modelling of charged adsorbates on doped graphene
In this study, we combine first-principles DFT, continuum screening theory and model TB Hamiltonians within a novel parameter-free multiscale approach in order to overcome the limitations of each method in isolation. We demonstrate the power of our approach by applying it to the case of an isolated Ca atom adsorbed on doped graphene. This is a case of particular interest due to ongoing work on the detailed experimental characterization of the doping-dependent screening response [28], and the observation of supercritical states for small clusters of Ca adatom dimers [12].
Ca donates electrons to the graphene resulting in an impurity charge of Z = 1.6 e (with e denoting the proton charge). We compute the graphene local density of states (LDOS) in the vicinity of the charged Ca atom and observe that it decays to the bulk graphene value on a length scale of several nanometres. The decay length depends on the chemical potential which influences graphene's screening response to the impurity potential.
The remainder of this article is structured as follows: in section 2 we describe our general multiscale approach that bridges large-scale first-principles electronic structure calculations, continuum models and TB simulations at the mesoscale, and details of how we apply this approach to the specific case of a Ca adsorb-ate on a graphene substrate; in section 3 we present and discuss our results; and we make some concluding remarks in section 4. Figure 1 shows the multiscale approach we developed to model the properties of adsorbates on doped graphene. The method consists of three key steps. First, we perform first-principles DFT calculations of the pristine graphene and the graphene with adsorbed species. From the former, we derive the pristine graphene lattice parameter a and the Fermi velocity v F , which we use in the continuum and TB Hamiltonians of the second and third steps. The periodic supercell for the graphene with adsorbate must be sufficiently large to converge chemical interactions between the adsorbate and the substrate [29]. To sample different values of the carrier density and, thereby, chemical potential μ, we perform a series of calculations with different numbers of electrons in the supercell.

Outline of multiscale approach
Secondly, we parametrize a continuum screening model of the graphene with the adsorbed species. In this approach, the adsorbate is treated as a point impurity charge Z located above the graphene sheet, which is described by the Dirac Hamiltonian (i.e. without atomistic detail). To determine the charge transfer between the adsorbate and graphene, we compute the selfconsistent (or screened) potential of the continuum model in the same supercell that was used for the DFT calculations and, for each value of μ, we fit the value of Z to achieve optimal agreement with the DFT selfconsistent potential. Different continuum screening models may be used in this procedure, corresponding to different treatments of electron-electron interactions. In this work we employ non-linear Thomas-Fermi theory [25] (NLTF), but also include interband transitions through an effective dielectric function. The screened potential at a position r r in the graphene plane resulting from an impurity at ( ) = z R R 0, 0, imp is given by Zv n n v r r r r r r r r r r r r ; d , z scr denote the electron density of graphene without and with the adsorbate, respectively, and v r r z ( ) is the Coulomb interaction screened by graphene interband being the interband contribution to the graphene dielectric function in the random-phase approximation [30][31][32]. In the above equation, the first term captures the bare impurity potential and the second term is the induced potential arising from electronic screening. We denote this model, which reduces to the well-known randomphase approximation in the linear-response limit [30][31][32], as NLTF + inter. We further included scr ( ) due to the adsorbate is calculated using a continuum screening model; the model depends on the charge transfer Z, which is obtained by fitting the model potential to the first-principles self-consistent potential extracted from DFT simulations of the graphene with adsorbate system. additional screening due to graphene σ bands [33] and exchange-correlation effects [20], but found that they have a small effect and can be neglected.
Thirdly, we perform a TB calculation in a supercell that is sufficiently large to capture the decay of the screened impurity potential. We find that the supercells needed to converge the screened potential are several orders of magnitude larger than the supercells needed to converge chemical interactions. The screened potential calculated from the continuum model, with the value of Z that has been fitted to DFT, is used as an on-site term in the TB Hamiltonian. The TB calculation then yields the local density of states (LDOS), which can be compared to scanning tunneling experiments [10,28].

Application to calcium adatoms on graphene
In this section, we describe the application of the multiscale approach to isolated Ca adatoms on doped graphene.

Structure
All calculations use hexagonal supercells containing × n n graphene unit cells (2n 2 C atoms) and a single Ca adatom. We use a graphene lattice constant a of 2.47 Å , obtained from plane-wave DFT calculations of a single unit cell with an energy cutoff of 1000 eV, a × 28 28 k-point grid (centred at the Γ-point), norm-conserving pseudopotentials, 20 Å separation between periodic images of the graphene sheets and the PBE exchangecorrelation energy functional [34]. These calculations were carried out using the castep [35] code. For the Ca adatom we use the geometry found by a previous DFT study of the system [15], with the adatom 2 Å above the hollow site at the centre of a hexagon. The same study also reported negligibly small relaxations of the graphene atoms near the Ca adatom and we thus keep these atoms fixed in our calculations.

First-principles calculations
We perform DFT calculations in a × 56 56 graphene supercell (6272 C atoms and a single Ca adatom). While this supercell is not large enough to capture the decay of the screened impurity potential, it is sufficiently large to describe chemical interactions resulting from the coupling of adsorbate and graphene states, such as exchange interactions giving rise to spin-split states. Such calculations are already very large by the standards of DFT, and are made possible by using the onetep [36,37] DFT code whose computational cost scales only linearly with the size of the system, in contrast to the cubic scaling of conventional methods. We use normconserving pseudopotentials (including the 3s and 3p semi-core states for the Ca atom) and the PBE exchangecorrelation functional. The Brillouin zone is sampled at the Γ-point. Kohn-Sham states are expanded in a basis of non-orthogonal, atom-centred orbitals with a radius of 5.3 Å that are variationally optimized in situ, and described on a real-space grid corresponding to a plane-wave energy cutoff of 1000 eV. As we are dealing with a semi-metallic system, we do not truncate the density matrix and employ an ensemble-DFT [38] formalism with an electronic temperature of 50 K. The simulation employs periodic boundary conditions in all three directions, and we use a periodic image separation perpendicular to the graphene plane of 20 Å .
In experiments on graphene in a field-effect transistor setup the density of carriers can be modified by application of a gate voltage. To explore the effect of varying the carrier density in our simulations, we perform separate DFT calculations with different numbers of electrons in the supercell. Specifically, by adding and subtracting pairs of electrons we vary the total system charge ∆Q from −4 e to +8 e, corresponding to carrier densities n between − × 3 10 12 cm −2 and × 5 10 12 cm −2 . For calculations with a non-zero total charge ∆Q we use a compensating homogeneous background charge.
The self-consistent (screened) local potential from each DFT calculation, which is comprised of the ionic, Hartree and exchange-correlation terms, is locally averaged in a Voronoi cell around each C atom in the graphene plane in order to smooth out the large variations that occur on the scale of an interatomic spacing. This smoothed DFT potential ∆ V Q r r; scr DFT ( ) is then used to parametrize the continuum model.

Continuum screening calculations
The screened impurity potential from the continuum screening model is calculated with a custom in-house code with 2D periodicity [39]. The centre of the impurity potential is placed 2 Å above a hollow site at the centre of a hexagon of the graphene plane, as in the DFT geometry. For the non-linear models, a self-consistency cycle is used with a mixing scheme in which the screened potential of the next iteration is composed of 1% of the screened potential of the current iteration and of 99% of the screened potential from the previous iteration.
The continuum screening model depends on three parameters: Z, μ and the Fermi velocity v F . The substrate dielectric constant usually constitutes an additional parameter; here, however, we are interested in suspended graphene and so set it to unity. = × v 0.8 10 F 6 m s −1 is taken directly from the pristine graphene DFT band structure, and Z and μ are treated instead as fitting parameters of the model.
The accuracy of the model is evaluated against DFT using a fitness metric ( ) µ F Z, , which is defined as the root-mean-square difference between the model potential V scr cont of equation (1) and the DFT potential V scr DFT , integrated over the supercell area. We calculate V scr cont for a range of Z and μ values using the same supercell as in the DFT calculation (including periodic boundary conditions). The values of Z and μ that minimize F are those that give the optimal fit of the screening model to DFT (see figure S1 of the supplementary material (stacks.iop.org/TDM/4/025070/mmedia)).
To study the effect of different approximations for the electron-electron interactions, we also use three other continuum screening models: (i) TF theory within a linear-response formalism (LTF), which only includes intraband transitions between the Dirac bands of graphene; (ii) linear-response TF theory plus interband transitions (LTF + inter); (iii) non-linear TF theory (NLTF) [25] without interband transitions.

Tight-binding calculations
The large-scale TB simulations are performed in a × 168 168 supercell (56448 C atoms and a single Ca adatom); this is an order of magnitude increase in system size with respect to the underlying DFT simulations. We use a custom in-house TB code for graphene with 2D periodicity. The TB Hamiltonian includes the C p z orbitals and nearest-neighbour hoppings only with a hopping energy of 2.54 eV. A × 2 2 Monkhorst-Pack mesh (centred at the Γ-point) is used to sample the Brillouin zone.

Charge transfer
The charge transfer Z is an important quantity for characterizing the interaction between an adsorbate and the substrate. Standard approaches for calculating Z partition the continuous electron density obtained from DFT into an adsorbate and a substrate contribution. For example, Chan et al [15] calculated the charge transfer for different adatoms on graphene, but found large discrepancies between two different methods of partitioning the electron density (0.18 e and 0.78 e for Ca adatoms). As shown in the inset of figure 2 and figure S2 in the supplementary material, we obtain an even broader range of values for Z, from 0.2 e to 1.4 e, depending on the particular charge partitioning method used. While each method provides a different value for Z, all of them find the charge transfer to be independent of the supercell charge ∆Q, indicating that Ca states are neither filled nor depleted as the graphene doping level is varied.
In contrast to these approaches, no partitioning of the charge density is required in our multiscale method. Instead, Z is obtained by imposing that the screened impurity potential of the continuum screening model reproduces the self-consistent DFT potential. In the fitting procedure, we enforce a μ-independent value of Z by performing a global fit to all doping levels at once with equal weighting. This is justified by the charge stability of the impurity discussed previously. It should be noted that performing individual fits for each doping level indeed results in a fairly constant Z for all charge states except ∆ = Q 2 e, in which case the value is enhanced by up to ∼50%; we believe this to be because the graphene system is very close to being undoped and, hence, is the least well-described by our models as the screening is dominated in this case by the interband transitions [40]. Table 1 gives the value of Z obtained from the new multiscale approach using different approx imations for the description of electron-electron interactions. The NLTF + inter method, which offers the most accurate description of electron-electron interactions, yields Z = 1.6 e. The same value is obtained from the NLTF method indicating that inclusion of interband trans ition does not change the value of Z. In contrast, the linear-response LTF and LTF + inter approaches yield significantly smaller values for Z. Table 1 also gives the value of the fitness function F averaged over all DFT doping levels. It is interesting to note that there is almost no difference in this value between the different approaches 1 ; all of them reproduce the DFT screened potential with high accuracy (<10 meV), but with significantly different Z values. This observation was explained by Katsnelson [25] who pointed out that the screened potential of a charged impurity in doped graphene has the same functional form in linear and non-linear Thomas-Fermi theory, but with a different prefactor corresponding to different effective values of Z.
An alternative method for calculating Z which also does not require a partitioning of the charge density is shown in the main panel of figure 2. For each value of the supercell charge ∆Q, we determine μ from the spectrum of Kohn-Sham energies (open squares) and then fit a linear-bands model of doped graphene of the f o r m ( Q Z is the total free carrier charge and A is the area of the supercell. This method yields = ± Z 1.6 0.3 e, significantly higher than previously suggested values [15], but in excellent agreement with our multiscale method. Figure 3 shows the LDOS obtained from the largescale TB simulation with the screened potential from the NLTF + inter model. The main panel shows a representative example for the case of n-doped graphene (at a chemical potential µ = 0.12 eV).

Local density of states
As observed in previous studies employing an unscreened impurity potential [10,19], the LDOS near the adatom exhibits a broken electron-hole symmetry. This behavior can be understood within a simple TFbased picture, in which the LDOS in the presence of the adsorbate is given by that of unperturbed graphene, but shifted in energy by the value of the screened impurity potential (shown by the dotted line in figure 3). While the shifted LDOS agrees well with the full result for energies more than ± ∼0.5 eV from the Dirac point, there are significant deviations within this range. In particular, the position of the Dirac point remains pinned at its unperturbed value, which is a consequence of the linear dispersion of the bands from graphene's chiral Dirac fermions.
The inset of figure 3 shows the decay of the LDOS at a fixed energy towards its unperturbed value far from the adatom (corresponding to experimental / I V d d linescans [10]). The decay length decreases with increasing µ . This is in agreement with very recent experimental measurements [28]. The chemical potential only enters the simulation through its effect on the screened potential; this demonstrates the impact of screening on measured properties of the system and highlights the importance of an accurate description of electron-electron interactions. 1 However, the four models result in somewhat lower overall values of F than fits based on empirical forms for the screened potential, discussed in figure S3 in the supplementary material. Figure 3. LDOS of free-standing, n-doped graphene at different distances from a single Ca adatom, calculated with the firstprinciples multiscale approach described in the text. The dotted line shows the LDOS of unperturbed graphene shifted in energy by the average value of the on-site screened potential at a distance of 10 Å from the adatom. The inset shows the LDOS at 0.1 eV above the Dirac point as a function of the distance r from the adatom (normalized by the value for unperturbed graphene), and as the graphene doping is increased (the chemical potential μ is given in the key).

Conclusions
In conclusion, we have introduced an accurate and efficient multiscale approach for the theoretical description of charged impurities on doped graphene. Starting from linear-scaling DFT simulations of a Ca adatom in a × 56 56 graphene supercell, we parametrize a continuum screening model by fitting the screened potential of the model to the DFT selfconsistent potential. The continuum screening model is then evaluated in a much larger × 168 168 supercell and the resulting screened impurity potential is used as input for large-scale TB simulations, which yield observable quantities, such as the local density of states.
Although the continuum screening model is nonanalytic and requires a self-consistent calcul ation, its computational cost is negligible, many orders of magnitude less than the TB simulation. The TB simulation itself is also computationally inexpensive: the × 168 168 supercell with TB is ∼0.1% of the cost of the × 56 56 supercell with DFT, allowing us to access much larger systems than would be possible with DFT alone. As shown in figure S4 in the supplementary material, the limiting factor in the accuracy of the large-scale simulation is the TB approx imation itself rather than our fitting procedure.
This multiscale approach, which bridges local chemistry at the adsorbate site on the nanoscale and long-range screening effects on the mesoscale within a single, physically-motived framework, opens up the possibility of investigating more realistic and complex mescoscale configurations relevant for the understanding of real devices, e.g. disordered arrangements of different adsorbate species.
Data underlying the DFT calculations are available on figshare [41] and may be used under the creative commons Attribution licence. The code for calculating the screened potential of such systems with the NLTF + inter model is available online with a userfriendly web-based interface [39].