Ab initio quantum transport in AB-stacked bilayer penta-silicene using atomic orbitals

The current carried by a material subject to an electric field is microscopically inhomogeneous and can be modelled using scattering theory, in which electrons undergo collisions with the microscopic objects they encounter. We herein present a methodology for parameter-free calculations of the current density from first-principles using density functional theory, Wannier functions and scattering matrices. The methodology is used on free-standing AB-stacked bilayer penta-silicene. This new Si allotrope has been proposed to have a higher stability than any of its hexagonal bilayer counterparts. Furthermore, its semiconducting properties make it ideal for use in electronic components. We unveil the role of the pz orbitals in the transport through a three-dimensional quantum wire and present current density streamlines that reveal the locations of the highest charge flow. The present methodology can be expanded to accommodate many electron degrees of freedom, the application of electromagnetic fields and many other physical phenomena involved in device operation.


Introduction
Recent advances in fabrication techniques 1 and inkjet printing of 2D materials 2 have brought the principles of mesoscopic physics into the foreground in order to explain and predict electron device behaviour in a wide range of occasions. In resistive switching memory devices (RRAM), where a conductive lament is formed under the application of a voltage at the electrodes, 3,4 conductance quantization has been observed 5 and the SET operation has been modelled using the Landauer formula for electron tunnelling. 6 'Filaments' have also been studied in relation to the atomic structure of graphene sheets theoretically 7 while the transition from quantum to classical regime with the use of weighted phonon self-energies has been modelled using networks of sites and the Keldysh Green function formalism. 8 These methods can serve as a valuable tool for explaining experimental topography 9 and scanning tunnelling potentiometry. 10 Silicon has historically been the most widely used semiconducting material in the electronics industry both due to its abundance in nature and due to the low-defect interface that it forms with the insulating SiO 2 used in transistors and other electronic components. Hexagonal silicene 11,12 was one of the rst twodimensional sensations with a plenitude of novel properties, however, when compared to its Dirac cone counterpart, graphene, silicene has a much smaller elastic constant 13 and this poses difficulties in its manipulation in the lab.
Other forms of silicon derivatives have been the focus of studies using systematic materials search, revealing novel allotropes with direct band-gaps 14 and even high mobility values. 15 Penta-silicene is a new form of monolayer silicon that has been observed recently. 16,17 Pentagons are rare atomic congurations that are less frequently studied than hexagonal structures mainly due to their difficulties imposed in their fabrication. Using rst principles calculations, bilayer penta-silicene, whose layers are stacked in 90 angle between them has been found to be more stable than the most stable form of bilayer hexagonal silicene. 18 It is made up of pentagonal rings of Si atoms, similar to pentagraphene 19 and other two-dimensional penta-structures. 20 Twisting angles between the layers of few-layer materials has been proved to be an efficient method for tuning their properties. 21,22 Contrary to the bilayer form with no twisting angle, the AB stacking conguration of penta-silicene has semiconducting properties. This makes it ideal for many different congurations of electronic devices, stemming from the recent surge in heterostructure fabrication. 23 In the structures that we examine herein, parts of the system lose their symmetry in favour of stability. Situations such as this, call for the more intuitive chemical description of the system in terms of its valence bonding conguration and this can be very well described with the use of maximally localized Wannier functions (WFs). The latter, acting as alternative representations of Bloch orbitals, have recently gained popularity for examining a plethora of mesoscopic phenomena due to advancements in localization techniques 24,25 but also when used as tools for high-throughput screening of topological materials. 26 Even more interesting for this study, is the possibility to derive model Tight Binding (TB) Hamiltonians in the WF basis, with specic orbitals involved (i.e. excluding higher conductions bands), arbitrary number of neighbor terms and cut-off values to adjacent cells, and therefore increased accuracy and reduced computational cost, suited to multi-scale modelling of 2D semiconductor transistors. 27,28 Therefore, compared to previous work using Kohn-Sham states and mesoscopic transport using Green's functions, 29 the use of Wannier functions provides more computational advantages. Finally, the method of scattering matrices is both fast and stable 30 and can be expanded to include, among others, spin and other degrees of freedom, 31 electromagnetic elds and time-dependent effects. 32 2 Results & discussion Fig. 1a shows the unit cell used in the calculation. It contains two layers of Silicon with a total of 12 atoms (6 at the top and 6 at the bottom layer). 18 Each atom is labelled with a number (Si1-Si12) as shown in Fig. 1b. The computational details are shown in the Methods section. The well-known bandgap-problem of DFT becomes relevant when considering device characteristics, as current transport is mainly performed by the states around it. This can be overcome using hybrid exchange correlation functionals, Hartree-Fock methods or the GW approximation. 33 We have compared band structure results for two different exchange-correlation functionals (Fig. 2). The calculated bandgap is indirect with 0.321 eV for BLYP and 0.107 eV for the PBE functional, while the minimum of the conduction band is located between GZ in the rst case and between GM in the second. The maximum of the valence band in located at M in both cases. The PBE functional was used further in this work, however, we can arrive at similar conclusions using any of the aforementioned methods.
Each layer of the structure has tetragonal symmetry and each sublattice is twisted 90 to the other one as shown in Fig. 1a. Aierken et al. 18 showed that the structure was stable when the electrons on each of the surfaces of its free-standing form lose their symmetry (Fig. 1b). It is natural that when the material is fabricated on various substrates, depending on the method used and substrate material, bonding and hybridization will be adjusted. Cerdá et al. showed that the pentagons that are formed on Ag(110) host both sp 2 and sp 3 bonded Si atoms, in which the latter also bond with Ag substrate atoms 16 For hexagonal silicene nanoribbons, there has been increased interest in their edge states as doping and hydrogenation has shown to change their electronic and magnetic properties. 34,35 Using orbital projected calculations with DFT, the Lowdin charges 36 of the atoms in the structure were derived ( Table 1). All atoms, except for two at each surface (Si4, 6, 10 and 12) were found to have total charge approximately 4, of which roughly 30% resides in s states and 70% in p states. Then, on each of the surfaces, the Si6 and Si12 atoms lose part their charge to the Si4 and Si10 atoms equivalently, whose hybrid orbitals acquire 35% s and 65% p character.
Orbital projected Density of States (DoS) for three representative atoms is plotted next to the band structure in Fig. 2. Most of the atoms have contributions equivalent to that of Si2, except for the two surface atoms that were described previously, for which their p z orbitals show increased contribution around the Fermi level. The bottom of the conduction band is formed by the p z orbitals of Si6 and Si12, while the top of the valence band is formed by the p z orbitals of the Si4 and Si10.
Wannier90 (ref. 37) was used for the Wannierization. In order to derive the effective Hamiltonian from the atomic orbitals, 48 projections were used with atom-centered orbitals, namely, 8 sp 3 , 4 sp 2 and 4 p z , as they were found in the orbital-projected DFT   calculations. In Fig. 3, the plots of 5 WFs are given. The results correctly reect the sp 3 nature of the equivalent atoms and the remaining p z orbitals of the distorted symmetry atoms. Transport calculations on a free-standing quantum wire (Fig. 4) with a scattering region of dimensions 9 Â 2 Â 1 unit cells were performed using Kwant. 30 The length expands in the k x direction. The size of the leads is 3 Â 2 Â 1 with 1D translational symmetry in the directions away from the scattering region. Both the scattering region and the leads are of the same material type, as this allows us to concentrate solely on its properties. A generalization of this can include leads of different material type, or a scattering region with an interface between two different materials. The solution of the Schrodinger equation in the system now corresponds to the plane wave nature of the electrons and not the Bloch waves used in DFT. The Hamiltonian of the system is given by, where i and j denote Wannier sites with V i ¼ hf i |H|f i i the diagonal elements of the Hamiltonian matrix derived aer the Wannierization procedure, c i and c † i denote the electron creation and annihilation operators respectively, and t ¼ hf i |H|f j i the off-diagonal matrix elements equivalently. The on-site energies V i correspond to those of the 48 WFs of the penta-Si atomic orbitals and the hopping integrals t represent the possibility for an electron to jump from state f i to state f j . The whole system is Hermitian. Expanding to further degrees of freedom could be done by representing one site with one atom in the system, and each orbital being an element in the on-site Hamiltonian matrix of the atom. The on-site and hopping integral energies were extracted using TBModels. 26,38 Spin can also be added by solving the spinpolarized Kohn-Sham equations and then separating the up and down spin components from the Wannier90 results. 39 In the case examined here, the calculations are closed-shell, and therefore, we include a factor of two in all further units that include electron charge (e). We also note that all wavefunctions resulting from the atomic orbitals are one-body states.    The on-site energies of the orbitals are given in Fig. 5a. Shown in red are the on-site energies of the p z orbitals of the equivalent atoms, whose energies are higher than the rest of the atomic sp orbitals. Fig. 5b shows the average value of the hopping integrals for each orbital in the unit cell. Both the on-site energies and the hopping integrals show a slight increase in the atoms of the upper layer of the system. The asymmetry of the nal TB Hamiltonian is a problem that arises aer the disentanglement procedure and which has only recently been addressed. 40 The TB Hamiltonian matrix resulting from the Wannierization is cut-off at three nearest neighbour hoppings. In this context, the nearest neighbours are the hoppings that lie within the home unit cell, second nearest neighbours are the hoppings to adjacent cells in all directions etc. The maximum order of the nearest neighbour interactions is dependent on the k-points used in the Monkhorst-pack mesh during the DFT calculations, as they dene the number of periodic images that will appear when the system is translated in real space during Wannierization. As the number of nearest neighbours increases, a large number of neighbour interactions can signicantly increase computation time both for creating the system and solving it. However, in most cases, three nearest neighbour interactions provide sufficient accuracy (see ESI †).
The conductance of the wire is calculated from the scattering matrix of the system. 30 The dispersion relations within the leads that have 1D translational symmetry are shown together with the conductance plots in Fig. 6. The energies of the modes take on constant values in each direction due to the lack of a 3D translational symmetry. The magnitude of the conductance reects the probability of transmission of the modes at each electrochemical potential difference between the leads. A value of two shows the there are two modes propagating in the wire at this energy etc.
Local charge and current densities are accessible through solving for individual or pairs of sites in the system, giving access to computations of many properties that exist at surfaces and interfaces. Fig. 7a shows the averaged local density of states from each orbital type which is taken from the expectation value of the local density operator in Kwant using, Fig. 7 (a) Local density of states for each orbital type in the quantum wire (b) orbital-projected DoS derived from DFT calculations using the PBE functional. Fig. 8 (a) Vector plot showing the location of the hopping with the highest current density for each site. The colour scale shows the magnitude of the expectation value (2e h À1 ). The arrowhead shows the direction of the current at this hopping (b) interpolated current density and streamlines at two cuts perpendicular to k z at the vicinity of the highest concentration of local current density. Transport is in the k x direction (lead 1 towards lead 2). The locations of the atoms that exist at this z cut (with an approximate displacement of 1Å) are marked in blue for the first four unit cells. Both results are taken at where j is any propagating wavefunction in the scattering region at each energy, n runs over all expectation values resulting from all incoming and outgoing wavefunctions and k runs all sites representing each orbital type in the system. At E f À0.19 eV, an increased electron density is observed at Si7 sp 3 and Si12 p z orbitals. To explain this, the expectation value of the current density at this energy was also extracted for the propagating wavefunctions originating from lead 1 (le lead in Fig. 4). For a hopping from site l to site k, this is calculated from, X where n runs over all aforementioned propagating wavefunctions. One of the many advantages of using this method is that it gives us access to three-dimensional quantities. To visualize the results we have used both a vector plot and a 2D cut of the scattering region. Fig. 8a shows the magnitude and direction of the highest hop for each site in the system. Global maxima occur between Si7 sp 3 / Si12 sp 2 and Si8 sp 3 / Si12 p z . For a more clear view, Fig. 8b shows a 2D cut at the location of the highest current density in the system and the atoms that are located at its vicinity for the rst four unit cells. It is seen more clearly that transport occurs in the right most side of the wire.
In the absence of any potential or magnetic eld dened explicitly over the 3D region, the distribution of the current density is dictated by the connement effects in the x, y and z directions, 7 that is by the solution of the Schrodinger equation with the translation operator in the leads dening the energies of the propagating modes. Finally, the on-site energies as well as the hopping integral values dene the electrostatic landscape for the electrons to nd their route through the wire (Fig. 5). ‡ Taking into account collective effects from many orbitals, the effects of impurities and screening or the formation of dipoles at interfaces 41 can furthermore be examined, opening the road to versatile computational microscopic and topographic studies, while phonon effects can be added as site self-energies using standard procedures. 8 Progress in obtaining and engineering model TB Hamiltonians from rst principles is currently ourishing 42 and is expected to lead to signicant advances in predictions of device characteristics.

Conclusions
An ab initio multi-scale simulation approach has been presented for the calculation of the current density between adjacent atomic positions in a material using an effective Hamiltonian derived from Wannier functions. The methodology has been applied to a newly predicted material by the name bilayer pentasilicene, where we have observed an increased concentration of charge at specic orbitals in a free-standing quantum wire. This was found to be consistent with the expectation values of the DC local current between its atomic orbitals, which revealed the locations of the highest ow of charge.
This methodology presents many advantages for the examination of electron device operation as it allows many microscopic details of the quantum transport to be revealed using realistic values and the addition of disorder, spin, phonon and nite-temperature effects. Bias calculations are also possible by summing the propagating modes around the Fermi level in the scattering region. 43

DFT
Density functional theory calculations were performed using the Quantum Espresso package, 44 norm-conserving Goedecker/ Hartwigsen/Hutter/Teter pseudopotential with BLYP § and PBE{ exchange correlation functionals. 96 bands were included in the calculation for the 12 atoms of the unit cell, each with 4 valence electrons (3s 2 3p 2 ). The vacuum between periodic images in the z direction was set to 22Å. The plane wave cut-off energy was set to 37 Ry and a Monkhorst-pack k-point mesh of 9 Â 9 Â 1 was used for the relaxation and the band structure calculations. These settings resulted in a lattice constant of 5.29Å for the case of the BLYP and 5.21Å for PBE.

Wannier calculations
Wannier90 (ref. 37) calculations were performed with PBE functional DFT results given as input. A frozen window that included the 26 valence and 22 conduction bands was dened. The tolerance for the gauge invariant term of the spread of the WF was set to the really low value of 10 À7Å2 as the bands towards higher energies are highly entangled. The density of the k-point mesh was explicitly optimized to 15 Â 15 Â 1 k-points. Imaginary/real ratios of the order 10 À3 to 10 À4 were achieved for the WFs.

Conflicts of interest
There are no conicts to declare. through the National Grid Infrastructures NGI GRNET, Hellas-GRID as part of the SEE Virtual Organisation. Data presented in this paper are available from DOI: 10.5281/zenodo.1438840.