Multi-scale approach to first-principles electron transport beyond 100 nm

TL;DR

This paper introduces a multi-scale computational method combining DFT and tight-binding models to simulate electron transport in large 2D devices over 100 nm, capturing atomistic details and complex device features.

Contribution

The paper presents a novel multi-scale approach that integrates DFT and tight-binding models for large-scale 2D electron transport simulations beyond 100 nm.

Findings

Successfully applied to pristine, defected, and nanoporous graphene devices

Enabled atomistic-level current propagation analysis in large devices

Addressed practical implementation challenges of multi-scale modeling

Abstract

Multi-scale computational approaches are important for studies of novel, low-dimensional electronic devices since they are able to capture the different length-scales involved in the device operation, and at the same time describe critical parts such as surfaces, defects, interfaces, gates, and applied bias, on a atomistic, quantum-chemical level. Here we present a multi-scale method which enables calculations of electronic currents in two-dimensional devices larger than 100 nm$^2$, where multiple perturbed regions described by density functional theory (DFT) are embedded into an extended unperturbed region described by a DFT-parametrized tight-binding model. We explain the details of the method, provide examples, and point out the main challenges regarding its practical implementation. Finally we apply it to study current propagation in pristine, defected and nanoporous graphene devices, injected by chemically accurate contacts simulating scanning tunneling microscopy probes.

Figures

Equations

$$\{{\bar{\varphi }}_{\alpha }\}=P\{{\varphi }_{\alpha }\}.$$

$${\mathcal{T}}_{ij}=\mathrm{Tr}\left[{\mathbf{G}}_D{\mathbf{\Gamma }}_i{\mathbf{G}}_D^{†}{\mathbf{\Gamma }}_j\right]$$

$${\mathbf{\Gamma }}_j=i\left({\mathbf{\Sigma }}_j-{\mathbf{\Sigma }}_j^{†}\right)$$

$${\mathbf{G}}_D={\left[\mathbf{S}(E+i\eta )-\mathbf{H}-\sum _i^N{\mathbf{\Sigma }}_i\right]}^{-1}$$

$$\mathbf{H}+\mathbf{V}=\left[\begin{array}{cc}{\mathbf{H}}_{1,1}& 0\\ 0& {\mathbf{H}}_{2,2}\end{array}\right]+\left[\begin{array}{cc}0& {\mathbf{V}}_{1,2}\\ {\mathbf{V}}_{2,1}& 0\end{array}\right],.$$

$$\left[z1-\mathbf{H}-\mathbf{V}\right]\mathbf{G}(E)=1,$$

$$\mathbf{G}={\mathbf{G}}^0+{\mathbf{G}}^0\mathbf{VG},$$

$${\mathbf{G}}_{1,1}$$

$${\mathbf{G}}_{1,1}^0$$

$${\mathbf{G}}_{2,1}$$

$${\mathbf{G}}_{1,1}$$

$$={\left[z\mathbf{I}-{\mathbf{H}}_{1,1}-{\mathbf{V}}_{1,2}{\mathbf{G}}_{2,2}^0{\mathbf{V}}_{2,1}\right]}^{-1}$$

$$={\left[z\mathbf{I}-{\mathbf{H}}_{1,1}-{\mathbf{\Sigma }}^{\prime }(E)\right]}^{-1},$$

$${\mathbf{\Sigma }}^{\prime }(E)={\mathbf{V}}_{1,2}{\mathbf{G}}_{2,2}^0(E){\mathbf{V}}_{2,1},$$

$${\mathbf{\Sigma }}^{\prime }(E)=({\mathbf{V}}_{1,2}-E{\mathbf{S}}_{1,2}){\mathbf{G}}_{2,2}^0(E)({\mathbf{V}}_{2,1}-E{\mathbf{S}}_{2,1}).$$

$${\mathbf{G}}_D={\left[\mathbf{S}(E+i\eta )-\mathbf{H}-\sum _i{\mathbf{\Sigma }}_i-\delta \mathbf{\Sigma }\right]}^{-1}$$

$${\mathbf{\Sigma }}_{\mathrm{avg}}=⟨\mathbf{S}{⟩}_k(E+i\eta )-⟨\mathbf{H}{⟩}_k-⟨\mathbf{G}{⟩}_k^{-1},$$

Full text

I Introduction

Developing high-performance computational strategies to simulate electronic devices is a fundamental asset for prototype design and research planning in basically any technological context.Hofer2018 ; Marzari2016 ; Thijssen2007 Being able to model nanometer-scale devices with atomic resolution has become particularly crucial for e.g. novel low-dimensional materials, molecular junctions, or, generally, ballistic quantum systems which have appeared on the electronics horizon in the last decades.Aprojanz2018 ; Shulaker2013 ; Ferrari2015 ; Novoselov2016 ; Geim2013 ; Martinez-Blanco2015 ; Moreno2018 ; Jia2016 ; Cao2018 Quantum-chemical details are often critical for describing local electronic structure, chemical contacts (e.g. electrodesKretz2018 ; Palacios2008 ) and electrostatics in realistic nano-electronic devices. This is especially clear in a context of one-atom-thick two-dimensional (2D) devices,Fiori2014 ; Novoselov2016 ; Iannaccone2018 where details on the atomic scale govern the electronic behavior.Schedin2007 ; Wehling2008 ; Chang2012 ; Caridad2018 ; talbotNPG These details can be accessed using various ab-initio methods,Szabo1989 ; Levine2014 ; Thijssen2007 in particular density functional theory (DFT),Kohn1965 ; Hohenberg1964 which can be used to simulate up to thousands of atoms.Goedecker1999 ; Mohr2015 ; Lin2014 This is typically enough to simulate isolated portions of a realistic device, such as bulk regions, interfaces or locally perturbed areas, but it is not suitable to simulate all the different length-scales involved in the operation of realistic devices.

A promising solution which is rapidly growing along many lines of research is to perform hybrid multi-scale simulations.Hofer2018 ; Hofer2015 ; Troisi2006 ; Kubar2010 ; Zhang2018 ; Shen2016 ; Fiori2013 ; AufDerMaur2011 ; Fediai2016 ; Fediai2016b In this context an interesting approach is to treat crucial parts of a relatively large device using full quantum-chemical detail, while reducing the number of degrees of freedom to the bare essentials elsewhere.Hofer2018 ; Hofer2015 ; Troisi2006 ; Kubar2010 ; Zhang2018 ; Shen2016 ; Ng2008 ; Rapacioli2011 A very popular example, coming from a biomolecular context, is the quantum mechanics/molecular mechanics (QM/MM) technique,Warshel1976 which has also been generalized to study solid-state surfaces and interfaces.Hofer2015 The key problem of multi-scale approaches is to partition the system into a number of subregions and, most importantly, to ensure a smooth, physically sound, transition among them.Hofer2018

Electronic structure methods employing linear combinations of atomic orbitals (LCAO) provides an intuitive route for elaborating multi-scale approaches analogous to QM/MM for simulating electron transport, as the accuracy and scalability of LCAO-based approaches can be “tuned” via clever approximations of the Hamiltonian. One common case is the Density Functional Tight-Binding (DFTB), which is based on a Taylor series expansion of the Kohn-Sham DFT total energy.Elstner1998 ; Elstner2014 Once a reliable element or compound-specific parametrization has been obtained, either comparing to experiments or higher-level ab-initio calculations, DFTB can be two to three orders of magnitude faster than DFT,Gaus2013 overall enabling simulations with several thousand atoms without need for any massive parallelization.Elstner2014 Another useful approximation is the Wannier Tight-Binding (WTB), where maximally localized Wannier functions (MLWF) are constructed ad-hoc to reproduce electronic bands from plane-wave DFT within optimized energy windows.Marzari1997 ; Marzari2012 ; Bruzzone2014 ; Marin2018 ; Lv2015 ; Pizzi2016 Once a system-specific parameterization is found, and good initial trial orbitals are defined,Gresch2018 one can use WTB to access samples containing hundreds of thousands of atoms.Rudenko2015 The latter two tight-binding (TB) like approximations can be ultimately simplified by reducing the number of orbitals per atom to a minimum, and limiting the interaction range to only a few nearest neighbors in the lattice, fitting parameters by hand to experiments or ab-initio calculations.Martin2004 ; Ashcroft ; Liu2015 ; Beconcini2016 ; Calogero2018 Multi-scale approaches for electron transport based on these models have been proposed and include combinations of Classical Molecular Dynamics (CMD) and DFT with Langevin Dynamics,Troisi2006 combinations of CMD, DFTB, and wave function propagation,Kubar2010 or combinations of TB models with patched Green’s function techniques.Settnes2015

In this article we present a multi-scale scheme based on seamless integration of a number of perturbed regions described by DFT and an unperturbed extended region represented by a simpler representative LCAO model, namely TB. We parametrize the TB Hamiltonian directly from DFT, such that the resulting multi-scale model is defined without any fitting parameter or by hand adjustments. The obtained models can be combined with the non-equilibrium Green’s function (NEGF) formalismTaGuWa2001 ; Brandbyge2002 ; Datta2000 ; Papior2017 to enable current simulations of over $100\mathrm{nm}\,\times 100\,\mathrm{nm}$ large devices with local quantum-chemical detail. The manuscript is divided into two main sections. In the first section we develop the general formalism, didactically providing concrete examples based on simple graphene two-electrode devices. In the second section we apply the multi-scale approach for imaging real-space “far-field” currents (i.e. far from the source) in pristine, defected and nanoporous graphene, injected by chemically accurate contacts simulating Scanning Tunnel Microscopy (STM) probes.

II Method

II.1 TB models from orbital-projected DFT

Our starting point is a pristine structure of a periodically repeated unit cell, illustrated in Fig. 1, and a DFT Hamiltonian describing this system with a localized basis. We will here consider a LCAO basis for the DFT calculation, $\phi_{\alpha}$, but it could equally well be a basis obtained from e.g. maximally localized Wannier states. We will from this construct a smaller TB-like basis, which only describe the bands in a region around the Fermi energy, $E_{F}$, set via a projection $P$,

[TABLE]

In the simplest case, used here, this involves selecting a particular subset of the original DFT basis functions, however one may imagine more involved projections. Our prime example is shown in Fig. 1b, where we consider the bandstructure of graphene obtained with DFT and a LCAO single-$\zeta$ plus polarization basis (i.e. a single orbital for each $s$, $p$ states and polarization $d$ orbitals). A TB-like model can be readily constructed by projecting on to the $p_{z}$ orbitals of the DFT LCAO Hamiltonian (and overlap), by extracting all rows and columns associated to $p_{z}$ orbitals. This cheaper model is sufficient to capture the $\pi$ bands of graphene (Fig. 1c). Likewise, the $\sigma$ bands further away from the Dirac point energy can be captured using a TB model parametrized from $s$, $p_{x}$ and $p_{y}$ DFT orbitals (Fig. 1d). Importantly, as highlighted (dashed lines) in Fig. 1c, a partial $p_{z}$-projection where only couplings among nearest neighbors are retained results in significant band-misalignment and rescaling. Analogous but less dramatic deviations are observed when using a standard nearest neighbor model with orthogonal basis and hopping $t=2.7\,\mathrm{eV}$.

II.2 Multi-scale approach

In the following we omit energy and $\mathbf{k}$ dependence, highlighting it only where necessary.

The multi-scale method presented in this work is based on the Non-Equilibrium Green’s Function (NEGF) transport formalism.TaGuWa2001 ; Brandbyge2002 ; Datta2000 ; Papior2017 In the NEGF framework transmission between any two leads $i$ and $j$ of a $N$-electrode device with Hamiltonian $\mathbf{H}$ and overlap $\mathbf{S}$ is given by

[TABLE]

where the device Green’s function $\mathbf{G}_{D}$ is given by

[TABLE]

and $\boldsymbol{\Sigma}$ are the so-called self-energies of the semi-infinite electrodes.

It is the self-energies that will play a pivoting role in connecting TB and DFT models. The self-energy is nothing but the effect of degrees of freedom not accounted explicitly for in the equation but have been eliminated as can be done for variables in linear equations. The self-energies can represent the semi-infinite electrodes as well as a finite number of degrees of freedom in the case of a finite region. In the following we show that it is in practice possible to locally replace small perturbed regions of a large TB device with DFT-precision models by simply including one or more “special” self-energies in the sum of Eq. (4).

II.2.1 Self-energy for a partitioned system

To understand the construction of these special self-energies it is instructive to recall the definition of self-energy connecting two subregions of a generic binary system. We do this by simply generalizing the derivation of self-energy for a simple system divided into two parts.PapiorThesis

Let us consider the system in Fig. 2a, with Hamiltonian $\mathbf{H}$, and overlap $\mathbf{S}\equiv\mathbf{1}$ (we use an orthogonal basis without loss of generality),

[TABLE]

Here we have explicitly indicated with $\mathbf{V}$ the coupling of region 1 ($\mathbf{H}_{1,1}$) to the neighboring region 2 ($\mathbf{H}_{2,2}$). In the following we are only interested in the Green’s function in region $\mathbf{G}_{1,1}$, and thus a derivation of this matrix is required. The full Green’s function is

[TABLE]

with $z=E+i\eta$ and $\eta\to 0^{+}$. By using the Dyson equation

[TABLE]

we can express the Green’s function in only one region

[TABLE]

and thus

[TABLE]

Here the term

[TABLE]

is the self-energy describing how region 1 is perturbed by the coupling to the degrees of freedom in region 2. In case of a non-orthogonal basis set ($\mathbf{S}\nequiv\mathbf{1}$) Eq. (14) simply becomes

[TABLE]

The real part of the self-energy is an energy renomarlization/shift, while the, possibly finite, imaginary part corresponds to a finite life-time or broadening of the energy levels in region 1. Additionally the self-energy $\boldsymbol{\Sigma}^{\prime}$ is independent of $\mathbf{H}_{1,1}$ and may thus be used in any other system $\mathbf{H}_{\rm ext}$ so long as the coupling matrix is unchanged, i.e. $\mathbf{V}_{1,2}=\mathbf{V}_{{\rm ext},2}$.

In Fig. 2a, we see the above exemplified using an extended system with a hole perturbing the local potential (dark colors). After two layers the perturbation goes to zero and the coupling between neighbouring cells turns to the bulk value($\mathbf{V}$). At this point one can calculate $\boldsymbol{\Sigma}^{\prime}$ for region 2, insert it into an external pristine system (equal to region 1 in Fig. 2b) and, in turn, reproduce the correct Green’s function in region 1, as though the new system was connected with the hole.

II.2.2 Example: graphene device with DFT-DFT connection

In the following we illustrate with a concrete example where the two models, connected by $\boldsymbol{\Sigma}^{\prime}$, are constructed using the same method (here DFT with a SZ basis set, i.e.  4 orbitals per C atom), then effects of the perturbation can be propagated from one model to the other with 100% accuracy. Consider the basic DFT two-electrode graphene device illustrated in Fig. 2c (upper inset). This has periodic boundary conditions along $y$, semi-infinite electrodes $L$ and $R$ along $\pm x$ and a hole in the scattering region. Total transmission $\mathcal{T}_{LR}$ across the device can easily be computed by using Eq. (2) and is plotted in Fig. 2c. One can choose to split this device into two sections, a bulk part (region 1) and a perturbed region 2 (hole). Since the potential of the hole is fully screened in region 1, we have that the coupling between $\mathbf{V}_{2,1}$ equals the bulk coupling, $\mathbf{V}$. Hence one can construct the self-energy $\boldsymbol{\Sigma}^{\prime}_{R^{\prime}}$ which contains the effects of region 2 propagated into region 1 using Eq. (14). This can then be incorporated as a new electrode on the right side of an external DFT device which has no holes in the scattering region (lower inset) and verify that $\mathcal{T}_{LR^{\prime}}(E)=\mathcal{T}_{LR}(E)$ for all $E$ within the numerical accuracy.

II.3 Self-energy for DFT-TB connections

We now turn to the more interesting situation where the perturbation and the external unperturbed system hosting its self-energy $\boldsymbol{\Sigma}^{\prime}$ are modelled using different basis sets, e.g. DFT and TB. We have already anticipated with the example in Fig. 1 that the electronic structure of a DFT system within a particular energy window often can be reproduced by a small TB model obtained as a projection of the DFT Hamiltonian. The drastically reduced number of orbitals in a TB model makes it potentially very convenient for generating unperturbed external host systems with large dimensions, normally inaccessible by DFT calculations. By generalizing the definition of the self-energy, $\boldsymbol{\Sigma}^{\prime}$, it becomes possible to incorporate DFT-precision perturbations in TB models. However determining a self-energy for a DFT-TB connection is not as trivial as for the basic DFT-DFT connection described in the previous section.

Consider now the system as shown in Fig. 3a. Instead of having both regions (1 and 2) using the same basis set, we change the basis set in region 1 to be a TB parameterized basis set, $P\{\phi_{\alpha}\}=\{\bar{\phi}_{\alpha}\}$, as discussed in Fig. 1. Here the approximation lies in $\mathbf{V}_{\mathrm{DFT-TB}}$. By the same arguments outlined in the theory section, we create the subset of $\mathbf{V}_{\mathrm{DFT-TB}}$ such that the rows correspond to the full basis set $\{\phi_{\alpha}\}$ and the columns correspond to the projection orbitals $P\{\phi_{\alpha}\}$. Using $\mathbf{V}_{\mathrm{DFT-TB}}$ in Eq. (14) results in a projection of the self-energy onto the parameterized orbitals. Importantly, the fact that some elements of the original DFT coupling, $\mathbf{V}$, are now missing in $\mathbf{V}_{\mathrm{DFT-TB}}=P\mathbf{V}$ inevitably leads to some scattering at the DFT-TB boundary.

II.3.1 Example: graphene device with DFT-TB connection

We have demonstrated that a TB model fully parametrized from DFT can reproduce the $\pi$ bands of graphene with only minimal deviations (Fig. 1c). As a result DFT results can be reproduced, within accuracy, using the multi-scale DFT-TB approach. Let us consider again the graphene device with a hole, as shown in Fig. 3a. The only difference between DFT-DFT and DFT-TB is that here $\mathbf{V}_{1,2}$, in Eq. (14), is chosen such that the columns of $\mathbf{V}_{1,2}$ project onto the parameterized $p_{z}$ orbitals. This yields a self-energy only existing on the $p_{z}$ orbitals in the TB model. However, one may additionally play with the number of neighbors each atom connects with, e.g. 1, 2, 3 or all neighbors (with “all” defined by the DFT basis).

In Fig. 3c the transmission for the hole system is shown at four levels of precision compared to the full DFT calculation. We find that if the TB parametrization contains only DFT couplings among nearest neighbors then the transmission spectrum undergoes significant scattering $\sim$0.1 eV away from the Dirac point. When the range of interaction is extended to 2, 3 or 4 neighbors the transmission spectrum becomes comparable to DFT in a wider energy range, except for a shift at positive energies. As soon as couplings among 5 nearest neighbors are included we obtain an almost perfect agreement with DFT.

II.4 Self-energy of isolated perturbations

In all the examples considered so far periodic boundary conditions was employed in the DFT calculations from which the connecting self-energy $\boldsymbol{\Sigma}^{\prime}_{\mathrm{DFT-TB}}$ was constructed. This means that the perturbation, e.g. the hole, is periodically repeated along the transverse direction $y$. However, if the potential is screened to the bulk value one may replace the surrounding periodic images by other environments with the same potential and Hamiltonian. In this way one may effectively change the boundary conditions.

The origin of this versatility lies in the fact that the two regions in which the perturbed DFT model is divided can have arbitrary size, shape and periodicity. For example, with reference to Fig. 3, one could choose to define the pristine DFT region 1 as the outermost frame-shaped area of the cell surrounding the perturbation, while treating all the orbitals enclosed by it as region 2. The same formalism discussed above can be readily applied to this case, with the only difference that periodic boundary conditions are removed from the DFT Hamiltonian before constructing the self-energy $\boldsymbol{\Sigma}^{\prime}_{\mathrm{DFT-TB}}$.

An example is illustrated in Fig. 4. The top geometry is the DFT region which is used in calculating $\boldsymbol{\Sigma}^{\prime}_{\mathrm{DFT-TB}}$. The bottom geometry is the TB parameterization and the red atoms indicate the overlay region where the self-energy is transferred from region 1 (DFT) to region 2 (TB). We point out that this particular approach provides the basis for “modular” multi-scale simulations, where multiple DFT-precision perturbations, modules, are incorporated into the same large TB device. We will present concrete applications of this in the second part of the work.

II.5 Challenges and implementation

Despite reflecting a conceptually simple theory, the protocol to embed DFT-precision perturbations locally into an unperturbed TB region parametrized from DFT implies a number of critical issues when it comes to its practical implementation. Let us briefly go through the procedure to include non-periodic DFT perturbations into large TB models step-by-step:

Generate a DFT model of the perturbation. This can be represented by either a bulk (periodic) or device (semi-infinite electrodes) setup and must be large enough to ensure the potential turns constant in its outermost areas; 2. 2.

Define two regions in the DFT system, one containing all orbitals involved in the perturbation and the other involving the projection orbitals. The latter needs to be as far as possible from the perturbation so as to ensure $\mathbf{V}_{\mathrm{DFT-TB}}=\mathrm{constant}$; 3. 3.

Represent the DFT Hamiltonian and overlap in real-space (without periodic boundary conditions); 4. 4.

Compute the self-energy $\boldsymbol{\Sigma}^{\prime}$ using equation Eq. (14); 5. 5.

Generate a DFT model of the same system, but without the perturbation; 6. 6.

Construct a large TB model using parameters from the projection of the unperturbed DFT calculation. 7. 7.

Incorporate the self-energy $\boldsymbol{\Sigma}^{\prime}$ locally into the TB model using Eq. (4).

To begin with, most of the described steps require flexible manipulation of the DFT Hamiltonian and overlap matrices. Extraction of selected on- and off-diagonal elements exclusively associated to a subset of orbitals is crucial to compute $\mathbf{G}_{2,2}$, extract $\mathbf{V}_{\mathrm{DFT-TB}}$ or to create TB models parametrized from DFT. Secondly, to ensure the correct couplings at the DFT-TB interface when incorporating $\boldsymbol{\Sigma}^{\prime}$ into the TB model, it is important that region/atoms selected to host $\boldsymbol{\Sigma}^{\prime}$ in the large TB geometry have a one-to-one correspondence with those in region 1 of the DFT geometry. User-friendly tools to access and compare lattice coordinates of the various involved geometries are therefore highly desirable.

Furthermore, the self-energy, $\boldsymbol{\Sigma}^{\prime}$, needs to be stored into external files with general and compact format, readily accessible for usage in NEGF calculations. Being able to flexibly input/output self-energies into external host models is especially crucial to carry out multi-scale calculations where several DFT-precision regions are accounted for in the same large TB device.

Overall, the computation of $\boldsymbol{\Sigma}^{\prime}$ through Eq. (14) and NEGF transmissions needs to be optimized both time and memory-wise, so as to efficiently handle dense DFT matrices and large sparse TB ones, as well as calculation of the self-energy for a fine $k$ grid and several energy values.

Computational tools

We tackle all the issues described in the previous section with open-source tools TBtrans and sisl.Papior2017 ; Papior2017a TBtrans is distributed with the Siesta/TranSiesta software package.Soler2002 ; Papior2017 TranSiesta enables high-performance DFT+NEGF self-consistent calculations of large, multi-terminal systems under various electrostatic conditions (e.g. gating.Papior2015 ) TBtrans is a post-processing NEGF code which provides a flexible interface to DFT or TB Hamiltonians and enables large-scale calculations of spectral physical quantities, interpolated $I$-$V$ curves and/or orbital/bond-currents for systems easily exceeding millions of orbitals on few-core machines. sisl is a Python package used to create and/or manipulate DFT and TB models for arbitrary geometries, with any number of orbitals and any periodicity.Calogero2018 The device Green’s function in TBtrans is generally implemented as:

[TABLE]

where $\mathbf{S}$ and $\mathbf{H}$ are DFT or TB-modeled overlap and Hamiltonian, $\boldsymbol{\Sigma}_{i}$ are the electrodes self-energies and $\delta\boldsymbol{\Sigma}$ is a user-defined, optional, perturbative term. The flexibility of using sisl and TBtrans is appreciated since no code extension of TBtrans is required as everything being done is added using $\delta\boldsymbol{\Sigma}$ and $\boldsymbol{\Sigma}_{i}$ terms. Implementation of these two terms in Eq. (16) is done only in Python and is thus much simpler than a full fortran implementation. Secondly, both sisl and TBtrans are scaling easily to millions of orbitals in the TB method.

All of the calculations described in this article have been carried out using a single high-performance computing node with 20 cores (Core Intel Xeon E5-2660v3, 2.6GHz, 128GB RAM), parallelizing TBtrans and either MPI or OMP parallelization. The transmission spectra in Fig. 2 and Fig. 3 were calculated in TBtrans using 100 energy points, 201 k-points in the periodic direction and a total of 540 atoms. The computationally most expensive calculation presented so far, i.e. , that involving a DFT-DFT connection, required 3 hours and 50 minutes to compute/store the connecting self-energy $\boldsymbol{\Sigma}^{\prime}_{\mathrm{DFT-TB}}$, and 10 minutes (53MB RAM, 113kB disk) to calculate transmission with TBtrans.

III Applications

In the second part of this work we show how the multi-scale approach can be used to study the far-field behavior of electrons injected from atomic-scale contacts into large-scale graphene-based devices.

Being able to interpret and predict the behavior of electrons over large-scale devices is of importance for graphene-based electronics, especially in the growing field of electron optics. Real-space visualization of charge (or spin) transport can be achieved experimentally using various quantum imaging techniques, including probe microscopy,Hiraoka2017 ; Kolmer2017 superconducting interferometryAllen2015 and magnetometry with diamond-NV centers.Casola2018 It has been shown that defects and contacts with tips at the atomic scale yield strong spatial variations of current flow in graphene devices.Bhandari2016 ; Tetienne2017

Modeling realistic atomic-scale contacts to inject currents in graphene devices, while simultaneously accessing current flow far away from them, is a challenge for state-of-the-art atomistic transport calculations. Large graphene devices can be accessed using the TB approximation. Despite allowing for good scalability, in this context the complex chemical nature of defects or tip contacts is often subject to drastic simplifications. For instance a STM probe is usually modeled as a constant, on-site, level broadening ($i\Gamma$) self-energy term in the Green’s function,Aprojanz2018 ; Settnes2014 a localized effective force or electrostatic field,Eder2013 ; Mark2017 or even as a narrow semi-infinite in-plane electrode.Calogero2018 On the other hand, DFT-based models have often demonstrated to successfully corroborate probe microscopy measurements on graphene, even when involving complex tip functionalizations or inelastic effects.Jelinek2017 ; Frederiksen2011 ; Palsgaard2015 ; Heijden2016 This is due to an accurate description of tip structure and orbital symmetries, as well as charge and potential variations in the sampled regions.

The essential need for both scalability and accuracy in this problem calls for the multi-scale approach described above. In the following we will provide a bird’s eye view on carrier injection from point sources into pristine, defected and nano-porous graphene (NPG). Furthermore, In order to further emphasize the versatility of the method, we will show an example where multiple DFT-precision regions from various DFT calculations are included in the same large-scale TB transport calculation.

III.1 Far-field currents in graphene

Here we use the multi-scale approach to predict the far-field behavior of electrons injected by a DFT-precision gold tip into a large-scale TB model of pristine graphene. Using wave-packet dynamics simulations to inject currents from an effective electrostatic field model of STM tip, Mark et al. reported anisotropic electron current in graphene.Mark2012 The authors predicted that a well-defined six-fold asymmetry current pattern occurs when the wave packet is injected in the middle of a hexagon in the graphene lattice.Mark2012 Below we confirm this anisotropic signature by looking at the far-field currents flowing up to 10 nm away from the injection. The system is modeled by TB and is coupled to a small DFT-precision injection region through a self-energy $\boldsymbol{\Sigma}_{\mathrm{DFT-TB}}$. We will limit our focus to a single representative energy value where this six-fold anisotropy is clearly visible.

We start by setting up a three-electrode DFT model device (see Fig. 5) consisting of a semi-infinite gold tip placed $\approx 2.0\,\textup{\AA}$ above the center of a graphene hexagon and two semi-infinite graphene electrodes along $y$. The tip structure is chosen so as to ensure a flat local density of states on the tip apex in a wide energy interval, $[-0.75,\,1]\,\mathrm{eV}$. We optimize the tip apex atom and the nearest $\sim 20$ C atoms in Siesta until forces are less than $0.01\,\mathrm{eV}/\textup{\AA}$, using periodic boundary conditions along $x$ and $y$, a $3\times 3$ $k$-point Monkhorst-Pack grid, the GGA-PBE exchange-correlation functional,Perdew1996 a SZP basis set and an energy cutoff of 300 Ry. The tip is in chemical contact to graphene, therefore we assume van der Waals interactions to be negligible and do not include them in the calculations. We then calculate bond-currents flowing across the three-electrodes device using the tip as a source, using $30$ $k$-points to sample the $x$ direction in TBtrans. In order to avoid artificial crosstalk with periodic images of the tip when computing the injected bond-currents, we substitute the self-energies for the two graphene electrodes with a new self-energy term $\boldsymbol{\Sigma}_{\rm avg}$ in Eq. (4),

[TABLE]

where $\langle\cdot\rangle_{k}$ represents an average over all $k$-points along $x$, i.e. the quantities in the principal unit cell in real-space.Papior2019 In effect, this is equivalent to setting up a new drain electrode in the border region of the cell.

The results for $E-E_{\rm F}=0.8\,\mathrm{eV}$ are summarized in Fig. 6. We find that the total transmission from tip to $\boldsymbol{\Sigma}_{\rm avg}$ is $\mathcal{T}=0.686$. This is determined mainly by three transmission eigenchannelsPaulsson2007 . The two with largest contribution (37% and 34% of the total transmission) are mostly delocalized over the whole graphene structure, whereas the third one (contributing with 29%) exhibits a preferential propagation along the six armchair lattice directions departing from the probed graphene hexagon (Fig. 6a). We find that this six-fold anisotropy dominates the bond-current pattern injected by the tip at this energy (Fig. 6b).

Then, using the approach discussed in Fig. 4, we insert the DFT-precision injection region inside a larger TB model parameterized from DFT, using a complex absorbing potentialCalogero2018 ; Xie2014 (CAP) to absorb currents at the cell boundaries. As shown in Fig. 6c we find that the six-fold anisotropic propagation can still be observed far away from the source, although collimation is rapidly lost. As a check we estimate transmission outgoing from the tip by summing over positive bond-currents which cross a circle of radius $R$ centered at the tip $(x,y)$ coordinates. The result, $\mathcal{T}=0.672$, is only slightly lower than the value obtained by summing bond-currents in the smaller all-DFT system, $\mathcal{T}=0.681$, regardless the choice of radius $R$.

Importantly, we find that the result does not change if the shape of the DFT region in which $\boldsymbol{\Sigma}_{\mathrm{DFT-TB}}$ is computed is modified from a rectangular to a circular one (Fig. 7a).

Effects of single N dopant on far-field currents

One major advantage of the multi-scale approach is that complex situations, e.g. where different perturbations simultaneously affect the same region, can be described with chemical accuracy. An illustrative example is depicted in Fig. 7b, where the multi-scale method has been used to study how a N substitutional dopant in proximity of the Au tip contact affects the far-field currents. In order to unravel possible effects of spin polarization we carry out this calculations for its two possible spin configurations. We find that in both cases, due to the local doping induced by the N atom, electrons injected by the tip are prevented from propagating towards the N atom (Fig. 7c).

TB parameterization

We point out that, compared to the simple case discussed in Fig. 3, the presence of a metal contact in these systems alters the graphene work-function, and hence its potential and couplings far from the tip. As a result, parameterizing the TB model directly from a DFT calculation of pristine graphene would induce some degree of coupling mismatch at the interface between the two models. The solution we adopt here to ensure constant coupling $\mathbf{V}_{\mathrm{DFT-TB}}$ is to parameterize the TB model from the pristine-like DFT elements associated to the atoms far away from the tip, namely those further than $\approx 17\,\textup{\AA}$ from it (Fig. 7c).

III.2 Far-field currents in nanoporous graphene

The principles behind the multi-scale method are general and not limited to the basic graphene systems presented above. For instance, let us now consider nanoporous graphene (NPG), which is inherently semiconducting and anisotropic.Pedersen2008 In particular we consider the NPG structure shown in Fig. 8a, which has recently been fabricated through bottom-up synthesis with unprecedented sample sizes and quality.Moreno2018 In this system an electron wave is forced to channel into coupled 1D parallel pathways, in analogy to light waves propagating through an array of optical wave-guides. In a recent worktalbotNPG we have exploited the capabilities of the multi-scale approach to tackle two essential questions in this context, namely whether these transmission channels interfere with each other and how the wave profile can be tuned and controlled over long propagation distances. We have demonstrated that an electron wave injected into NPG from a model STM tip spreads over the nanomesh according to the same equations which govern the so-called Talbot optical interference effect (Fig. 8b-c). Further information can be found in Ref. [21].

Application of the multi-scale method in this context was vital, as the typical distances accessible by DFT are not long enough to capture the characteristic interference fringes (see Fig. 8b). Importantly, contrary to the graphene systems considered in the previous section, in this study we have ensured constant coupling $\mathbf{V}_{\mathrm{DFT-TB}}$ by introducing doping via an electrostatic gate $15\,\textup{\AA}$ below the NPG plane, such that $\pm 10^{13}\,e^{-}/\mathrm{cm}^{2}$ carriers are induced in the NPG.Papior2015 This ensures that the Fermi level is pinned throughout the whole DFT+TB device, effectively avoiding any possible artificial mismatch between the models involved in the calculation of $\boldsymbol{\Sigma}_{\mathrm{DFT-TB}}$.

III.3 Large-scale TB with N>1 DFT-precision regions

Another advantage of the multi-scale method is that one can accommodate any number of perturbed regions in the TB model, by simply including a corresponding number of self-energy terms in Eq. (4). This is exemplified in Fig. 9, where we show results from multi-scale TB+DFT calculations where multiple DFT-precision regions have been set up within a single large-scale TB model. Fig. 9a illustrates far-field currents from a multi-scale transport calculation of a $55\times 35\,\mathrm{nm}^{2}$ TB graphene device with seven self-energy terms. Two of these model regular graphene electrodes semi-infinite along $y$. The third one simulates a CAP in the outermost $5\,\mathrm{nm}$ of the cell along $x$. The remaining $N=4$ represent DFT-precision regions, namely two Au tips in contact to graphene ($\approx 2\,\textup{\AA}$ above an hexagon) and two epoxy groups, which preserve the $sp^{2}$ nature of grapheneSljivancanin2013 and are of particular interest for e.g. engineering the thermal conductivityZhao2015 or catalytic activitySinthika2015 of carbon/graphene-based systems. Both DFT systems are described using a SZP basis set and have been optimized with force threshold $0.01\,\mathrm{eV}/\textup{\AA}$. Self-energies are computed once for each of the two systems, while a bottom gate induces $-10^{13}\,e^{-}/\mathrm{cm}^{2}$ into graphene. Bond-currents injected at $E-E_{\rm F}=0.8\,\mathrm{eV}$ by the upper Au tip into the large TB model clearly show how the six-fold anisotropic electron wave scatters off the defects and the second tip. We point out that, once the corresponding self-energies have been computed and stored, the DFT-precision regions can be moved around the large-scale TB model very efficiently. In the Supplementary Material we show how the bond-currents landscape looks like for different positions of an epoxy defect relative to two STM tip contacts.

Another interesting situation, already discussed in Ref. [21] is depicted in Fig. 9b, with regard to the Talbot interference pattern discussed in the previous section. By inserting a second STM tip into the large-scale TB model it is possibletalbotNPG to map out the interference pattern observed in the bond-currents by performing a dual-probe experiment where a second STM tip is scanned over the NPG structure.

The computationally most expensive calculation presented in this section, i.e., the one illustrated in Fig. 8, required (for 150 energy points, $\Gamma$ only)  110 minutes to compute/store the connecting self-energy $\boldsymbol{\Sigma}^{\prime}_{\mathrm{DFT-TB}}$, and  4 hours (55GB RAM, 15GB disk) to calculate bond-currents with TBtrans.

IV Conclusions

In conclusion, we have presented a multi-scale method which enables calculations of local currents in devices larger than 100 $\mathrm{nm}^{2}$, by linking a perturbed region described by DFT to an unperturbed large-scale region described by an effective TB model parametrized from DFT. We have introduced the theory behind the method using basic concepts in the Green’s function framework, provided didactic examples and pointed out the main difficulties connected with its implementation. By applying the method to study realistic current injection by STM probes into pristine, defected and nanoporous graphene devices we have highlighted versatility, efficiency and generality of the method. Similar to hybrid QM/MM techniques, combining the advantages of DFT and TB methodologies using this scheme provides an adequate framework for embedding regions where accuracy is necessary into regions where size matters more. The largest system size accessible using this scheme obviously depends on the available computational resources and, most importantly, on the parametrization of the TB model, which varies with the material under study. We have illustrated the method in the simple case of carbon-based 2D materials, where the mapping from DFT to TB is a simple projection. Further studies may unveil whether other projectors are better or more generalizable. Overall, the most critical point is to define the TB to be computationally manageable, while coupling the pristine TB and DFT regions such that the interface scattering is much smaller than the scattering mechanisms under study. We further anticipate that the basic concepts of this method can be generalized to investigate thermal transport, where multi-scale approaches are often crucial.Mortazavi2017 ; Mortazavi2014 ; Nozaki2011 For this purpose the open-source PHtrans packagePapiorThesis (a TBtrans variant) may turn out to be very helpful.

Acknowledgments

Financial support by Villum Fonden (00013340), Danish research council (4184-00030) and EU-Horizon 2020 research and innovation programme (766726) is gratefully acknowledged. The Center for Nanostructured Graphene (CNG) is sponsored by the Danish Research Foundation (DNRF103).

V Supplementary Material

(Available online) Large-scale TB transport calculations of a $60\times 35\,\mathrm{nm}^{2}$ graphene device with N=3 DFT-precision regions, namely one epoxy defect and two source and drain Au tips, both in contact to graphene $\approx 2\,\textup{\AA}$ above an hexagon. Two regular graphene electrodes are used along $y$ and CAP is set at the outermost $5\,\mathrm{nm}$ on cell boundaries along $x$ to absorb currents. The file shows bond-currents injected by one of the two tips into the system, scattering off the drain tip and the epoxy defect, for different positions of the epoxy defect with respect to the two tips.