{'fraction_non_alphanumeric': 0.03000072817301391, 'fraction_numerical': 0.011723585523920483, 'mean_word_length': 4.560323886639676, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modified Dirac Hamiltonian for Efficient Quantum Mechanical Simulations of Micron Sized Devices

Representing massless Dirac fermions on a spatial lattice poses a potential challenge known as the Fermion Doubling problem. Addition of a quadratic term to the Dirac Hamiltonian circumvents this problem. We show that the modified Hamiltonian with the additional term results in a very small Hamiltonian matrix when discretized on a real space square lattice. The resulting Hamiltonian matrix is considerably more efficient for numerical simulations without sacrificing on accuracy and is several orders of magnitude faster than the atomistic tight binding model. Using this Hamiltonian and the Non-Equilibrium Green's Function (NEGF) formalism, we show several transport phenomena in graphene, such as magnetic focusing, chiral tunneling in the ballistic limit and conductivity in the diffusive limit in micron sized graphene devices. The modified Hamiltonian can be used for any system with massless Dirac fermions such as Topological Insulators, opening up a simulation domain that is not readily accessible otherwise.

Two-dimensional materials have attracted considerable recent attention for their superior electronic characteristics and potential for electronic, opto-electronic and spintronic applications 1,2 . Among them, graphene, topological insulators and transition metal dichalcogenides stand out in particular. Experimental progress on such materials are advancing rapidly 3 and modeling carrier transport to capture the new physics of these materials has become quite crucial. NEGF based numerical calculation is widely used nowadays to accurately model nanoscale materials 4 . The formalism is extremely powerful in solving any quantum transport problem accurately including sophisticated contact-channel effects and various forms of scatterings within a self-energy correction to the Hamiltonian. Despite these considerable strengths, NEGF has so far been considered mainly as a ballistic quantum transport simulation platform, since it becomes computationally prohibitive to model diffusive systems. This is especially true for experimentally relevant device dimensions which are often in the hundreds of nanometers to µm regime.
In materials such as graphene and topological insulators etc., electrons behave as massless Dirac fermions described by the Dirac Hamiltonian. Quantum transport simulations using the NEGF formalism require a real space representation of the Hamiltonian matrix to fully describe the channel material. While the tight-binding representation is valid for most practical purposes, it is computationally expensive. To expedite the calculation, a discretized version of the Dirac Hamiltonian can be used. It has however been shown that representing the Dirac Fermions on a spatial lattice poses a problem com-monly known as the Fermion Doubling 5,6 problem, where the numerical discretization of the Hamiltonian creates additional branches within the Brillouin zone (Fig. 1). One way to solve this problem is to add a quadratic term with the Hamiltonian 5,7 .
In this letter, we show that the additional term not only solves the Fermion Doubling problem but also has a profound implication for numerical simulations. We show that the modified Dirac Hamiltonian results in a spatial Hamiltonian matrix which is several orders of magnitude smaller than the atomistic tight binding Hamiltonian 8 while preserving the same level of accuracy over the relevant energy range. As a result, bandstructure and transport simulations using the modified Hamiltonian are several orders of magnitude faster than the tight binding, allowing us to carry out several large scale simulations such as angle dependent transmission in 1µm wide graphene (Fig. 2), magnetic electron focusing in a 0.4µm×0.4µm multi-electrode graphene (Fig. 3) and diffusive transport in a 1µm × 0.5µm graphene device (Fig. 4) . Our results show very good agreement with recently reported experimental results 9 . Other methods such as tight-binding make it almost impossible to simulate devices at such dimension.
The modified effective k.p Hamiltonian for graphene at low energy is, where, v F is the Fermi velocity, k = k xx + k yŷ is the wave vector, σ's are Pauli matrices representing the pseudospins, andh is the reduced Planck's constant. The last term, β(k 2 x + k 2 y )σ z serves two purposes: (i) it circumvents the well known Fermion doubling problem 6 and (ii) it allows us to generate a computationally efficient Hamiltonian using a course grid with reasonable accuracy as shown below. The k-space Hamiltonian in Eq. 1 is transformed to a real-space Hamiltonian by replacing k x with differential operator −i∂/∂x, k 2 x with −∂ 2 /∂x 2 and so on. The differential operators are then discretized in a square lattice using the finite difference method to obtain, (2) , a is the grid spacing and α ≡ β/a. When α = 0, Eqs. (1) and (2) reduce to the unmodified k.p Hamiltonian. The eigen-energy calculated from Eq. 2 becomes a sine function of wavevector k, and three extra Dirac cones appear inside the first Brillouin zone as shown in Fig. 1(a). This is known as the Fermion doubling problem 6 . The last term in Eq. 1 gets rid of this problem by opening bandgaps for each of these Dirac cones. The resulting band structure contains only one Dirac cone as shown in Fig. 1(b). The effect of the extra term in Eq.
(1) is more clearly illustrated in Fig. 1(c). When α = 0 and a = 20A o , the band structure along k y = 0 is a sine function and an extra Dirac cone appears at the zone boundary. This Dirac cone is removed when α = 1.165 and the band structure closely follows the ideal band structure of graphene. The band structure calculated using α = 1.165 with a = 20Å is accurate over a larger window than that is calculated using α = 0 and a = 10Å. Thus, the α parameter not only circumvents the Fermion Doubling problem but also enables us to use coarser lattice grid without losing the accuracy. With a = 20Å and α = 1.165, the Hamiltonian size given by Eq.
(2) for a 1µm × 1µm graphene sheet is 5E5 × 5E5 compared to ∼ 38E6 × 38E6 in atomistic tight binding model. The band structure calculated using this discretized Hamiltonian closely resembles the ideal linear band structure within |E| ∼ 0.6eV. To obtain the same level of accuracy, the recently proposed scaled graphene model 10 requires a Hamiltonian of size ∼ 9E6 × 9E6.
The accuracy of the proposed Hamiltonian is demonstrated by the band structure of a 200nm graphene ribbon shown in Fig 2. The bands calculated using the one p z orbital tight-binding model (in black) and the modified Dirac Hamiltonian (in blue) are in good agreement in the low energy limit. For the tight binding model, the calculation takes about 1 hour and 30 minutes. In comparison, the modified Hamiltonian model takes about 2 seconds using the same number of cores. This shows that our proposed model is three orders of magnitude faster with excellent accuracy for band structure calculation. With the modified Dirac Hamiltonian, we employ the standard Recursive Green's Function (RGF) algorithm 11 and employ it for large scale simulations, in both ballistic and diffusive regimes. Fig. 2b shows the conductance calculated using this model (with α = 1.16 and a = 16A 0 ) for a 1µm wide graphene sheet and compared with the conductance from the linear graphene E − K relationship. We see a good agreement until E F = 0.6eV. In Fig. 2c, we show the conductance of a electrostatically doped split-gated graphene pn junction, showing excellent agreement with the exact analytical solution 12,13 . That means the angle dependent transmission, an important property of graphene originating from its chiral nature 14 , is undistorted in the modified Hamiltonian. It can be shown that the pseudospins of the modified Hamiltonian are of the form, ψ = ψ 1 ψ 2 = 1 exp iθ /f (β, k F ) where f (β, k F ) = βk F + 1 + β 2 k 2 F , independent of the angle (θ = tan −1 k y /k x ) and therefore does not distort the chiral properties. Transverse magnetic field (TMF)s have been used in the past to study various transport phenomena and characterize surfaces and interfaces 15 .
In a recent experiment 9 , a TMF is used to focus electrons in a monolayer graphene device. The device geometry and biasing scheme used in our simulation are shown in Fig. 3(a). Electrons are injected from contact with a current source a and collected at contact b. In presence of a magnetic field B, an electron follows a circular path inside the graphene channel with cyclotron radius r c and is directly focused from contact a to contact c when 2r c = L, where L is the distance between contacts a and c. This is the resonance condition where the voltmeter registers a large voltage. Thus, the magnetic field required for the focusing is where n is the density of electrons in the graphene channel and m is an integer. For m > 1 the electron reaches contact c after skipping along the edges through multiple specular reflections. NEGF simulations for the multi-terminal device are shown in Fig. 3. The contacts are modeled using the self-energies of semi-infinite graphene ribbons. The current at contact α is calculated using the multi-terminal Landauer-Büttiker formalism 4, 16 , where f is the Fermi function, µ is the electro-chemical potential of the contact and T αβ is the total transmission between contact α and contact β. The transmission is calculated using the Fisher-Lee formula 4 , where G R is the retarded Green's function and Γ α,β = i(Σ α,β − Σ † α,β ) is the broadening from contacts (α, β) with Σ α,β being the corresponding energy dependent selfenergy matrices. For computational efficiency, the retarded Green's function G R was obtained using the recursive Green's function algorithm 11 and the self-energies were calculated using the decimation method 17 . The real FIG. 4. Modeling the transition from ballistic to diffusive transport in a micron long graphene, (a-b) sample potential landscape for various impurity concentrations, (c) total conductivity as a function of channel carrier density from ballistic (σ ∝ n) to diffusive (σ ∝ √ n).
space Hamiltonian matrix for the channel and the contacts were obtained using Eq.
(2) and the effect of magnetic field was included using the Peierls substitution 12 .
To calculate the voltage at contact c, we set I a = −I b = I, I c = 0 and µ a = E F then solve Eq. (4) for µ c and µ b where E F is the Fermi level obtained from the electron density n. Then, V = µ c − µ b and resistance R ≡ V /I. The whole procedure was repeated for each n and B. Fig. 3(c) shows the Resistance R as a function of the electron density n and magnetic field B. The color map was generated using the quantum mechanical (NEGF) approach and the the dashed lines were computed using the semi-classical formula Eq. 3. The dashed lines and the bright bands represent the focusing of electrons to contact c. These results agree well with the experiment 9 .
As a final example, we show how our model can interpolate between the ballistic and diffusive limits. We model the transport with impurity scattering by using a sequence of Gaussian potential profiles for the scat-tering centers 18 , U (r) = Nimp n=1 U n exp (−|r − r n | 2 /2ζ 2 ) that specifies the strength of the impurity potential at atomic site r. r n are the positions of the impurity atoms and ζ is the screening length (≈ 3nm for long range scatterers). The amplitudes U n are random numbers following a Gaussian profile 19 , N imp is the impurity concentration. With U added to H (potential landscape shown in Fig. 4a-b), we study the evolution of electron transport in graphene from ballistic to diffusive (for varying N imp . Fig. 4c shows the graphene conductivity at various channel carrier densities for several impurity concentrations. This time we keep the two ends of the graphene device at constant doping to capture the contact induced doping in graphene. This produces electron-hole asymmetry in the ballistic limit due to formation of pn junction near the contact. At high impurity concentration, the contact resistance becomes less dominant compared to the device resistance and the electron-hole asymmetry washes out, similar to what is seen in experiments 20 . Furthermore at the diffusive limit, σ becomes proportional to n for a sample dominated by long range scatterers and can be fitted with a carrier density (n) independent mobility. As we lower the impurity concentration towards the ballistic limit, the graphene conductance (G) becomes proportional to the effective doping E F , G = 4q 2 W E F T /(πhv F ) leading to a sub-linear σ ∝ √ n), where the transmission T is determined by the metal-graphene contact. Such an evolution of electron transport in graphene has been verified in experiments 21 .
In conclusion, we have shown that for massless Dirac fermions, an additional quadratic term in the Dirac Hamiltonian not only circumvents the Fermion doubling problem in a spatial lattice but also has a huge computational advantage over the atomistic tight binding model. In particular, we have shown that the modified Hamiltonian results in an extremely small matrix on a real space square lattice. As a result, the Hamiltonian is orders of magnitude faster than the tight binding Hamiltonian when used in band structure and quantum transport simulations. We applied this Hamiltonian for micron scaled graphene devices to study magneto-transport and electron transport in ballistic and in diffusive limit. Although only graphene is considered here, it is applicable to any other Dirac materials like topological insulators and can be used to calculate the spin current 22 as well.