Mexican-hat-shaped quartic dispersion manifests itself in certain families of single-layer two-dimensional hexagonal crystals such as compounds of groups III–VI and groups IV–V as well as elemental crystals of group V. A quartic band forms the valence band edge in various of these structures, and some of the experimentally confirmed structures are GaS, GaSe, InSe, SnSb, and blue phosphorene. Here, we numerically investigate strictly one-dimensional and quasi-one dimensional (Q1D) systems with quartic dispersion and systematically study the effects of Anderson disorder on their transport properties with the help of a minimal tight-binding model and Landauer formalism. We compare the analytical expression for the scaling function with simulation data to distinguish the domains of diffusion and localization regimes. In one dimension, it is shown that conductance drops dramatically at the quartic band edge compared to the quadratic case. As for the Q1D nanoribbons, a set of singularities emerge close to the band edge, suppressing conductance and leading to short mean-free-paths and localization lengths. Interestingly, wider nanoribbons can have shorter mean-free-paths because of denser singularities. However, the localization lengths sometimes follow different trends. Our results display the peculiar effects of quartic dispersion on transport in disordered systems.

## I. INTRODUCTION

Several families of two-dimensional (2D) materials have been identified since the first isolation of a graphene monolayer, and they have received considerable attention due to their unusual and unique properties. They are not only interesting for fundamental science but are also promising candidates for next-generation devices.^{1–12} Many of the novel properties are related to their electronic band structures, some of which are unprecedented. Mexican hat-shaped (MHS) quartic energy dispersion is one of them. Graphene like honeycomb lattices of group V elements and group III–VI and group IV–V compounds in hexagonal $P 6 \xafm2$ symmetry display MHS quartic dispersion in their valance bands (VBs). In some of these structures, the valence band maximum is formed by the MHS quartic band.^{13–46} Owing to their unusual band structure, these materials are reported to have peculiar magnetic and ferroelectric phases, very high thermoelectric efficiencies, and the potential to be useful for applications like water splitting, optoelectronics, photonics, and thermoelectrics.^{21,30–36,45,47–51} On the experimental side, single and few layers of quartic materials such as gallium based compounds (GaS, GaSe, GaTe),^{30,52–56} InSe,^{57–59} blue phosphorene,^{60} and SnSb^{61} have been synthesized and their distinctive properties have been reported. It is worth noting that similar band dispersion has been studied in the context of optical lattices, with implications about chiral spin liquid, and that the tight binding approach has been used to show the absence of Bose condensation on these lattices.^{62–64}

A key feature in the electronic structure of MHS quartic dispersion is the strong (inverse-square-root) Van Hove singularity with divergent DOS at the VB edge. Such strong singularity and divergent DOS is uncommon in two dimensions. In addition to the strong singularity, the pristine transmission spectrum displays a step-like behavior at the band edge. A stepwise transmission spectrum is a characteristic feature of 1D systems, and in two dimensions, it is observed only in quartic materials. These peculiarities of quartic dispersion, together with finite amounts of disorder, are potentially responsible for strong scatterings that could alter the electronic properties significantly. The effects of disorder on the transport properties of quartic materials have not been addressed in the literature before. Here, using nonequilibrium Green’s function (NEGF) methodology, we study quantum transport properties around the quartic VB edge. We consider large-scale quartic systems by taking advantage of a minimal tight-binding (TB) model,^{21,62,63} and thus, we are able to determine the transport length scales of these systems in the presence of the Anderson disorder. The effects of the Anderson disorder on conductance ( $G$) are numerically investigated in the case of strictly one-dimensional (1D) and quasi-one dimensional (Q1D) systems at the weak disorder limit. Our results show that there is a dramatic suppression of transmission at the quartic band edge, which cannot be understood using the conventional effective mass approaches or any multiple band models.

## II. SYSTEMS AND METHODS

^{21}

^{65}which can be expressed as $ E \xb1( k)=\xb1 t 1 3 + f ( k )\u2212 t 2f( k)$, with $f( k)=2cos( 3 k ya)+4cos( 3 k ya/2)cos(3 k xa/2),$ $ k$ and $a$ being the wave vector and the lattice constant, respectively. Different from the half-filled bands of $\pi $-orbitals of graphene, we are interested in completely filled bands (e.g., replacing carbon with nitrogen, see Ref. 21). In this sense, valance band maximum (VBM) is set to that of the highest occupied electronic state by assuming two electrons per atom in this study. For $\xi $ = $ t 2/ t 1>1/6$, the above dispersion yields a quartic band around the VB edge, such that

^{16,21,27}That is, $ t 1=6.1$ eV and $ t 2=1.27$ eV are chosen, which corresponds to $\xi =0.21$. In Fig. 1(a), the quartic band for nitrogene as obtained from our TB approximation is shown to have a good agreement with that obtained from density functional theory (DFT) calculations, especially for low energy holes.

^{21}For the sake of simplicity, we refer the energy values ranging from the minimum of the Mexican hat at the center of the BZ to its maximum, namely, the VBM, as the MHS energy region. This is not meant to exclude the energies close to these values from being quartic in dispersion, but to emphasize that in this energy range, there are four solutions at a given direction in $k$-space for the 2D structure. In the case of hexagonal nitrogene, the MHS energy region lies within $\u22120.47$ eV $<E\u2212 E VBM<0$.

In this study, starting from a toy model, namely, monatomic chain, we investigate the electronic and transport properties of the Q1D nitrogene nanoribbons (NRs) with zigzag and armchair edges (ZNR and ANR) in the presence of uncorrelated disorders. The hexagonal lattice structure and the unit cells of these Q1D systems are illustrated in Fig. 1(b). The puckered geometry of the lattice does not play a role for the purposes of this study, therefore, it is disregarded. As for short-range disorder, the Anderson disorder is introduced by adding the term $ H A$ = $ \u2211 i \u03f5 i c i \u2020 c i$ to the Hamiltonian. Here, $ \u03f5 i$ represents the on-site energy, which randomly fluctuates in the energy interval [ $\u2212W$/2, $W$/2] with a fixed disorder strength $W=25$ and $250$ meV for the chain and ribbon geometries, respectively. The overall contribution of the on-site potential energies is set to zero, that is, $ \u2211 i \u03f5 i=0$. In all cases, Eq. (1) with additional on-site terms is numerically solved for various system lengths, $L$.

^{66}Scattering events are allowed to take place only in the $C$ region for which the retarded Green function can be expressed as $ G C C r(E)= [ ( E + i 0 + ) I \u2212 H C C \u2212 \Sigma L r \u2212 \Sigma R r ] \u2212 1$, with $ 0 +$ being an infinitesimal positive number, $I$ the identity matrix, $ H C C$ is Hamiltonian of the $C$ region, $ \Sigma L / R r= H C L / C R G L L / R R r , 0 H L C / R C$ are the self-energy matrices with the free Green’s functions $ G L L / R R r , 0= [ ( E + i 0 + ) I \u2212 H L L / R R ] \u2212 1$ of the isolated $L$ and $R$ reservoirs. The level-broadening caused by system-electrode coupling is $ \Gamma L / R=i ( \Sigma L / R \u2212 \Sigma L / R \u2020 )$, and the transmission amplitude is given by

^{67,68}

System sizes should be large enough for a reliable analysis of the localization regime. In order to be able to simulate large systems a recursion–decimation algorithm is implemented.^{69} As for the statistical average of conductance, simple $\u27e8G\u27e9$ averaging over a set of samples may not converge toward any meaningful value since $G$ does not follow Gaussian statistics in the localized regime.^{70–73} Therefore, a logarithmic average $\u27e8ln\u2061G\u27e9$ over 100 samples at each system length $L$ is considered in this study. In this way, statistical average of conductance is obtained via $ G av=exp\u2061\u27e8ln\u2061G\u27e9$, which is expected to give a better statistical result.^{72,73}

^{74}Interrelation between transport length scales in 1D systems is set by the Thouless relation based on the random matrix theory, which conditions that $ \u2113 loc$ = ( $\eta $[ $ N ch$-1]+2)/2 $ \u2113 mfp$,

^{73}where $\eta $ = 1 in the absence of external magnetic field, simplifying to

## III. RESULTS

### A. Monatomic chain with quartic dispersion

For investigating the effects of disorder, we restrict ourselves to the energies around the VB edge. In the chosen energy interval $(E\u2212 E VBM)/ t 1\u2208[\u22120.25,0]$, $ G av$ is plotted as a function of $E$ and $L$ in Fig. 3(a) for the quadratic (top panel) and quartic (bottom panel) chains. $ G av$ is computed for lengths ranging from 0 up to $5\xd7 10 4a$. At short distances, $ G av$ decreases very fast at the quartic edge and the step exhibits a rounded shape. At energies away from the VBM, conductance decreases relatively slowly, in fact slower than the corresponding energies of the quadratic band. Simulation data of the quadratic and quartic chains are fitted with Eq. (4), and $ \u2113 mfp$ and $ \u2113 loc$ are shown in Figs. 3(b) and 3(c), respectively. It should be noted that the simulation data above $ G ballistic/2$ are used for the $ \u2113 mfp$ fitting processes of both chains, which will be discussed in further detail in Sec. III D. As it is shown in Fig. 3(b), $ \u2113 mfp$ and $ \u2113 loc$ of quadratic chain are equal, in agreement with the Thouless relation, i.e., $ \u2113 loc/ \u2113 mfp=1$ for $ N ch=1$. In a similar fashion, $ \u2113 mfp$ and $ \u2113 loc$ curves for the quartic chain are displayed in Fig. 3(c). As expected, $ \u2113 mfp$ of the quartic chain equals to $ \u2113 loc$ except for the quartic edge with $ N ch=2$, where the ratio turns out to be $ \u2113 loc$/ $ \u2113 mfp$ $\u223c$ 1.5 as suggested by the Thouless relation. Interestingly, $ \u2113 mfp$ ( $ \u2113 loc$) of the quartic chain is approximately six (four) times smaller than that of the quadratic chain at the quartic edge [see insets of Figs. 3(b) and 3(c)]. In the MHS energy region with a disorder strength of $W=25$ meV, we have $ \u2113 mfp=700a$ and $ \u2113 loc=1010a$ at $E= E VBM\u22120.05 t 1$ [see the inset of Fig. 3(c)]. In comparison, we compute transmission at $E= E VBM\u2212 t 1$ with a much stronger disorder ( $W=250$ meV) and find $ \u2113 mfp=1983a$ and $ \u2113 loc=1996a$, which are several times larger than those in the MHS energy region with weaker disorder. $ \u2113 mfp$ and $ \u2113 loc$ are almost equal because $ N ch=1$ at these energies. This comparison in another illustration of the role of strong singularity in transport properties of quartic systems.

^{75}Assuming that the Bloch wave has equal probability weight on each atom in the unit cell, one can write

### B. Nanoribbons with quartic dispersion

The electronic and transport properties of the Q1D quartic systems are investigated for ZNR and ANR of nitrogene with the TB parameters of $ t 1=6.1$ eV and $ t 2=1.27$ eV. Widths are chosen as $N=10$, i.e., $ N atom=20$ in a unit cell for both cases. The full spectra of considered bands, $\rho (E)$ per area, and $ G ballistic$ as a function of dimensionless energy $(E\u2212 E VBM)/ t 1$ are shown in the upper panels of Fig. 4. Valence band edges of the corresponding spectra are displayed in the lower panels. In the band structure of ZNR, a family of quartic bands appears around the VBM, where some of them exhibit crossings (see the bottom-left of the left panel of Fig. 4). In two dimensions, the critical wave-vectors that form the band edge are degenerate and nonisolated, which give rise to a strong Van Hove singularity. In Q1D NRs, the critical wave-vectors are isolated and nondegenerate, whereas Van Hove singularities are still strong because of the reduced dimension. In addition, there exists numerous strong singularities in the entire spectrum, which are characteristic in Q1D structures. A distinguishing character of $\rho $ in Q1D quartic materials from those of other materials is the emergence of excessively dense singularities at the MHS energy region (see Fig. 4). Both ANR and ZNR display rapid changes in their $ G ballistic$ close to the VBM (Fig. 4). In ZNR, top bands are more dispersive than in ANR, which is because of the narrower Brillouin zone of ANR. As a result, $ G ballistic$ values are larger for ZNR in these energies. Interestingly, we observe formation of nearly flatbands in ZNRs, close to the band edge. These nearly flatbands are formed by hybridization of quartic bands due to the edges. Flattening gives rise to even stronger peaks in the DOS.

It can be expected that the increase in the NR’s width causes both the DOS and $ G ballistic$ spectra to evolve into those of the 2D system. Figure 5 exhibits the evolution of $\rho (E)$ per area and $ G ballistic$ for different system widths, $N$, in the case of pristine ZNRs. Here, the corresponding 2D spectra are given by black curves. It can be observed in Fig. 5-top that $\rho (E)$ approaches to that of the 2D quartic system as $N$ is systematically increased. Similarly, $ G ballistic$ of NRs approach to that of 2D quartic systems with increasing $N$, which is shown in the bottom panel of Fig. 5.

The edge shape is a determinant factor for the electronic structure and hence for the transport regimes. The effects of short-range Anderson disorder are studied for the Q1D systems for which the numerical calculations are performed for the ribbon widths of $N=10$ for both edge shapes. The evolution of conductance is examined for various device lengths at energies close to the VBM ( $(E\u2212 E VBM)/ t 1\u2208[\u22120.12,0]$) for a relatively weak realization of disorder ( $W=250$ meV). Length dependent $ G av$ values of ZNR and ANR are plotted in Fig. 6(a). Although $ G ballistic$ values are close within this range, length dependent $ G av$ decays much faster for the MHS energy region. This behavior is the same for both ANR and ZNR. The origin of this behavior can be understood by examining the pristine DOS at these energies. As it is shown in Figs. 6(b) and 6(c), the pristine DOS (red curves) has multiple singularities at the energies, where $ G av$ drops suddenly. The computed $ \u2113 mfp$ values are shown in the same plots (blue curves). We note that simulation data around $ G ballistic/2$ are used for fitting the diffusion formula, cf. Eq. (4). As it is expected, the dips in $ \u2113 mfp$ correspond to the peaks in the DOS. It is much probable for the particle to scatter due to the higher DOS. In the MHS energy region, where the DOS values are high, $ \u2113 mfp$ of ZNR (ANR) always remains below 300 $a$ (400 $a$). In the first conduction step, it is lower than $10a$ and at the band edge, $ \u2113 mfp$ converges to zero as $(E\u2212 E VBM)/ t 1\u21920$. On the contrast, $ \u2113 mfp\u22483200a$ at energies outside the MHS energy region, namely, for $(E\u2212 E VBM)/ t 1\u223c\u22120.1$ for ANR.

The $ \u2113 loc$ values are obtained by fitting the $ G av$ at $L\u226b \u2113 mfp$. The consistency of $ \u2113 mfp$ and $ \u2113 loc$ is checked by using the Thouless relation [Eq. (5)]. Within the MHS energy region, the localization lengths are short, and, therefore, they can be obtained within reasonable sizes of simulated devices. However, for lower energies, this is not the case. Figure 6(d) demonstrates $ \u2113 loc$ for the ZNR (red solid curve) and the ANR (blue solid curve) in which $ \u2113 loc$ approaches zero for both NRs as the energy goes to $ E VBM$. It is clear from Fig. 6(d) that for both NRs, Thouless relations (black dashed curves) are in good agreement with the fitted data, especially within the MHS energy region.

In order to reveal the effects of device width $N$ on the transport length scales, the Q1D systems with $N=10$ are compared to those with $N=20$ for both edge shapes. Figure 7(a) shows pristine DOS per area for the ZNRs with $N=10$ (red solid curve) and $N=20$ (blue solid curve) at the MHS energy region. The spectra explicitly display that $N=20$ case has more singularities compared to the ZNR with $N=10$ at the band edge. Multiple singularities exist around $ E VBM$ for both sizes. Similarly, DOS for the ANR with $N=10$ (pink solid curve) and $N=20$ (orange solid curve) are displayed in Fig. 7(d). It is observed that the density of singularities increase with width. On the other hand, when the number of atoms in the unit cells are equal, the number of singularities in ZNRs is larger than in ANRs simply because the number of bands within the MHS energy region is larger in ZNRs (see Fig. 4). Dense singularities in the DOS give rise to shorter $ \u2113 mfp$ independent of the edge termination.

In quasi-1D structures, the system width affects $ \u2113 mfp$ by means of more than one mechanisms. First, considering the ribbon as a dimensionally reduced form of the two-dimensional structure, edges are sources of scatterings, hence $ \u2113 mfp$ is expected to reduce with ribbon width. Indeed, it was reported that $ \u2113 mfp$ increases with width in edge disordered graphene nanoribbons.^{76,77} Similarly, in graphene NRs with oxygen functionalization, $ \u2113 mfp$ was shown to increase with the ribbon width.^{78} However, in quartic NRs, we do not observe such an increase in $ \u2113 mfp$ with width. On the contrary, $ \u2113 mfp$ can be greater for narrower ribbons at some energy values. This can be understood in terms of another factor that determines $ \u2113 mfp$, namely, the width dependence of DOS at the MHS energy region. Recalling Fermi’s golden rule, the scattering rate increases with the DOS, which is the reason of the dips in the computed $ \u2113 mfp$ [see Figs. 6(b), 6(c), and 7(b)–7(e)]. In quartic NRs, the number of Van Hove singularities depends on the ribbon width, namely, narrower the ribbon, less the number of singularities. As a result, $ \u2113 mfp$ of wider ribbons are suppressed due to denser singularities. The trade-off between these two mechanisms is the main factor that determines width dependence of $ \u2113 mfp$.

The localization length also depends on these factors. Additionally, the number of channels has an important role as it was described in Eq. (5). The number of channels increases with width, therefore, unlike $ \u2113 mfp$, $ \u2113 loc$ increases with width for both ANRs and ZNRs [Figs. 7(c)–7(f)]. Comparing $ \u2113 loc$ of ANRs and ZNRs, those of ZNRs are considerably longer, whereas their $ \u2113 mfp$ are similar. This difference is also rooted in the difference between the number of channels. ANRs have wider unit cells, and narrower Brillouin zones; therefore, the quartic bands are less dispersive and the number of channels are less than ZNRs. Further details in the energy dependent $ \u2113 loc$ of quartic NRs can be understood by investigating the details in the DOS and $ G ballistic$.

### C. 2D limit

^{75}As tube radius increases, DOS converges to that of the 2D structure. In the MHS region, it is given as

^{79}

^{21}and $ \u2113 mfp$ increases linearly with $|E|$. At the same disorder strength as those in quasi-1D structures, $ \u2113 mfp$ of the 2D quartic system is plotted in Figs. 7(b)–7(e) for comparison. Owing to the high DOS near $ E VBM$, $ \u2113 mfp$ of the 2D system mostly stays below those of NRs and goes to zero at the quartic band edge.

^{80}

^{81}

### D. Scaling analysis

^{82–84}is used. In the case of disordered systems, the metal-insulator transition (MIT) can be studied through $\beta $ as a function of ln( $g$), where the dimensionless conductance is $g= G av/ G o$. There is a fixed point above (below) which the sytem scales toward a metal (an insulator) as the system size is increased. However, the system dimension plays a crucial role in our conceptual understanding of MIT. The 2D system is the marginal case,

^{83}but the systems with a dimension smaller than $d<2$ undergo to the Anderson localization.

^{82}In addition, the Q1D quartic systems have the multiple strong singularities at the quartic band edge, pointing to the strong localization regime. For these reasons, we make use of the $\beta $-function to reveal the points where the diffusion regime ends and the localization regime starts. The two different forms of the $\beta $-function given below render easy to distinguish the transport regimes.

^{85}

In Figs. 8(a) and 8(b), we compare $\beta $ with $ \beta data$ for the strictly 1D systems, namely the quadratic and quartic chains, as functions of $L$ for $(E\u2212 E VBM)/ t 1=\u22120.017$. The solid curves represent the $\beta $ function ( $ \beta dif$ and $ \beta loc$), which are obtained using $ \u2113 mfp$ and $ \u2113 loc$ [which were obtained by fitting simulation data to Eq. (4)]. The data points represent $ \beta data$ ( $ \beta dif data$ and $ \beta loc data$), which do not rely on any fittings but are obtained solely from the simulation data. At short lengths, there is a perfect agreement between $ \beta dif$ and $ \beta dif data$ [Figs. 8(a)–8(b), red curves and red dots]. This also verifies that the $ \u2113 mfp$ fit is satisfactory. The $ \u2113 mfp$ values are indicated by vertical dashed curves ( $ \u2113 mfp=1903a$ for quadratic chain and $ \u2113 mfp=308a$ for quartic chain). For $L> \u2113 mfp$, the agreement between $ \beta dif$ and $ \beta dif data$ is lost, namely, the simulation data do not obey the diffusion equation any more. This suggests that the fitting procedure for $ \u2113 mfp$ should use only the data from short distances (where $ G av\u2243 G ballistic/2$), as was done in this work. At long distances ( $L\u226b \u2113 mfp$), we observe that the $ \beta loc data$ converges to $ \beta loc$ (blue curves and dots), which confirms that the fitted $ \u2113 loc$ is correct. We should note that the $ \beta loc data$ converges to $ \beta loc$ at distances longer than $ \u2113 loc$. The intermediate distances where $ \beta data$ is not in agreement with $ \beta dif$ or $ \beta loc$ are considered the crossover region. The disagreement does not imply that $ \u2113 mfp$ or $ \u2113 loc$ are incorrect. The validity of $ \u2113 loc$ is confirmed by the fact that the slope of $ \beta loc$ matches perfectly with $ \beta loc data$ at long distances.

The same scaling analysis is performed for quartic NRs as well. In Figs. 8(c) and 8(d), $\beta $ and $ \beta data$ are shown for ZNR( $N=10$) and ANR( $N=10$) at $(E\u2212 E VBM)/ t 1=\u22120.0023$, respectively. Transport length scales are found to be $ \u2113 mfp=12a$ ( $28.8a$) and $ \u2113 loc=842.26a$ ( $51.43a$) for the ZNR (ANR). As it was the case in strictly 1D systems, $ \beta dif data$ is in very good agreement with $ \beta dif$ at short $L$ up to $ \u2113 mfp$. At $L\u226a \u2113 mfp$, $ \beta loc data$ converges to $ \beta loc$, showing that obtained $ \u2113 loc$ is consistent. We again observe that diffusion regime survives only at relatively short distances in quartic NRs. Although $ N ch$ can be large, $ \u2113 loc$ is also relatively short because of multiple strong singularities close to the VBM. These findings are confirmed by the $\beta $-function analysis of simulation data as well.

## IV. CONCLUSION

In summary, a minimal TB model is combined with Landauer formalism, which is used for studying the transport properties of disordered quartic systems. The scaling theory is utilized in a way to distinguish the diffusion and localization regimes and to find the relevant length scales. A comparison of a strictly 1D quartic chain against its quadratic counterpart at a given disorder strength shows that $ \u2113 mfp$ and $ \u2113 loc$ of the quartic chain are much smaller in the MHS energy region. In quasi-1D quartic ribbons, $ \u2113 mfp$ and $ \u2113 loc$ are considerably short because of multiple strong Van Hove singularities, which are denser and stronger in quartic systems compared to quadratic ones. Interestingly, $ \u2113 mfp$ can be longer in narrower NRs compared to the wider NRs, which is because the number of singularities increases with width in the MHS energy region. Exploration of the 2D limit suggests very long localization lengths and requires further numerical studies. Transport length scales of randomly distributed defects can also be studied using the same methodology. Dependence of conductivity on the sample quality and the doping rate could be interpreted using the results presented in this study.

## ACKNOWLEDGMENTS

This work was supported by The Scientific and Technological Research Council of Turkey (TÜBİTAK) under 1001 Grant Project No. 119F353. Support from the Air Force Office of Scientific Research (AFOSR, Award No. FA9550-21-1-0261) is also acknowledged.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Mustafa Polat:** Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Hazan Özkan:** Investigation (supporting); Validation (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). **Hâldun Sevinçli:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

## REFERENCES

*et al.*,

*et al.*,

*et al.*,

*Electronic Transport in Mesoscopic Systems*

*Theory of Quantum Transport at Nanoscale: An Introduction*