Quantum transport regimes in quartic dispersion materials with Anderson disorder

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 ¯m2$ 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,⁶⁰ and SnSb⁶¹ 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.

The rest of this paper is organized as follows. In Sec. II, the computational models of interest are introduced. The effects of disorder on conductance are separately discussed in Sec. III for both 1D and Q1D systems. Section IV summarizes our conclusions.

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$ = $∑ i ϵ i c i † c i$ to the Hamiltonian. Here, $ϵ i$ represents the on-site energy, which randomly fluctuates in the energy interval [ $−W$/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, $∑ i ϵ i=0$⁠. In all cases, Eq. (1) with additional on-site terms is numerically solved for various system lengths, $L$⁠.

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.⁶⁹ As for the statistical average of conductance, simple $⟨G⟩$ 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 $⟨ln⁡G⟩$ 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⁡⟨ln⁡G⟩$⁠, which is expected to give a better statistical result.^72,73

For investigating the effects of disorder, we restrict ourselves to the energies around the VB edge. In the chosen energy interval $(E− E VBM)/ t 1∈[−0.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× 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 $ℓ mfp$ and $ℓ 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 $ℓ 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), $ℓ mfp$ and $ℓ loc$ of quadratic chain are equal, in agreement with the Thouless relation, i.e., $ℓ loc/ ℓ mfp=1$ for $N ch=1$⁠. In a similar fashion, $ℓ mfp$ and $ℓ loc$ curves for the quartic chain are displayed in Fig. 3(c). As expected, $ℓ mfp$ of the quartic chain equals to $ℓ loc$ except for the quartic edge with $N ch=2$⁠, where the ratio turns out to be $ℓ loc$/ $ℓ mfp$ $∼$ 1.5 as suggested by the Thouless relation. Interestingly, $ℓ mfp$ (⁠ $ℓ 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 $ℓ mfp=700a$ and $ℓ loc=1010a$ at $E= E VBM−0.05 t 1$ [see the inset of Fig. 3(c)]. In comparison, we compute transmission at $E= E VBM− t 1$ with a much stronger disorder (⁠ $W=250$ meV) and find $ℓ mfp=1983a$ and $ℓ loc=1996a$⁠, which are several times larger than those in the MHS energy region with weaker disorder. $ℓ mfp$ and $ℓ 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.

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, $ρ(E)$ per area, and $G ballistic$ as a function of dimensionless energy $(E− 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 $ρ$ 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 $ρ(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 $ρ(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− E VBM)/ t 1∈[−0.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 $ℓ 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 $ℓ 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, $ℓ 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, $ℓ mfp$ converges to zero as $(E− E VBM)/ t 1→0$⁠. On the contrast, $ℓ mfp≈3200a$ at energies outside the MHS energy region, namely, for $(E− E VBM)/ t 1∼−0.1$ for ANR.

The $ℓ loc$ values are obtained by fitting the $G av$ at $L≫ ℓ mfp$⁠. The consistency of $ℓ mfp$ and $ℓ 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 $ℓ loc$ for the ZNR (red solid curve) and the ANR (blue solid curve) in which $ℓ 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 $ℓ mfp$ independent of the edge termination.

In quasi-1D structures, the system width affects $ℓ 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 $ℓ mfp$ is expected to reduce with ribbon width. Indeed, it was reported that $ℓ mfp$ increases with width in edge disordered graphene nanoribbons.^76,77 Similarly, in graphene NRs with oxygen functionalization, $ℓ mfp$ was shown to increase with the ribbon width.⁷⁸ However, in quartic NRs, we do not observe such an increase in $ℓ mfp$ with width. On the contrary, $ℓ mfp$ can be greater for narrower ribbons at some energy values. This can be understood in terms of another factor that determines $ℓ 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 $ℓ 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, $ℓ 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 $ℓ 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 $ℓ mfp$⁠, $ℓ loc$ increases with width for both ANRs and ZNRs [Figs. 7(c)–7(f)]. Comparing $ℓ loc$ of ANRs and ZNRs, those of ZNRs are considerably longer, whereas their $ℓ 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 $ℓ loc$ of quartic NRs can be understood by investigating the details in the DOS and $G ballistic$⁠.

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

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

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 $ℓ mfp$ and $ℓ loc$ of the quartic chain are much smaller in the MHS energy region. In quasi-1D quartic ribbons, $ℓ mfp$ and $ℓ loc$ are considerably short because of multiple strong Van Hove singularities, which are denser and stronger in quartic systems compared to quadratic ones. Interestingly, $ℓ 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.

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.

The authors have no conflicts to disclose.

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).

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