Application of the WKB Theory to Investigate Electron Tunneling in Kek-Y Graphene

Abstract

In this paper, we have constructed a WKB approximation for graphene having a Y-shaped Kekulé lattice distortion and a special folding of the K and ${\mathit{K}}^{\prime }$ valleys, which leads to very specific linear energy dispersions with two non-equivalent pairs of subbands. These obtained semi-classical results, which include the action, electron momentum and wave functions, are utilized to analyze the dynamics of electron tunneling through non-square potential barriers. In particular, we explore resonant scattering of an electron by a potential barrier built on Kekulé-distorted graphene. Mathematically, a group of consecutive equations for a semi-classical action have been solved by following a perturbation approach under the condition of small strain-induced coupling parameter ${\Delta }_0\ll 1$ (a good fit to its actual value ${\Delta }_0$∽ 0.1). Specifically, we consider a generalized model for Kek-Y graphene with two arbitrary Fermi velocities. The dependence of the electron transmission amplitude on the potential profile $V\left(x\right)$ and band parameters of Kekulé-patterned graphene has been explored and analyzed in detail.

1. Introduction

Kekulé-patterned graphene represents an unusual and technologically very promising modification of conventional graphene energy band structure due to coupling of its orbital with spin degrees of freedom. In fact, a periodic Y-shape modification of the bond strengths of graphene, as illustrated in Figure 1, leads to chiral electronic states as well as two nonequivalent coupled Dirac cones with different Fermi velocities [1,2,3]. Practically, such a unique feature can be realized by depositing graphene on a special copper-based substrate [4]. Here, we focus on Kek-Y graphene, which is further distinguished from Kek-O graphene by a finite bandgap [4,5,6].

The most crucial feature of Kek-Y graphene is its valley-momentum locking, i.e., two eigenstate modes of electrons are chiral (helical) and have zero bandgap, while both sublattice pseudospin and isospin are locked to the same group velocity [1]. Recently, there have been a number of crucial research publications aimed at addressing these unusual electronic, optical and transport properties of Kekulé-distorted graphene [3,7,8,9,10]. Even though the mismatch between slopes of two Dirac cones appears fixed in Kekul e-patterned graphene, i.e., $v_{F,1}/v_{F,2}\backsim 0.9$, we still consider a general case for Kek-Y graphene with two nonequivalent Dirac cones and having two arbitrary Fermi velocities $v_{F,1}$ and $v_{F,2}$.

Kekulé-patterned graphene is the newest member in the family of Dirac-cone materials, including graphene [11,12], silicene [13,14,15] and transition-metal dichalcogenides [16,17], as well as the most unusual group of $\alpha$-${\mathcal{T}}_3$ materials [18,19,20]. The most uncommon feature of the $\alpha$-${\mathcal{T}}_3$ model is the presence of a flat band in its energy dispersions, giving rise to the most unusual electronic and optical properties [21,22,23,24,25,26,27,28,29,30,31,32,33,34]. On the other hand, the magnetic properties and magnetic-quantization behavior of $\alpha$-${\mathcal{T}}_3$ materials have also attracted a lot of attention [35,36,37,38,39], including a phase transition from a diamagnetic to a paramagnetic susceptibility by continuously varying the parameter $\alpha$ [40]. In fact, the low-energy dispersion of the $\alpha$-${\mathcal{T}}_3$ model looks somewhat similar to that of the generalized Kek-Y graphene model if the Fermi velocity of the lower Dirac cone is greatly reduced compared to the larger one (i.e., $v_{F,2}\ll v_{F,1}$).

Apart from Kek-Y graphing, there is another famous modification of Kekule-patterned graphene, which is referred to as Kek-O. In the latter case, Kekule dimerization takes place in benzene rings of hexagons in graphene.

The energy spectrum of Kek-O graphene demonstrates a large bandgap. Therefore, the case of Kek-O graphene could be considered similar to regular graphene in the presence of circular polarized irradiation discussed above. A substantial bandgap in the energy dispersions of Kek-O graphene evidently results in a tremendous loss of electron transmission over square barrier, as well as significant drop of its optical and Boltzmann conductivity.

The Wentzel–Kramers–Brillouin (WKB), or the semi-classical, theory provides an intuitive but realistic description for dynamical behaviors of electrons in a variety of complex quantum systems under the condition of a high-kinetic energy. In this case, the properties of such electrons are not much different from their classical counterparts and could be effectively described by classical dynamics. However, in contrast to a perturbation theory, the interaction between a particle and an external potential can be very strong, implying that in a number of cases, the WKB description can still be valid even as other types of approximation fail [41,42,43,44,45,46]. Therefore, from this point of view, it is essential and crucial to develop a semi-classical WKB theory for each newly discovered Dirac lattice, such as $\alpha$-${\mathcal{T}}_3$ or Kekulé-patterned graphene. Mathematically, the WKB approximation is regarded as a solution technique for a differential equation with the highest-order derivative term multiplied by a small parameter. In quantum physics, this small parameter is just Planck’s constant ℏ, so that the solution for a WKB wave function can be sought as a power series of ℏ. However, the WKB for Dirac electrons in graphene is found drastically different from that in the case of a regular Schrödinger electron, and therefore, we expect that the WKB solution for each novel Dirac material will be quite unique compared with all other known solutions. Moreover, the sought wave function for a new material will usually appear as a quickly oscillating function, which is further modulated by a much slower variation from an external potential $V\left(x\right)$.

Among multiple applications of the WKB theory, we would also like to calculate electron transmission, or tunneling probability, and investigate the various bound states of an electronic system [44,45,47]. In contrast to previous studies on tunneling and transport of electrons in Kek-Y graphene [48,49] through complicated numerical computations, WKB methods, on the other hand, provide reliable estimates based on closed-form analytical expressions, which can be easily used to interpret experimental measurements on electron transport.

Studying electron tunneling, often associated with calculating a transmission amplitude, becomes one of the most important issues for each of the newly discovered Dirac materials, mostly because of the so-called Klein paradox, i.e., a full tunneling of chiral Dirac electrons through a square potential barrier, independent of its width and height if the incoming electron moves in the direction perpendicular to the barrier boundary [12,50,51,52]. For a finite-angle collision, however, a complete transmission also appears for an oblique (not a normal) incidence due to presence of Fabry–Perot resonances. Such a Klein-tunneling behavior was investigated for all Dirac-cone materials and nano-ribbons [49,53,54], including transport in Kekulé-distorted graphene nanoribbons [49]. Specifically, a magic transmission was demonstrated for $\alpha$-${\mathcal{T}}_3$ materials and dice lattice [55,56,57,58], in addition to the unusual tunneling behavior of anisotropic Dirac electrons [59,60,61].

The unique properties of electron tunneling and the existence of the Klein paradox will depend on a few factors, and the most crucial one is the energy bandgap between the valence and conduction bands [62,63,64,65,66]. The two-dimensional materials with a gap do not show Klein tunneling. Such a bandgap in graphene and other Dirac lattices could be created by applying an off-resonance circularly-polarized dressing field, which also turns these material into Haldane-like Chern insulators [67,68]. However, linearly polarized irradiation normally does not affect the gap, although it can induce or modify the existing anisotropy [69,70]. Interestingly, the effects of these external fields with different polarizations could further mix with anisotropic and tilted energy dispersions of lattices, such as 1T′MoS₂ [71,72,73].

One of the most crucial applications of WKB equations is investigating the electron tunneling for non-trivial potential barriers, such as trapezoidal barriers [74] or linearly increasing potentials without a sharp boundary. It was demonstrated for all previously known Dirac materials that the Klein paradox persists in such potential profiles. Verification of this property for Kekulé-patterned graphene is one of the goals of the present paper. Since Dirac electrons moving under a linearly increasing potential can acquire a very large kinetic energy, applying WKB approximation to such a case seems to be the most adequate and efficient approach.

The remaining part of our paper is organized as follows. We begin with a brief review of electronic states and their energy dispersions of Kek-Y graphene in Section 2. In Section 3, we formulate a full WKB approximation for a general case of Kek-Y graphene, which includes finding the semi-classical action, longitudinal electron momentum and WKB wave functions. Section 4 focuses on addressing the electron tunneling and Klein paradox in Kekulé-patterned graphene by utilizing our derived Wentzel–Kramers–Brillouin equations. Moreover, we also consider different methods for calculating electron transmissions through different types of potential barriers, and demonstrate how these methods could be applied to Kekulé-patterned graphene. Finally, we discuss and compare the obtained results and present concluding remarks in Section 5.

2. Energy Band Structure and Electronic States in Kek-Y Graphene

First, we briefly review the electronic band structure and eigenfunctions in Kek-Y graphene. For the type of distortion displayed in Figure 1, the Dirac cones will fold on each other, as seen in Figure 2, which results in two states with no bandgap, but with different Fermi velocities.

The low-energy Hamiltonian for Kek-Y graphene can be approximately written as [1,3]:

$${\mathcal{H}}_Y\left(\mathit{k}\right|{\Delta }_0)=\hslash v_F{\widehat{\mathcal{T}}}_0^{\left(2\right)}\otimes \left(\mathit{k}·{\widehat{\mathsf{\Sigma }}}^{\left(2\right)}\right)+\hslash v_F{\Delta }_0\left(\mathit{k}·{\widehat{\mathcal{T}}}^{\left(2\right)}\right)\otimes {\widehat{\mathsf{\Sigma }}}_0^{\left(2\right)},$$

(1)

where $\mathit{k}=\{k_x,k_y\}$ is a wave vector of electrons, $v_F$ is the Fermi velocity in regular non-strained graphene and ${\Delta }_0$ is a dimensionless strain-induced coupling parameter. As ${\Delta }_0\to 0$, we recover the regular graphene Hamiltonian for unfolded K and ${\mathit{K}}^{\prime }$ valleys, i.e., keeping only the first term on the right hand side of Equation (1). We indicate that Equation (1) holds only for a small strain-coupling parameter ${\Delta }_0$, and it may be modified as the stain strength increases. In most cases, however, ${\Delta }_0$ is always assumed small in order to retain graphene signature of the system. Here, we choose ${\Delta }_0$∽0.1 in our calculations below but derive a generalized model by assuming a substantial difference between two Fermi velocities $v_F(1\pm {\Delta }_0)$.

Additionally:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\Sigma }}}_x^{\left(2\right)}& =& {\widehat{\mathcal{T}}}_x^{\left(2\right)}=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right],\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\Sigma }}}_y^{\left(2\right)}& =& {\widehat{\mathcal{T}}}_y^{\left(2\right)}=\left[\begin{array}{cc}0& -i\\ i& 0\end{array}\right],\hfill \end{array}$$

(3)

which act in the spin- and pseudo-spin subspaces, respectively. Correspondingly, ${\widehat{\mathsf{\Sigma }}}_0^{\left(2\right)}$ and ${\widehat{\mathcal{T}}}_0^{\left(2\right)}$ in Equation (1) are two $2\times 2$ unit matrices in the spin- and pseudo-spin subspaces, written as:

$$\begin{array}{ccc}\hfill {\widehat{\mathsf{\Sigma }}}_0^{\left(2\right)}& =& {\widehat{\mathcal{T}}}_0^{\left(2\right)}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right].\hfill \end{array}$$

(4)

Equation (1) is presented by using a Kronecker product of two $2\times 2$ matrices, defined as:

$$\overleftrightarrow{A}\otimes \overleftrightarrow{B}=\left[\begin{array}{cc}{a}_{11}\overleftrightarrow{B}& a_{12}\overleftrightarrow{B}\\ a_{21}\overleftrightarrow{B}& a_{22}\overleftrightarrow{B}\end{array}\right].$$

(5)

More specifically, for the case of Pauli matrices, we have:

$${\widehat{\mathsf{\Sigma }}}_z^{\left(2\right)}\otimes {\widehat{\mathcal{T}}}_z^{\left(2\right)}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& 1\end{array}\right],$$

(6)

$${\widehat{\mathsf{\Sigma }}}_0^{\left(2\right)}\otimes {\widehat{\mathcal{T}}}_z^{\left(2\right)}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& -1\end{array}\right].$$

(7)

Therefore, any block-diagonal (b–d) matrix can be built by following this procedure, namely:

$${\overleftrightarrow{M}}_{\mathrm{b}-\mathrm{d}}={\widehat{\mathsf{\Sigma }}}_0^{\left(2\right)}\otimes \overleftrightarrow{A}=\left[\begin{array}{cccc}{a}_{11}& a_{12}& 0& 0\\ a_{21}& a_{22}& 0& 0\\ 0& 0& a_{11}& a_{12}\\ 0& 0& a_{21}& a_{22}\end{array}\right].$$

(8)

As a result, Equation (1) can take a simple matrix form:

$${\widehat{\mathcal{H}}}_Y\left(\mathit{k}\right|{\Delta }_0)=\hslash v_F\left[\begin{array}{cc}\mathit{k}·{\widehat{\mathsf{\Sigma }}}^{\left(2\right)}& {\Delta }_0(k_x-i{k}_y){\widehat{\mathsf{\Sigma }}}_0^{\left(2\right)}\\ {\Delta }_0(k_x+i{k}_y){\widehat{\mathsf{\Sigma }}}_0^{\left(2\right)}& \mathit{k}·{\widehat{\mathsf{\Sigma }}}^{\left(2\right)}\end{array}\right],$$

(9)

or explicitly,

$${\mathcal{H}}_Y\left(\mathit{k}\right|{\Delta }_0)=\hslash v_F\left[\begin{array}{cccc}0& k_x+i{k}_y& {\Delta }_0(k_x-i{k}_y)& 0\\ k_x-i{k}_y& 0& 0& {\Delta }_0(k_x-i{k}_y)\\ {\Delta }_0(k_x+i{k}_y)& 0& 0& k_x+i{k}_y\\ 0& {\Delta }_0(k_x+i{k}_y)& k_x-i{k}_y& 0\end{array}\right].$$

(10)

The energy dispersions of Kekulé-Y graphene are simply obtained as eigenvalues of the matrix defined in Equation (10), yielding:

$${\epsilon }_Y(s_1,s_2|k,{\Delta }_0)=\hslash v_F(s_1+s_2{\Delta }_0)k,$$

(11)

which represent two Dirac cones with different Fermi velocities $v_F(1\pm {\Delta }_0)$, as displayed in Figure 2, where $s_1=\pm 1$ is a band index while $s_2=\pm 1$ corresponds to two different Fermi velocities due to strain coupling. These two types of subbands are often referred to as “fast” ($s_2=+1$) and “slow” ($s_2=-1$) Dirac cones, and importantly, photo-induced electronic transitions are allowed between these two.

In correspondence with energy dispersions in Equation (11), the normalized wave functions are calculated as:

$${\mathsf{\Psi }}_Y(s_1,s_2|\mathit{k},{\Delta }_0)=\frac{1}{2}\left[\begin{array}{c}{e}^{-i{\theta }_{\mathbf{k}}}\\ s_1\\ s_2\\ s_1{s}_2{e}^{+i{\theta }_{\mathbf{k}}}\end{array}\right],$$

(12)

where ${\theta }_{\mathbf{k}}={tan}^{-1}(k_y/k_x)$ is the angle associated with the electron wave vector k with respect to the x axis. The wave function in Equation (12) consists of two subspinors corresponding to standard graphene eigenstates, namely:

$${\mathsf{\Psi }}_Y(s_1,s_2|\mathit{k},{\Delta }_0)=\frac{1}{2}\left[\begin{array}{c}\hfill {\mathsf{\Psi }}_g\left(s_1\right|\mathit{k})\\ \hfill s_2{\mathtt{e}}^{+i{\theta }_{\mathbf{k}}}{\mathsf{\Psi }}_g\left(s_1\right|\mathit{k})\end{array}\right]=\frac{1}{2}\left[\begin{array}{c}\hfill \left[\begin{array}{c}1\\ s_1{\mathtt{e}}^{+i{\theta }_{\mathbf{k}}}\end{array}\right]\\ \hfill s_2{\mathtt{e}}^{+i{\theta }_{\mathbf{k}}}\left[\begin{array}{c}1\\ s_1{\mathtt{e}}^{+i{\theta }_{\mathbf{k}}}\end{array}\right]\end{array}\right].$$

(13)

Quantum-mechanically, the velocity operators ${\widehat{V}}_x$ for an electron can be determined by the Heisenberg equation, that is:

$${\widehat{V}}_x\equiv \frac{d\widehat{x}}{dt}=\frac{1}{i\hslash }\left[\widehat{x},{\widehat{\mathcal{H}}}_Y\left(\mathit{k}\right|{\Delta }_0)\right],$$

(14)

which relates to an intraband probability current density ${\widehat{J}}_x(s_1,s_2|\mathit{k},{\Delta }_0)$, where ${\mathcal{H}}_Y\left(\mathit{k}\right|{\Delta }_0)$ is presented in Equation (9). Semi-classically, on the other hand, we define another particle velocity $v_x$, determined by:

$$v_x=\frac{1}{\hslash }\frac{\partial {\langle {\mathcal{H}}_Y\left(\mathit{k}\right|{\Delta }_0)\rangle }_{\mathrm{av}}}{\partial k_x},$$

(15)

which relates to an intraband transport current density $j_x=n_0{v}_x$, where $n_0$ is the areal density of electrons in our system, while ${\langle \cdots \rangle }_{\mathrm{av}}$ denotes a quantum-mechanical average with respect to the wave function ${\mathsf{\Psi }}_Y(s_1,s_2|\mathit{k},{\Delta }_0)$ presented in Equation (13).

Generally, from the Schrödinger equation, we know that the intraband probability-current density can be calculated by:

$$J_x(s_1,s_2|\mathit{k},{\Delta }_0)\approx \left(\frac{v_F}{k_F}\right)\mathrm{Im}\left[{\mathsf{\Psi }}^{*}\frac{\partial }{\partial x}\mathsf{\Psi }\right]=v_x{\left|{\mathsf{\Psi }}_Y(s_1,s_2|\mathit{k},{\Delta }_0)\right|}^2,$$

(16)

where we treat a Kekulé particle as a free one. The definition in Equation (15) can be employed for calculating the semi-classical transport current density, yielding:

$$\begin{array}{}\mathrm{(17)}& \hfill j_x(s_1,s_2|\mathit{k},{\Delta }_0)& =& n_0{v}_x=n_0{v}_F\left[{\widehat{\mathcal{T}}}_0^{\left(2\right)}\otimes {\widehat{\mathsf{\Sigma }}}_x^{\left(2\right)}+{\Delta }_0)\left(\mathit{k}·{\widehat{\mathcal{T}}}_x^{\left(2\right)}\right)\otimes {\widehat{\mathsf{\Sigma }}}_0^{\left(2\right)}\right]\hfill \mathrm{(18)}& & =& n_0{v}_F\left[\begin{array}{cccc}0& 1& {\Delta }_0& 0\\ 1& 0& 0& {\Delta }_0\\ {\Delta }_0& 0& 0& 1\\ 0& {\Delta }_0& 1& 0\end{array}\right].\hfill \end{array}$$

An expression similar to Equation (17) should also be obtained for transport current ${\widehat{j}}_y(s_1,s_2|\mathit{k},{\Delta }_0)$.

For a Kekule particle, besides the intraband probability current density $\widehat{J}(s_1,s_2|\mathit{k},{\Delta }_0)$, there exists another interband probability current density due to ${\Delta }_0$ (strain) coupling between two Dirac-cones [75]. By including this interband probability current, the total particle current should be continuous and satisfies the following particle-number conversation law:

$$\nabla ·\mathit{J}(s_1,s_2|\mathit{k},{\Delta }_0)+\frac{\partial }{\partial t}|{\mathsf{\Psi }}_Y(s_1,s_2|\mathit{k},{\Delta }_0){|}^2=0.$$

(19)

Since there are no external sources for the current and each observable quantity is time-independent, every term in Equation (19) equals to zero, leading to:

$$\frac{d}{dx}{J}_x(s_1,s_2|\mathit{k},{\Delta }_0)=0⟶J_x(s_1,s_2|\mathit{k},{\Delta }_0)=\mathrm{const}.$$

(20)

This implies that even under a non-uniform potential $V\left(x\right)$, such as a potential barrier, the total current across the barrier must be conserved. This unique feature plays a crucial role in electron tunneling and transmission calculation. For our considered system, we always assume a translational symmetry along the y-direction, and then $J_y(s_1,s_2|\mathit{k},{\Delta }_0)=$ a constant.

3. Wentzel–Kramers–Brillouin Approximation

In this section, we aim to derive the semi-classical Wentzel–Kramers–Brillouin (WKB) approximation applied to Kek-Y graphene. We begin with the Hamiltonian formula in Equation (1) and generalize it for the case with a spatially non-uniform potential $V\left(x\right)$, including the special case with a square barrier $V\left(x\right)=\mathsf{\Theta }\left(x\right)\mathsf{\Theta }(W_B-x)$, where $W_B$ is the width of a square barrier.

By taking into account the fact that our system is no longer translationally invariant in the x-direction, we take $k_x\to -i\partial /\partial x$ and obtain from Equation (1) that:

$${\mathcal{H}}_Y(x,k_y|{\Delta }_0)=i{v}_F\left[\begin{array}{cccc}V\left(x\right)& -\hslash \partial /\partial x-p_y& {\Delta }_0(-\hslash \partial /\partial x-p_y)& 0\\ -\hslash \partial /\partial x-k_y& V\left(x\right)& 0& {\Delta }_0(-\hslash \partial /\partial x-p_y)\\ {\Delta }_0(-\hslash \partial /\partial x+p_y)& 0& V\left(x\right)& -\hslash \partial /\partial x+p_y\\ 0& {\Delta }_0(-\hslash \partial /\partial x+p_y)& -\hslash \partial /\partial x-p_y& V\left(x\right)\end{array}\right],$$

(21)

which is independent of y, where $p_y=\hslash k_y=\mathrm{const}$ remains conserved.

The key idea of the WKB approximation originates from a series expansion in powers of a small parameter ℏ. Since we will consider a semi-classical particle with a large kinetic energy $\epsilon \left(k\right)$, it is necessary to estimate the order of magnitude for each of the considered quantities. By following Refs. [43,45], it is reasonable to rewrite the Hamiltonian formula in Equation (21), as well as the corresponding eigenvalue equation, in a form with only dimensionless quantities, namely:

$$x⟶\frac{x}{W_B}\equiv \xi \mathrm{and}\epsilon \left(k\right)⟶\frac{\epsilon \left(k\right)}{V_0},$$

(22)

where $V_0$ represents the strength of a potential, e.g., the height of a square-barrier potential. Correspondingly, we also introduce other dimensionless quantities, such as:

$$V\left(x\right)⟶\frac{V\left(x\right)}{V_0}\mathrm{and}{p}_y⟶\frac{v_F{p}_y}{V_0}.$$

(23)

In this way, the spatial derivative now becomes:

$$\partial /\partial x\to \partial /(W_B\partial \xi ).$$

(24)

Such a scheme allows us to modify the Planck constant ℏ, which is our chosen small parameter for a series expansion, leading to:

$$\hslash \to \left(\frac{v_F}{W_B{V}_0}\right)\hslash ,$$

(25)

Consequently, the new energy scale of an incoming particle is large and will not be limited by the Fermi energy of electrons in graphene.

In this study, we follow an established WKB solution for a Dirac equation and expand our sought-after wave function as a series, given by:

$$\begin{array}{cc}\hfill {\mathsf{\Psi }}_Y(s_1,s_2|x,k_y)=& \mathtt{exp}\left[\frac{i}{{\hslash }_0}\mathbb{S}\left(x\right)\right]\sum _{\lambda =0}^{\infty }{(-i{\hslash }_0)}^{\lambda }{\mathsf{\Psi }}_{\lambda }(x,k_y)\hfill \\ \hfill =& \mathtt{exp}\left[\frac{i}{{\hslash }_0}\mathbb{S}\left(x\right)\right]\left[{\mathsf{\Psi }}_0(x,k_y)-i{\hslash }_0{\mathsf{\Psi }}_1(x,k_y)-{\hslash }_0^2{\mathsf{\Psi }}_2(x,k_y)+\cdots \right],\hfill \end{array}$$

(26)

where $\mathbb{S}\left(x\right)$ stands for a semi-classical action in the WKB approximation. Here, our goal is finding a differential equation with respect to x, which connects consecutive terms in the expansion introduced in Equation (26). Explicitly, we would seek the form:

$${\widehat{\mathbb{D}}}_x^{\left(4\right)}\left({\Delta }_0\right)\left[\frac{\partial }{\partial x}{\mathsf{\Psi }}_{\lambda }(s_1,s_2|x,k_y)\right]-{\widehat{\mathbb{O}}}_T(s_1,s_2|x,p_y,{\Delta }_0){\mathsf{\Psi }}_{\lambda +1}(s_1,s_2|x,k_y)=0,$$

(27)

where

$$\begin{array}{cc}\hfill {\widehat{\mathbb{D}}}_x^{\left(4\right)}\left({\Delta }_0\right)=& {\widehat{\mathcal{T}}}_0^{\left(2\right)}\otimes {\widehat{\mathsf{\Sigma }}}_x^{\left(2\right)}+{\Delta }_0{\widehat{\mathcal{T}}}_x^{\left(2\right)}\otimes {\widehat{\mathsf{\Sigma }}}_0^{\left(2\right)}\hfill \\ \hfill =& \left[\begin{array}{cccc}0& 1& {\Delta }_0& 0\\ 1& 0& 0& {\Delta }_0\\ {\Delta }_0& 0& 0& 1\\ 0& {\Delta }_0& 1& 0\end{array}\right]\hfill \end{array}$$

(28)

is proportional to the probability-current operator in Equation (17), $\lambda =0,1,2,3,\cdots$, and ${\mathsf{\Psi }}_{\lambda =-1}(s_1,s_2|x,k_y)\equiv 0$. Furthermore, the transport operator ${\widehat{\mathbb{O}}}_T(s_1,s_2|x,p_y,{\Delta }_0)$ introduced in Equation (27), which connects each two consecutive terms of expansion in Equation (26), is as follows:

$$\begin{array}{cc}& {\widehat{\mathbb{O}}}_T(s_1,s_2|x,p_y,{\Delta }_0)\hfill \\ \hfill =& \xi \left(x\right){\widehat{\mathcal{T}}}_0^{\left(2\right)}\otimes {\widehat{\mathsf{\Sigma }}}_0^{\left(2\right)}+{\widehat{\mathcal{T}}}_0^{\left(2\right)}\otimes \left[{\widehat{\mathsf{\Sigma }}}^{\left(2\right)}·\left\{\frac{\partial }{\partial x}\mathbb{S}\left(x\right),p_y\right\}\right]+{\Delta }_0\left[{\widehat{\mathcal{T}}}^{\left(2\right)}·\left\{\frac{\partial }{\partial x}\mathbb{S}\left(x\right),p_y\right\}\right]\otimes {\widehat{\mathsf{\Sigma }}}_0^{\left(2\right)}\hfill \\ \hfill =& \left[\begin{array}{cccc}\xi \left(x\right)& \partial \mathbb{S}\left(x\right)/\partial x-i{p}_y& {\Delta }_0[\partial \mathbb{S}\left(x\right)/\partial x-i{p}_y]& 0\\ \partial \mathbb{S}\left(x\right)/\partial x+i{p}_y& \xi \left(x\right)& 0& {\Delta }_0[\partial \mathbb{S}\left(x\right)/\partial x-i{p}_y]\\ {\Delta }_0[\partial \mathbb{S}\left(x\right)/\partial x+i{p}_y]& 0& \xi \left(x\right)& \partial \mathbb{S}\left(x\right)/\partial x+i{p}_y\\ 0& {\Delta }_0[\partial \mathbb{S}\left(x\right)/\partial x+i{p}_y]& \partial \mathbb{S}\left(x\right)/\partial x-i{p}_y& \xi \left(x\right)\end{array}\right],\hfill \end{array}$$

(29)

where $\xi \left(x\right)=V\left(x\right)-\epsilon \left(p_y\right)$.

Specifically, for $\lambda =-1$, we obtain:

$${\widehat{\mathbb{O}}}_T(s_1,s_2|x,p_y,{\Delta }_0){\mathsf{\Psi }}_0(s_1,s_2|x,k_y)=0.$$

(30)

For the linear homogeneous equation in Equation (30), a non-trivial solution becomes possible only if its determinant is zero, namely:

$$\begin{array}{cc}& {\left[{\left(\frac{\partial \mathbb{S}\left(x\right)}{\partial x}\right)}^2+p_y^2\right]}^2-2{\Delta }_0^2\left[{\left(\frac{\partial \mathbb{S}\left(x\right)}{\partial x}\right)}^4-p_y^4\right]+{\Delta }_0^4{\left[{\left(\frac{\partial \mathbb{S}\left(x\right)}{\partial x}\right)}^2+p_y^2\right]}^2\hfill \\ & -2{\xi }^2\left(x\right)\left[{\left(\frac{\partial \mathbb{S}\left(x\right)}{\partial x}\right)}^2+p_y^2\right]\left(1+{\Delta }_0^2\right)+{\xi }^4\left(x\right)=0.\hfill \end{array}$$

(31)

As a quartic equation, Equation (31) could be solved analytically, but in general, its solutions are too tedious to write down. Instead, we use the fact that ${\Delta }_0\ll 1$. Therefore, as ${\Delta }_0⟶0$, Equation (31) turns into:

$${\left[{\left(\frac{\partial \mathbb{S}\left(x\right)}{\partial x}\right)}^2+p_y^2\right]}^2-2\left[{\left(\frac{\partial \mathbb{S}\left(x\right)}{\partial x}\right)}^2+p_y^2\right]{\xi }^2\left(x\right)+{\xi }^4\left(x\right)=0$$

(32)

with a set of doubly degenerate solutions for incident energy ${\epsilon }_0$ of an electron, given by:

$$p_x\left(x\right)\equiv \frac{\partial \mathbb{S}\left(x\right)}{\partial x}=\pm \sqrt{{\left[{\epsilon }_0-V\left(x\right)\right]}^2-p_y^2},$$

(33)

which results in:

$$\mathbb{S}\left(x\right)-\mathbb{S}\left(x_0\right)=\underset{x_0}{\overset{x}{\int }}d\zeta \sqrt{{\left[{\epsilon }_0-V\left(\zeta \right)\right]}^2-p_y^2}.$$

(34)

For a finite but small ${\Delta }_0$, we are able to solve Equation (31) by using a perturbation method. Assuming:

$$p_x\left(x\right)=p_x^{\left(0\right)}\left(x\right)+{\Delta }_0^2{p}_x^{\left(1\right)}\left(x\right)+\cdots ,$$

(35)

where $p_x^{\left(0\right)}\left(x\right)=\pm \sqrt{{\left[{\epsilon }_0-V\left(x\right)\right]}^2-p_y^2}$ derived from the zeroth-order equation in Equation (32), we obtain the first-order correction $p_x^{\left(1\right)}\left(x\right)$ (on the order of $\backsimeq {\Delta }_0^2$) from the following equation:

$$\begin{array}{cc}& 4p_x^{\left(0\right)}\left(x\right)p_x^{\left(1\right)}\left(x\right)\left\{{\left[p_x^{\left(0\right)}\left(x\right)\right]}^2+p_y^2\right\}-2\left\{{\left[p_x^{\left(0\right)}\left(x\right)\right]}^4-p_y^2\right\}\hfill \\ & +2\left\{{\left[p_x^{\left(0\right)}\left(x\right)\right]}^2+2p_x^{\left(0\right)}\left(x\right)p_x^{\left(1\right)}\left(x\right)+p_y^2\right\}{\left[{\epsilon }_0-V\left(x\right)\right]}^2=0,\hfill \end{array}$$

(36)

given by:

$$p_x^{\left(1\right)}\left(x\right)=\frac{{\left[p_x^{\left(0\right)}\left(x\right)\right]}^4-p_y^4+{\left[p_x^{\left(0\right)}\left(x\right)\right]}^2{\xi }^2\left(x\right)+p_y^2{\xi }^2\left(x\right)}{2p_x^{\left(0\right)}\left(x\right)\left\{{\left[p_x^{\left(0\right)}\left(x\right)\right]}^2+p_y^2-{\xi }^2\left(x\right)\right\}}.$$

(37)

Equations (35)–(37) define the solution for the x-dependent momentum $p_x\left(x\right)$, while the semi-classical action for a given momentum $p_x\left(x\right)$ can be determined by Equation (34). This is one of the most crucial results in semi-classical theory for Kekulé-patterned graphene.

Now, the WKB wave function could be obtained from Equation (27). If we look for an eigenstate formally in the following form:

$${\mathsf{\Psi }}_{\mathrm{WKB}}(s_1,s_2|x,p_y)=\left[\begin{array}{c}{\mathsf{\Psi }}_A\left(x\right)\\ {\mathsf{\Psi }}_B\left(x\right)\\ {\mathsf{\Psi }}_A^{\prime }\left(x\right)\\ {\mathsf{\Psi }}_B^{\prime }\left(x\right)\end{array}\right],$$

(38)

the relations between its components become:

$$\begin{array}{cc}& \xi \left(x\right){\mathsf{\Psi }}_A\left(x\right)+{\Delta }_0\left[p_x\left(x\right)-i{p}_y\right]{\mathsf{\Psi }}_A^{\prime }\left(x\right)+\left[p_x\left(x\right)-i{p}_y\right]{\mathsf{\Psi }}_B\left(x\right)=0,\hfill \\ & \left[p_x\left(x\right)-i{p}_y\right]{\mathsf{\Psi }}_A\left(x\right)+\xi \left(x\right){\mathsf{\Psi }}_B\left(x\right)+{\Delta }_0\left[p_x\left(x\right)-i{p}_y\right]{\mathsf{\Psi }}_B^{\prime }\left(x\right)=0,\hfill \\ & {\Delta }_0\left[p_x\left(x\right)+i{p}_y\right]{\mathsf{\Psi }}_A\left(x\right)+\xi \left(x\right){\mathsf{\Psi }}_A^{\prime }\left(x\right)+{\Delta }_0\left[p_x\left(x\right)-i{p}_y\right]{\mathsf{\Psi }}_B^{\prime }\left(x\right)=0,\hfill \end{array}$$

(39)

where $p_x\left(x\right)$ is determined from Equations (35)–(37). It is important to point out that the equations in Equation (39) includes strain-induced coupling ${\Delta }_0$ between two sets of wave functions ${\mathsf{\Psi }}_{A,B}\left(x\right)$ and ${\mathsf{\Psi }}_{A,B}^{\prime }\left(x\right)$, giving rise to a similar interband probability current [75] within the semi-classical frame.

In general, the solution for a WKB wave function is quite complicated, but we can model a reliable solution by using the fact that the wave functions in Kek-Y graphene is a combination of two graphene wave function and does not depend on ${\Delta }_0$. For this purpose, we first introduce a notation:

$$\mathsf{\Theta }(x,p_y)=\frac{1}{\xi \left(x\right)}\left[{\pi }_x\left(\xi \right)-i{\pi }_y\right]=-\mathtt{exp}\left[-i\tau {\theta }_{\mathbf{k}}\left(x\right)\right].$$

(40)

where ${\theta }_{\mathbf{k}}\left(\xi \right)={tan}^{-1}[k_x\left(\xi \right)/k_y]$ is the angle between wave vector $k=\{k_x\left(\xi \right),k_y\}=(1/\hslash )\{{\pi }_x\left(\xi \right),{\pi }_y\}$ and the x-axis. Making use of this notation, the wave function now takes a simple form:

$${\mathsf{\Psi }}_Y(s_1,s_2|x,p_y)=\left[\begin{array}{c}\mathsf{\Theta }(x,p_y)\\ s_1\\ s_2\\ s_1{s}_2\mathsf{\Theta }*(x,p_y)\end{array}\right]{\phi }_H^{\left(0\right)}\left(x\right).$$

(41)

Finally, in analogy with graphene study, we obtain the following expression:z

$${\phi }_H^{\left(0\right)}\left(x\right)={\left[\frac{p_x^2\left(x\right)+p_y^2}{p_x^2\left(x\right)}\right]}^{1/2}={\left\{1+{\left[\frac{p_y}{p_x\left(x\right)}\right]}^2\right\}}^{1/2},$$

(42)

and now the wave function in Equation (12) is fully determined.

Up to now, we have eventually obtained the semi-classical action in Equation (34), the x-dependent momentum $p_x\left(x\right)$ and the WKB wave function in Equation (12). This concludes constructing a semi-classical theory for Kek-Y graphene.

4. Electron Tunneling in Kekulé-Patterned Graphene

The primary application of the obtained WKB approximation is a good numerical estimate of the electron transmission amplitude for a fast-moving electron. This method is easily applicable to the cases of non-square potential barriers and a wide class of potential profiles as long as this profile $V\left(x\right)$ is smooth and there is no abrupt change in the particle’s momentum $p\left(x\right)$. At the same time, WKB calculations were used to predict or verify the existence of the Klein paradox for a normal incident particle on a square barrier in graphene and other Dirac cone materials, i.e., they could be formally applied to a square barrier as well.

The transmission amplitude is defined as the ratio of the electron current passing through the barrier to the incoming current. Since the total current in the direction perpendicular to a potential barrier is conserved, the reflect current could also be easily calculated as the difference between them. The biggest loss of a transmitted current occurs in the classically forbidden regions (CFR), in which the longitudinal momentum $p_x\left(x\right)$ becomes purely imaginary and the moving particle is represented by an evanescent wave, whose amplitude is exponentially decreasing along the x-axis. The transmission amplitude $\mathbb{T}\left(p_y\right|a,{\Delta }_0)$ in the WKB theory is estimated by the integral of the purely imaginary momentum $p_x\left(x\right)$ over each of the classically forbidden regions along the x-direction. This results in the following equation:

$$\mathbb{T}\left(p_y\right|a,{\Delta }_0)=\mathtt{exp}\left[-\frac{2}{\hslash }\underset{\mathrm{CFR}}{\int }\left|p_x\left(x\right)\right|dx\right]$$

(43)

as long as the condition ${\displaystyle \underset{\mathrm{CFR}}{\int }\left|p_x\left(x\right)\right|dx\gg \hslash }$ is satisfied. Based on Equation (31), squared longitudinal momentum ${\mathrm{\Pi }}^{\left(S\right)}\left(x\right)=p_x^2\left(x\right)$ is determined from:

$${\left[{\mathrm{\Pi }}^{\left(S\right)}\left(x\right)+p_y^2\right]}^2-2{\Delta }_0^2\left\{{\left[{\mathrm{\Pi }}^{\left(S\right)}\left(x\right)\right]}^2-p_y^4\right\}-2{\left[{\epsilon }_0-V\left(x\right)\right]}^2\left[{\mathrm{\Pi }}^{\left(S\right)}\left(x\right)+p_y^2\right]\left(1+{\Delta }_0^2\right)+{\left[{\epsilon }_0-V\left(x\right)\right]}^4=0.$$

(44)

The calculated locations of CFRs are presented in Figure 3 for ${\Delta }_0=0$ and in Figure 4 for ${\Delta }_0\ne 0$. We see that the former case is pretty much equivalent to graphene and the boundaries of the CFR’s are linear and symmetric with respect to $x=0$. The width of the classically inaccessible region along the x-axis obviously depends on the slope of the potential profile, i.e., a higher slope a corresponds to reduced width of a CFR.

When ${\Delta }_0\ne 0$, the equation with classically forbidden regions is much more complicated and the region itself is extended in a non-trivial way, as seen from Figure 4. However, the boundaries of CFR are still linear and symmetric with respect to $x=0$. This is in stark contrast with gapped $\alpha$-${\mathcal{T}}_3$ materials in which the classically inaccessible regions are very complicated, asymmetric and even involve curved boundaries.

Our results for the transmission are presented in Figure 5. First, we find that the expansion of the classical and accessible region leads to lower transmission so that it is more suppressed for a smaller slope a in the potential profile $V\left(x\right)$, as found in Figure 5a. Similarly, the suppression of transmission $\mathbb{T}\left(p_y\right|a,{\Delta }_0)$ increases for a finite ${\Delta }_0$ due to the appearance of extended CFRs, However, its dependence on ${\Delta }_0$ is much smaller compared to that on the potential profile slope a, as can be verified by Figure 5a.

It is important to mention that the WKB approximation provides a fairly satisfactory estimate for the electron transmission in the cases as standard techniques are either not applicable or would lead to tremendously complicated equations. This is exactly the case for Kek-Y graphene in which similarly to tunneling of Rashba electrons, there are two non-equivalent solutions for wave vector k at a given energy [76]. However, in the case of Kek-Y graphene, we would need to solve a quartic equation so as to obtain the k value for a given energy ${\epsilon }_0$, and then use these expressions to find the transmission, which seems hardly possible.

As mentioned above, the main application of the WKB theory in any existing material is investigating electron tunneling over various potential barriers or regions in which the potential $V\left(x\right)$ of the incoming particle is larger than its kinetic energy so that a classical particle could not even be present in such regions. These locations are also referred to as classically forbidden or classically inaccessible regions. Another possible application of WKB is investigating electronic bound states with quantized momentum, which are rarely considered for the Dirac lattices.

A general condition for applicability of WKB theory for graphene or another Dirac cone material is given in terms of the momentum $p\left(x\right)$ as:

$$\hslash \left|\frac{dp\left(x\right)}{dx}\right|\ll p^2\left(x\right),$$

(45)

which is evidently equivalent to a slow change of the amplitude of a de Broglie wavelength $\left]lambda\right(x)$ of a particle:

$$\frac{p{\left(x\right)}^2}{p_x{\left(x\right)}^2}\left|\frac{d\lambda \left(x\right)}{dx}\right|\ll 1.$$

(46)

The change of the electron momenta is directly related to the slope of the potential energy profile, such as:

$$\left|\frac{dp\left(x\right)}{dx}\right|=\frac{m}{p\left(x\right)}\left|\frac{dV\left(x\right)}{dx}\right|$$

(47)

for a regular Schrodinger electron or

$$\left|\frac{dp\left(x\right)}{dx}\right|=\frac{1}{{\upsilon }_F}\left|\frac{dV\left(x\right)}{dx}\right|$$

(48)

for graphene, with Fermi velocity ${\upsilon }_F$, which implies two important limitations of any WKB approximation: a sharp, sudden change of the potential $V\left(x\right)$ and the regions in the vicinity of the turning points, where $p\left(x\right)\to 0$. At the same time, all the results for Klein tunneling over a square barrier with sharp edges and $V\left(x\right)=V_0\mathsf{\Theta }\left(x\right)\mathsf{\Theta }(W_B-x)$, where $W_B$ is constant, are faithfully reproduced by WKB theory, which appears to be valid much beyond its known range of applicability. At the time, WKB methods provide a straightforward and reliable estimate of the electron transmission amplitudes for the cases when traditional techniques are impossible or extremely complicated, such as for our case of Kek-Y-patterned graphene.

The validity of the WKB method is determined not by the size or height of the potential barrier. However, we should require the potential profile to be smooth and apply the approximation far from the turning points where the momentum of our particle becomes equal or close to zero.

5. Concluding Remarks

In conclusion, we have built a complete semi-classical WKB theory for recently discovered Kek-Y-patterned graphene with two coupled gapless Dirac cones with different Fermi velocities $v_F(1\pm {\Delta }_0)$. This includes finding the semi-classical action, x-dependent longitudinal momentum $p_x\left(x\right)$ and leading-order WKB wave function. Meanwhile, we have also derived and solved the transport equation, which connects every two successive orders of a sought wave function.

Our obtained WKB equations are further applied to calculate the electron transmission amplitudes based on classically inaccessible regions. We have demonstrated that the slope a of a linear potential $V\left(x\right)=ax$ can strongly affect electron transmission, while the coupling parameter ${\Delta }_0$ has a very limited effect on transmission, because it only slightly modifies the location and size of a classically forbidden region for a specified incoming electron.

Our obtained results have numerous applications to constructing and operating a new generation of electronic devices, since the electron transmission amplitude is directly related to the quantum transport in our material. Most of the known models for the electron conductivity, such as Kubo formalism, require the exact knowledge of electron transmission. Applied irradiation could induce or modify the existing bandgap and greatly affect the electron transmission, which is crucial for optoelectronic applications. We are strongly convinced that this work reveals a crucial step in advancing our knowledge on the electronic properties of newly discovered Dirac materials. Our obtained WKB theory provides a powerful tool for systematically investigating crucial electronic and transport properties of novel Dirac materials, which have numerous applications in modern electronic and opto-electronic devices, and especially in valleytronics.