# Developing a semiclassical Wentzel–Kramers–Brillouin theory for
$\alpha-\mathcal{T}_{3}$ model

Kathy Blaise1, Chinedu Ejiogu1, Andrii Iurov1111E-mail contact:
<EMAIL_ADDRESS><EMAIL_ADDRESS>, Liubov Zhemchuzhna1,2,
Godfrey Gumbs2,3, and Danhong Huang4,5 1Department of Physics and Computer
Science, Medgar Evers College of City University of New York, Brooklyn, NY
11225, USA
2Department of Physics and Astronomy, Hunter College of the City University of
New York, 695 Park Avenue, New York, New York 10065, USA
3Donostia International Physics Center (DIPC), P de Manuel Lardizabal, 4,
20018 San Sebastian, Basque Country, Spain
4Space Vehicles Directorate, US Air Force Research Laboratory, Kirtland Air
Force Base, New Mexico 87117, USA
5Center for High Technology Materials, University of New Mexico, 1313 Goddard
SE, Albuquerque, New Mexico, 87106, USA

###### Abstract

We have developed a complete semiclassical Wentzel-Kramers-Brillouin (WKB)
theory for $\alpha-\mathcal{T}_{3}$ model which describes a wide class of
existing pseudospin-1 Dirac cone materials. By expanding the sought wave
functions in a series over the powers of Planck constant $\hbar$, we have
obtained the leading order expansion term which is the key quantity required
for calculating the electronic and transport properties of a semiclassical
electron in $\alpha-\mathcal{T}_{3}$. We have derived the transport equations
connecting each two consecutive orders of the wave function expansion and
solved them to obtained the first order WKB wavefunction. We have also
discussed the applicability of the obtained approximation and how these
results could be used to investigate various tunneling and transport
properties of $\alpha-\mathcal{T}_{3}$ materials with non-trivial potential
profiles. Our results could be also helpful for constructing electronics
devices and transistors based on innovative flat-band Dirac materials.

##  I Introduction

The quantum-mechanical description of realistic electronic behavior in various
existing materials is predictably a lot more complicated than any classical
analysis of such phenomena. Therefore, it is important to build a
semiclassical approximation for high-energy and fast-moving electrons with a
simplified description which would still lead to a physically correct picture.
Therefore, creating a semiclassical WKB theory is one of the most important
ingredients for any new type of Hamiltonian and we believe that doing so for a
newly discovered $\alpha-\mathcal{T}_{3}$ and a dice lattice would be deemed
as a crucial advancement in addition to the current research being performed
in the field low-dimensional condensed matter physics. Vogl et al. (2017);
Vandecasteele et al. (2010); Zhang et al. (2012); Weekes et al. (2021);
Zalipaev et al. (2015); Zalipaev (2011); Berk et al. (1982)

Generally, the idea of a WKB approximation is based on finding an approximated
solution of a differential equation with spatially dependent coefficients in
the case when one of those them (which physically corresponds to the external
potential $V(x)$) is varying considerable slower than the others which related
to a the rate of change of the wavelength for our electron. Therefore, the
obtained solutions is interpreted as a rapidly oscillating a quantum state
modulated by a non-essential change of the external potential. The WKB
solution is obtained by expanding the wavefunction over the powers of Planck
constant $\hbar$ and in this respect is viewed as semiclassical. It was
previously obtained for all Dirac materials, specifically graphene and gapped
graphene. WKB theory has numerous applications, specifically, in calculating
the electron tunneling for nontrivial potential barriers, Sonin (2009); Anwar
et al. (2020) investigating into the resonant tunneling and resonant
scattering, trapped and localized electronic states. Ye et al. (2020); Xu and
Lai (2019); Gangaraj et al. (2020)

Specifically, the present work is intended to develop a complete semiclassical
WKB theory for the so-called $\alpha-\mathcal{T}_{3}$ model. Bercioux et al.
(2009); Dóra et al. (2011); Vidal et al. (1998, 2001); Illes (2017) This is a
general model for a large group of recently discovered and innovative low-
dimensional materials known by the presence of a flat band in their energy
dispersions in addition to a regular Dirac cone. This flat band is stable and
remains in the presence of a number of external perturbations, including
magnetic and electric fields, point defects and off-resonance dressing fields.
Iurov et al. (2019); Dey and Ghosh (2018, 2019); Goldman and Dalibard (2014);
Iurov et al. (2013); Kristinsson et al. (2016); Kibis (2010); Iurov et al.
(2017); Sandoval-Santana et al. (2020); Kibis et al. (2017); Kibis (2022)

The atomic structure of $\alpha-\mathcal{T}_{3}$ materials is represented by a
regular hexagon lattice similar to graphene with an additional atom in the
center of each hexagon. This buildup leads to the existence of extra electron
hopping coefficient associated with the hub atom which results in a
dispersionless band in the low-energy bandstructure of these materials. The
effect and the ”weight” of the electron transition from and into this flat
band is defined by a quantum phase $\phi$ and ranges from 0 (graphene with a
completely uncoupled flat band) and a dice lattice with $\phi=\pi/4$ and the
strongest possible effect of the flat band. These highly unusual energy
dispersions result in unique and truly fascinating electronic, Cunha et al.
(2022); Malcolm and Nicol (2016); Iurov et al. (2020a); Abranyos et al.
(2020); Islam and Dutta (2017); Huang et al. (2019); Cunha et al. (2021);
Iurov et al. (2021) transport, Urban et al. (2011); Illes and Nicol (2017);
Wang et al. (2020); Louvet et al. (2015) optical, Han and Lai (2022); Carbotte
et al. (2019); Iurov et al. (2020b); Kovács et al. (2017) and magnetic, Illes
et al. (2015); Illes and Nicol (2016); Raoux et al. (2014); Piéchon et al.
(2015); Biswas and Ghosh (2016, 2018); Balassis et al. (2020) properties of
these innovative materials, which are very different from those in graphene.
Neto et al. (2009) Together with tilted Dirac cone materials Tan et al.
(2020); Verma et al. (2017); Champo and Naumis (2019) and topological Dirac
semimetals Islam and Zyuzin (2019); Islam and Saha (2018),
$\alpha-\mathcal{T}_{3}$ materials are considered the most promising and
innovating low-dimensional lattices at the preset time.

At this moment, a large number of $\alpha-\mathcal{T}_{3}$ materials and dice
lattices have been fabricated. These include a trilayer SrTiO3/SrIrO3/SrTiO3,
Qiu et al. (2016); Wang and Ran (2011) Hg1-xCdx quantum well, Orlita et al.
(2014); Malcolm and Nicol (2016) Josephson arrays, Han and Lai (2022);
Ahmadkhani and Hosseini (2020) Leib and Kagome optical lattices and optical
waveguides. Jo et al. (2012); Ruostekoski (2009); Santos et al. (2004); Huang
et al. (2011); Mukherjee et al. (2015); Li et al. (2015); Vicencio et al.
(2015); Baba (2008); Romhányi et al. (2015), In0.53Ga0.47As/InP semiconducting
layers Franchina Vergel et al. (2020) and many-many other existing materials.
A comprehensive review of all Dirac materials with a flat band can be found in
Ref. [Leykam et al., 2018].

To sum up, we would like to refer a reader to a complete and timely review of
all existing flat-band structures published in 2019. Leykam et al. (2018)
Pseudospin-1 Dirac-Weyl Hamiltonian used to model the electronic behavior for
$\alpha-\mathcal{T}_{3}$ model describes the properties of all these materials
with a certain degree of similarity. Iurov et al. (2020c); Li et al. (2017)

Both theoretical and experimental investigation into electronic properties has
become an important direction of nowadays condensed matter physics, nearly all
possible issues related to the electronic phenomena in
$\alpha-\mathcal{T}_{3}$ have been now addressed. Thus, we feel that building
a semiclassical WKB theory for such new type of materials is one of the
crucial problems which still needs to be solved.

The rest of the present paper is organized in the following way. First, in
Sec. II we review the low-energy Hamiltonian, bandstructure and the electronic
states for various types of $\alpha-\mathcal{T}_{3}$ lattices, including
graphene and a dice lattice. Section III, which is the key section of this
paper, deals with calculating the semiclassical action and the longitudinal
electron momentum, finding the leading order wavefunctions, deriving the most
general system of transport equations which connect each to consequent orders
of the $\hbar$-expansion of the wave function, calculating the first order
wave function and verifying the applicability of the WKB approximation for our
model Hamiltonian. In Sec. IV, we briefly discuss the applications of the
obtained semiclassical wave function to studying the electron tunneling and
estimating the transmission coefficient for various non-trivial potential
potential profiles. Finally, we provide some concluding remarks and outlook in
Sec. V.

##  II Electronic states in pseudospin-1 Dirac cone materials

We are going to build and use a semiclassical WKB approximation for
$\alpha-\mathcal{T}_{3}$ model with the low-energy electronic states described
by a pseudospin-1 low-energy Dirac-Weyl Hamiltonian

$\hat{\mathcal{H}}_{\alpha}(\mbox{\boldmath$k$}\,|\,\tau,\phi)=\hbar
v_{F}\left[\begin{array}[]{ccc}0&k^{\tau}_{-}\,\cos\phi&0\\\
k^{\tau}_{+}\,\cos\phi&0&k^{\tau}_{-}\,\sin\phi\\\
0&k^{\tau}_{+}\,\sin\phi&0\end{array}\right]\,,$ (1)

where phase $\phi$ is obtained from the relative hopping parameter $\alpha$ as
$\alpha=\tan\phi$ and $k^{\tau}_{\pm}=\tau k_{x}\pm ik_{y}$ depend on the
valley index $\tau=\pm 1$ corresponding to $K$ and $K^{\prime}$ valleys.

For $\phi=0$, Hamiltonian (1) is obviously reduced to that of graphene. The
opposite limiting case $\phi=\pi/4$ defined as a dice lattice, corresponds to
the strongest effect of the presence of the additional hub atom and is always
selected as a separate case if of the highest importance. We will explicitly
provide all our results for the dice lattice.

If the electron is also subjected to a finite position-dependent potential
$V(x)$, such as a square barrier $V(x)=V_{0}\Theta(x)\Theta(W_{B}-x)$ or a
barrier with non-piecewise-uniform potential profile, the longitudinal
momentum $p_{x}(x)$ also depends on $x$ and Hamiltonian (1) could be rewritten
as

$\hat{\mathcal{H}}_{\alpha}(\mbox{\boldmath$k$}\,|\,\tau,\phi)=v_{F}\,\hat{\mbox{\boldmath$\Sigma$}}(\phi)\cdot\left\\{-i\hbar\mbox{\boldmath$\nabla$}_{\tau}\right\\}+V(x)\,\hat{\Sigma}^{(3)}_{0}\
,$ (2)

in terms of $\phi$-dependent $3\times 3$ Pauli matrices
$\hat{\Sigma}^{(3)}(\phi)=\left\\{\hat{\Sigma}^{(3)}_{x}(\phi),\,\hat{\Sigma}^{(3)}_{y}(\phi)\right\\}$

$\hat{\Sigma}_{x}(\phi)=\left[\begin{array}[]{ccc}0&\cos\phi&0\\\
\cos\phi&0&\sin\phi\\\ 0&\sin\phi&0\end{array}\right]\ ,$ (3)

and

$\hat{\Sigma}_{y}(\phi)=i\,\left[\begin{array}[]{ccc}0&-\cos\phi&0\\\
\cos\phi&0&-\sin\phi\\\ 0&\sin\phi&0\end{array}\right]\ .$ (4)

A $3\times 3$ unit matrix $\hat{\Sigma}^{(3)}_{0}$

$\hat{\Sigma}^{(3)}_{0}=\left[\begin{array}[]{ccc}1&0&0\\\ 0&1&0\\\
0&0&1\end{array}\right]\ ,$ (5)

which also enters Eq. (2), does not depend on phase $\phi$, and
$\mbox{\boldmath$\nabla$}_{\tau}=\\{\tau\partial/\partial
x,\,\partial/\partial y\\}$. For a dice lattice with $\phi=\pi/4$ matrices
Eqs. (3) and (4) are simplified to well-known $3\times 3$ Pauli matrices

$\hat{\Sigma}^{(3)}_{x}=\frac{1}{\sqrt{2}}\,\left[\begin{array}[]{ccc}0&1&0\\\
1&0&1\\\ 0&1&0\end{array}\right]\ .$ (6)

and

$\hat{\Sigma}^{(3)}_{y}=\frac{i}{\sqrt{2}}\,\left[\begin{array}[]{ccc}0&-1&0\\\
1&0&-1\\\ 0&1&0\end{array}\right]\ .$ (7)

For a finite bandgap the third Pauli matrix

$\hat{\Sigma}^{(3)}_{z}=\left[\begin{array}[]{ccc}1&0&0\\\ 0&0&0\\\
0&0&-1\end{array}\right]\ ,$ (8)

could be utilized. Gorbar et al. (2019, 2021)

The energy eigenvalues of Hamiltonian (1) are

$\varepsilon^{\gamma=\pm 1}_{\tau,\,\phi}(\mbox{\boldmath$k$})=\gamma\hbar
v_{F}k$ (9)

where $\gamma=-1$ ($\gamma=+1$) corresponds to the valance (conduction) band,
and $\gamma=0$ leads to a dispersionless solution

$\varepsilon^{\gamma=0}_{\tau,\,\phi}(\mbox{\boldmath$k$})=0\ .$ (10)

known as a ”flat band”. All three eigenenergies of (1) do not depend on phase
$\phi$.

for the remaining flat (or dispersionless) band. Here, all three bands in Eqs.
(9) and (10) do not show any dependence on phase $\phi$ (or parameter
$\alpha$).

The wave functions corresponding to the valence and conduction bands in Eq.
(9) are given as

$\Psi^{\gamma=\pm
1}_{\tau,\,\phi}(\mbox{\boldmath$k$})=\frac{1}{\sqrt{2}}\left[\begin{array}[]{c}\tau\cos\phi\,\,\texttt{e}^{-i\tau\theta_{\bf
k}}\\\ \gamma\\\ \tau\sin\phi\,\,\texttt{e}^{+i\tau\theta_{\bf
k}}\end{array}\right]\ ,$ (11)

where $\theta_{\bf k}=\arctan(k_{y}/k_{x})$ is the angle of wave vector
$\mbox{\boldmath$k$}=\\{k_{x},k_{y}\\}$ made with the $x$-axis. The remaining
wave function attributed to the flat band is

$\Psi^{\gamma=0}_{\tau,\,\phi}(\mbox{\boldmath$k$})=\left[\begin{array}[]{c}\sin\phi\,\,\texttt{e}^{-i\tau\theta_{\bf
k}}\\\ 0\\\ -\cos\phi\,\,\texttt{e}^{+i\tau\theta_{\bf k}}\end{array}\right]\
.$ (12)

Here, we would like to indicate that the energy bands in Eqs. (9) and (10), as
well as the wave functions in Eqs. (11) and (12), are obtained for a
spatially-uniform potential independent of position coordinates $x$ and $y$.

Figure 1: (Color online) Schematics for a Dirac electron moving in a varying
potential $V(x)$. The longitudinal momentum for such electron also depends on
its position $x$ while the transverse momentum $k_{y}$ is conserved. We also
demonstrate classically forbidden or classically inaccessible regions with
$p_{x}(x)^{2}<0$ and its boundaries known as the turning points. The crossing
point is a location of the particle for which $V(x)=E_{\gamma}$ so that the
transition between the electron and hole states is observed.

For a dice lattice, wave functions (11) and (12) are as follows:

$\Psi^{\gamma=\pm
1}_{\tau,\,\phi=\pi/4}(\mbox{\boldmath$k$})=\frac{1}{2}\,\left[\begin{array}[]{c}\texttt{e}^{-i\tau\theta_{\bf
k}}\\\ \sqrt{2}\,\tau\gamma\\\ \texttt{e}^{+i\tau\theta_{\bf
k}}\end{array}\right]\ ,$ (13)

and

$\Psi^{\gamma=0}_{\tau,\,\phi=\pi/4}(\mbox{\boldmath$k$})=\frac{1}{\sqrt{2}}\,\left[\begin{array}[]{c}\texttt{e}^{-i\tau\theta_{\bf
k}}\\\ 0\\\ -\texttt{e}^{+i\tau\theta_{\bf k}}\end{array}\right]\ ,$ (14)

in which two components have equal magnitudes.

##  III Semi-classical WKB model and solution

Here, we intent to derive the so-called transport equations connecting
different orders of the series expansion of the WKB wave function in powers of
$\hbar$, This is a set of the most important differential equation which allow
for calculation the wave function with any required precision. The most
important is the principal zero-order wave function which requires the
knowledge of the semi-classical action being the spatial derivative of the
position-dependent longitudinal electron momentum of our electron.

We also provide the general solutions for our transport equations, calculate
the next (first-order) expansion term of the electron wave functions and
demonstrate the applicability conditions for WKB approximation
$\psi_{1}\ll\psi_{0}$. This is the central section of our paper containing its
most crucial findings.

###  III.1 General solution and action

Our main consideration is an electron moving in $x$-dependent potential $V(x)$
so that the longitudinal momentum $p_{x}(x)$ is also non-uniform along the
$x$-direction. The wave function is presented as
$\Psi(x,y)\backsim\psi(x)\,\texttt{e}^{ik_{y}\,y}$ due to the fact that the
translational symmetry in the $y$-direction is still preserved.

Hamiltonian in Eq. (2) is now given by

$\displaystyle\hat{\mathcal{H}}_{\alpha}(x,k_{y}\,|\,\tau)=\left\\{\begin{array}[]{ccc}V(x)&v_{F}\,\cos\phi\,\Xi(x)_{\tau}^{-}&0\\\
v_{F}\,\cos\phi\,\Xi(x)_{\tau}^{+}&V(x)&v_{F}\,\sin\phi\,\Xi(x)_{\tau}^{-}\\\
0&v_{F}\,\sin\phi\,\Xi(x)_{\tau}^{+}&V(x)\end{array}\right\\}\,,$ (18)
$\displaystyle\Xi(x)_{\tau}^{\pm}=-i\hbar\,\tau\,\frac{\partial}{\partial
x}\pm ip_{y}\,,$

which reflect the remaining translational symmetry in the $y-$direction. Also,
we used a $p_{y}=\hbar k_{y}$ notation. We describe this situation in Fig. 1.
If the potential energy of the particle is increased, its longitudinal
momentum will become zero at some points which are defined as ”turning
points”. These points specify the classically forbidding regions - the areas
with imaginary longitudinal momentum $p_{x}(x)^{2}<0$. The point where
$E_{\gamma}=V(x)$ and the particle has zero energy is the crossing point
between the electron and hole states. In contrast to graphene, the electron
momentum in $\alpha-\mathcal{T}_{3}$ at such points cannot be determined due
to its infinite degeneracy for $E_{\gamma}-V(x)=0$.

Since we are planning to build up an approximation and estimate each of its
terms, it seems reasonable to rewrite our initial eigenvalue equation in a
dimensionless form. A similar approach was adopted in Refs. [Zalipaev et al.,
2015] and [Weekes et al., 2021]. First, we divide each of equations (18) by
the unit of energy $V_{0}$ as $E\longrightarrow E/V_{0}$ and
$V(x)\longrightarrow V(x)/V_{0}$. $V_{0}$ could be a maximum height of the
potential barrier or the potential energy of our particle at an any fixed
point. The length is renormalized over the barrier width $W_{B}$ as
$x\longrightarrow x/W_{B}$. Thus, the new dimensionless momentum is obtained
as $p_{x,y}\longrightarrow v_{y}p_{x,y}/V_{0}$ and, most importantly, our
expansion parameter the Planck constant is modified as

$\hbar\longrightarrow\left(\frac{v_{F}}{W_{B}\,V_{0}}\right)\ \hbar\ .$ (19)

The spatial derivative would change as $\partial/\partial
x\rightarrow\partial/(W_{B}\,\partial x)$. We feel reluctant to connect the
energy unit through the inverse length (a unit of momentum) is it would bring
another term dependent on a small parameter $\hbar$ into our calculation.
Thus, the kinetic energies of the incoming particles are assumed to be quasi-
classical, e.i., having much larger values than the Fermi energy of our
system.

Our main eigenvalue equation is now written as

$\displaystyle\hat{\mathcal{H}}_{\alpha}(x,p_{y}\,|\,\tau)\,\Psi_{\alpha}^{\,\gamma}(x,p_{y}\,|\,\tau)=E_{\gamma}\,\Psi^{\,\gamma}(x,p_{y}\,|\,\phi,\tau)$
(20)

where the wavefunction $\Psi_{\alpha}^{\,\gamma}(x,p_{y}\,|\,\tau)$ is
expressed as

$\displaystyle\Psi_{\alpha}^{\,\gamma}(x,p_{y}\,|\,\tau)=\psi_{\alpha}^{\,\gamma}(x\,|\,\phi,\tau)\,\texttt{exp}\left(\frac{i}{\hbar}\,p_{y}y\right)=\left[\begin{array}[]{c}\phi_{A}(x\,|\,\phi,\tau)\\\
\phi_{H}(x)\\\
\phi_{B}(x\,|\,\phi,\tau)\end{array}\right]\texttt{exp}\left(\frac{i}{\hbar}\,p_{y}y\right)\
.$ (24)

Hamiltonian in Eq. (20) now becomes

$\displaystyle\hat{\mathcal{H}}_{\alpha}(x,p_{y}\,|\,\tau)$ $\displaystyle=$
$\displaystyle\hat{\Sigma}^{(3)}_{0}\,V(x)+\hat{\Sigma}^{(3)}_{x}(\phi)\,\left[-i\hbar\,\tau\,\frac{\partial}{\partial
x}\right]+\hat{\Sigma}^{(3)}_{y}(\phi)p_{y}=$ (28) $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ccc}V(x)&v_{F}\,\Xi(x)_{\tau}^{-}\,\cos\phi&0\\\
v_{F}\,\Xi(x)_{\tau}^{+}\,\,\cos\phi&V(x)&v_{F}\,\Xi(x)_{\tau}^{-}\,\sin\phi\\\
0&v_{F}\,\Xi(x)_{\tau}^{+}\,\sin\phi&V(x)\end{array}\right\\}\,.$

A standard Wentzel–Kramers–Brillouin semiclassical approach is developed by
expanding the unknown wave function in a series over the powers of $\hbar$

$\displaystyle\Psi(\xi,p_{y}\,|\,\phi,\tau)=\texttt{exp}\left\\{\frac{i}{\hbar}\,\mathbb{S}(x)\right\\}\,\sum\limits_{\lambda=0}^{\infty}(-i\hbar)^{\lambda}\,\Psi_{\lambda}(x)=$
(29)
$\displaystyle=\texttt{exp}\left\\{\frac{i}{\hbar}\,\mathbb{S}(x)\right\\}\,\left[\Psi_{0}(x)-i\hbar\,\Psi_{1}(x)-\hbar^{2}\,\Psi_{2}(x)+\cdots\,\right]\
,$

where $\mathbb{S}(x)$ is the semi-classical action, which we still need to
calculate.

The problem of building a semiclassical approximation for an $x-$dependent
potential will be solve once we obtain the WKB wavefunctions or all the terms
of expansion (29). Therefore, we need to derive a set of differential
equations which connect all consecutive terms of expansion (29). These
equations, referred to as the transport equations, are the key relations in
any semiclassical theory.

Explicitly, the transport equation for $\alpha-\mathcal{T}_{3}$ is derived as

$\hat{\Sigma}^{(3)}_{x}(\phi)\,\left\\{\,\frac{\partial}{\partial
x}\,\Psi_{\lambda}(x,p_{y}\,|\,\phi,\tau)\,\right\\}-\frac{1}{\sqrt{2}}\,\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi\,\tau)\,\Psi_{\lambda+1}(x,p_{y}\,|\,\phi,\tau)=0\
,$ (30)

where $\lambda=0,\,1,\,2,\,3,\,\cdots$ and

$\Psi_{\lambda<0}(\xi,p_{y}\,|\,\phi,\tau)\equiv 0\,.$ (31)

Here, the transport operator
$\hat{\mathbb{O}}_{\,T}(\xi,p_{y}\,|\,\phi,\,\tau)$, connecting consequent
terms of expansion in Eq. (29), is easily found to be

$\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\tau)=\left\\{\begin{array}[]{ccc}\nu(x)&\Xi_{\,\mathbb{S},\tau}^{(-)}\,\cos\phi&0\\\
\Xi_{\,\mathbb{S},\tau}^{(+)}\,\cos\phi&\nu(x)&\Xi_{\,\mathbb{S},\tau}^{(-)}\,\sin\phi\\\
0&\Xi_{\,\mathbb{S},\tau}^{(+)}\,\sin\phi&\nu(x)\end{array}\right\\}\,,$ (32)

where

$\displaystyle\Xi_{\,\mathbb{S},\tau}^{(\pm)}=\tau\,\frac{\partial\mathbb{S}(x)}{\partial
x}\pm ip_{y}\,,$ (33) $\displaystyle\nu(x)=V(x)-E_{\gamma}\,.$

Specifically, by setting $\lambda=-1$ in Eq. (30), we obtain a simplified
equation for the zeroth order wavefunction

$\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\tau)\,\Psi_{0}(x,p_{y}\,|\,\phi,\,\tau)=0\
,$ (34)

which is a linear and homogeneous algebraic system; a non-trivial solution for
such such system could be found only if its determinant equals to zero
resulting in the following equation

$-\left[E_{\gamma}-V(x)\right]\,\left\\{\,\left(\frac{\partial\,\mathbb{S}(x)}{\partial
x}\right)^{2}+p_{y}^{2}-\left[E_{\gamma}-V(x)\right]^{2}\,\right\\}=0\ ,$ (35)

which is independent of $\phi$. For a non-uniform potential profile $V(x)$,
condition $E_{\gamma}=V(x)$ could not be always kept, apart from a finite
number of isolated points. Therefore, the classical action $\mathbb{S}(x)$

$\displaystyle\mathbb{S}(x)=\mathbb{S}(x_{0})+\int\limits_{x_{0}}^{x}p_{x}(\eta)\,d\eta\
,$ (36) $\displaystyle
p_{x}(x)=\pm\sqrt{\left[E_{\gamma}-V(x)\right]^{2}-p_{y}^{2}}$ (37)

is expressed through the position-dependent longitudinal momentum $p_{x}(x)$.
Transverse momentum $p_{y}$ always remains the same due to the translational
symmetry.

###  III.2 Leading order wave function

The principal step in building a semiclassical approximation is calculating
the leading order (zero-order) electron wavefunction
$\Psi_{0}(x,p_{y}\,|\,\phi,\tau)$ which has the following form

$\displaystyle\Psi_{0}(x,p_{y}\,|\,\phi,\tau)=\Phi_{0}(x,p_{y}\,|\,\phi,\tau)\,,$
(38)
$\displaystyle\Phi_{0}(x,p_{y}\,|\,\phi,\tau)=\left[\begin{array}[]{c}\varphi_{A}^{(0)}(x\,|\,\phi,\tau)\\\
\varphi_{H}^{(0)}(x)\\\ \varphi_{B}^{(0)}(x\,|\,\phi,\tau)\end{array}\right]\
.$ (42)

Using Eqs. (32) and (36), we can connect the components of wave function (42)

$\displaystyle\nu(x)\,\varphi_{A}^{(0)}(x)+\cos\phi\,\left[\tau
p_{x}(x)-ip_{y}\right]\,\varphi_{H}^{(0)}(x)$ $\displaystyle=$ $\displaystyle
0\ ,$ (43) $\displaystyle\sin\phi\,\left[\tau
p_{x}(x)+ip_{y}\right]\,\varphi_{H}^{(0)}(x)+\nu(x)\,\varphi_{B}^{(0)}(x)$
$\displaystyle=$ $\displaystyle 0\ ,$ (44)

which now accepts the following form

$\displaystyle\Psi_{0}(x\,|\,\phi,\,\tau)=\left[\begin{array}[]{c}\cos\phi\,\,\Theta(x\,|\,\tau)\\\
-1\\\
\sin\phi\,\,\Theta^{\star}(x\,|\,\tau)\end{array}\right]\,\varphi_{H}^{(0)}(x)\
,$ (48) $\displaystyle\Theta(x\,|\,\tau)=\frac{1}{\nu(x)}\,\left[\tau
p_{x}(x)-ip_{y}\right]=-\tau\,\texttt{exp}\left[-i\tau\,\theta_{\bf
k}(x)\right]\,,$ (49)
$\displaystyle\Theta^{\star}(x\,|\,\tau)\xRightarrow{i\rightarrow-i}\Theta(x\,|\,\tau)$

Angle $\theta_{\bf k}(x)=\tan^{-1}[k_{x}(x)/k_{y}]$ associated with the
wavevector ${\bf k}=\\{k_{x}(x),k_{y}\\}$ now depends on coordinate $x$, e.i.,
does not remain the same as the electron moves in potential profile $V(x)$.
Another $x-$dependent scalar function $\varphi_{H}^{(0)}(x)$ from (48) remains
undefined at the moment, we will need to determine it from general transport
equation (30).

Next, we connect the zero and the first order wavefunctions
$\Psi_{1}(x,p_{y}\,|\,\phi,\tau)$ and $\Psi_{1}(x,p_{y}\,|\,\phi,\tau)$ by
taking $\lambda=0$ in Eq. (30)

$\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\,\tau)\,\Psi_{1}(x,p_{y}\,|\,\phi,\tau)=\sqrt{2}\,\hat{\Sigma}^{(3)}_{x}(\phi)\,\frac{\partial}{\partial
x}\,\Psi_{0}(x,p_{y}\,|\,\phi,\tau)\ ,$ (50)

Wave function $\Psi_{1}(x,p_{y}\,|\,\phi,\tau)$, as any other vector in three-
dimensional Hilbert space could be decomposed as a linear combination of three
orthogonal basis vectors $\big{|}\mbox{\boldmath$v$}_{1}\bigr{>}$,
$\big{|}\mbox{\boldmath$v$}_{2}\bigr{>}$ and
$\big{|}\mbox{\boldmath$v$}_{3}\bigr{>}$.

The first basis vector $\big{|}\mbox{\boldmath$v$}_{1}\bigr{>}$ is equivalent
to a part of $\Psi_{0}(\xi\,|\,\phi,\tau)$ in Eq. (48)

$\big{|}\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\bigr{>}=\left[\begin{array}[]{c}\cos\phi\,\,\Theta(x\,|\,\tau)\\\
-1\\\ \sin\phi\,\,\Theta^{\star}(x\,|\,\tau)\end{array}\right]\ .$ (51)

This means that
$\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\tau)\big{|}\mbox{\boldmath$v$}_{1}\bigr{>}=0$.
We also choose two remaining orthogonal vectors
$\big{|}\mbox{\boldmath$v$}_{2}(\xi\,|\,\phi,\tau)\bigr{>}$ and
$\big{|}\mbox{\boldmath$v$}_{3}(\xi\,|\,\phi,\tau)\bigr{>}$ in the following
way

$\displaystyle\big{|}\mbox{\boldmath$v$}_{2}(x\,|\,\phi,\tau)\bigr{>}$
$\displaystyle=$
$\displaystyle\left[\begin{array}[]{c}\cos\phi\,\,\Theta(x\,|\,\tau)\\\ +1\\\
\sin\phi\,\,\Theta^{\star}(x\,|\,\tau)\end{array}\right]\ ,$ (55)
$\displaystyle\big{|}\mbox{\boldmath$v$}_{3}(x\,|\,\phi,\tau)\bigr{>}$
$\displaystyle=$
$\displaystyle\left[\begin{array}[]{c}\sin\phi\,\,\Theta(x\,|\,\tau)\\\ 0\\\
-\cos\phi\,\,\Theta^{\star}(x\,|\,\tau)\end{array}\right]\ .$ (59)

Vector $\Psi_{1}(x,p_{y}\,|\,\phi,\tau)$ is now expanded in terms of the basis
$\big{|}\mbox{\boldmath$v$}_{1}\bigr{>}$,
$\big{|}\mbox{\boldmath$v$}_{2}\bigr{>}$ and
$\big{|}\mbox{\boldmath$v$}_{3}\bigr{>}$

$\Psi_{1}(x,p_{y}\,|\,\phi,\tau)=\varphi_{H}^{(1,1)}(x)\,\big{|}\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\bigr{>}+\varphi_{H}^{(1,2)}(x)\,\big{|}\mbox{\boldmath$v$}_{2}(x\,|\,\phi,\tau)\bigr{>}+\varphi_{H}^{(1,3)}(x)\,\big{|}\mbox{\boldmath$v$}_{3}(x\,|\,\phi,\tau)\bigr{>}\
.$ (60)

Now, substituting Eq. (60) into Eq. (50), we find

$\displaystyle\bigr{<}\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\,\big{|}\,\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\tau)\,\Psi_{1}(x,p_{y}\,|\,\phi,\tau)\bigr{>}$
(61) $\displaystyle=$
$\displaystyle\bigr{<}\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\,\big{|}\,\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\tau)\,\big{|}\left\\{\big{|}\,\varphi_{H}^{(0)}(x)\,\big{|}\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\bigr{>}+\varphi_{H}^{(1,2)}(x)\,\big{|}\mbox{\boldmath$v$}_{2}(x,|\,\phi,\tau)\bigr{>}+\varphi_{H}^{(1,3)}(x)\,\big{|}\mbox{\boldmath$v$}_{3}(x\,|\,\phi,\tau)\bigr{>}\right\\}\bigr{>}=0\
\ \ \ \ \ \ \ $

since transport operator $\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\tau)$ is
Hermitian.

The first order transport equation Eq. (79) becomes

$\bigr{<}\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\,\big{|}\,\hat{\Sigma}^{(3)}_{x}(\phi)\,\frac{\partial}{\partial
x}\,\Psi_{0}(x,p_{y}\,|\,\phi,\tau)\,\bigr{>}=\bigr{<}\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\,\big{|}\,\hat{\Sigma}^{(3)}_{x}(\phi)\,\frac{\partial}{\partial
x}\,\left\\{\,\big{|}\,\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\bigr{>}\,\varphi_{H}^{\,(0)}(x)\,\right\\}\,\bigr{>}=0\
$ (62)

could be written as

$\displaystyle\left[\Theta(x\,|\,\tau)+\Theta^{\star}(x\,|\,\tau)\right]\,\frac{\partial\,\varphi_{H}^{(0)}(x)}{\partial
x}+\Gamma_{\alpha}(x\,|\,\tau)\,\varphi_{H}^{(0)}(x)=0\,,$ (63)

where

$\Gamma_{\alpha}(x\,|\,\tau)=\frac{\partial\,\Theta(x\,|\,\tau)}{\partial
x}+\sin^{2}\phi\,\frac{\partial}{\partial
x}\left[\Theta^{\star}(x\,|\,\tau)-\Theta(x\,|\,\tau)\right]\,.$ (64)

The phase-dependent term could be simply rewritten through relative hopping
parameter $\alpha$ as $\sin^{2}\phi=(1+\alpha^{2})^{-1}$. Also, it worth
noting that for a dice lattice with $\phi=\pi/4$ function
$\Gamma_{\alpha=1}(x\,|\,\tau)=\partial/\partial
x[\Theta(x\,|\,\tau)+\Theta^{\star}(x\,|\,\tau)]$ becomes exactly symmetric
for $\Theta(x\,|\,\tau)\leftrightarrow\Theta^{\star}(x\,|\,\tau)$. For
$\phi=0$ we see the same equation as was derived for graphene in Ref.
[Zalipaev et al., 2015].

Eq. (63) has a general solution

$\varphi_{H}^{(0)}(\xi,\phi)=c_{0}\,\texttt{exp}\left\\{-\int\left[\Theta(x\,|\,\tau)+\Theta^{\star}(x\,|\,\tau)\right]^{-1}\Gamma_{\alpha}(\zeta\,|\,\tau)\,d\zeta\right\\}\
.$ (65)

Now we need to directly evaluate the following quantities:

$\displaystyle\Theta(x\,|\,\tau)=\frac{p_{\tau,-}}{\nu(x)}=-\frac{\tau
p_{x}-ip_{y}}{p}\,,$
$\displaystyle\Theta(x\,|\,\tau)^{\star}=\frac{p_{\tau,+}}{\nu(x)}=-\frac{\tau
p_{x}+ip_{y}}{p}\,,$
$\displaystyle\left[\Theta(x\,|\,\tau)+\Theta(x\,|\,\tau)^{\star}\right]^{-1}=-\frac{\tau\,p(x)}{2p_{x}(x)}\,,$
$\displaystyle\frac{\partial\,\Theta(x\,|\,\tau)}{\partial x}=-\frac{i\tau
p_{y}}{p_{\tau,+}^{3/2}\,p_{\tau,-}^{1/2}}\,\frac{\partial\,p_{x}(x)}{\partial
x}\,,$ $\displaystyle\frac{\partial\,\Theta(x\,|\,\tau)^{\star}}{\partial
x}=\frac{i\tau
p_{y}}{p_{\tau,+}^{1/2}\,p_{\tau,-}^{3/2}}\,\frac{\partial\,p_{x}(x)}{\partial
x}$ (66)

Using identities (66), equation (65) is evaluated as

$\varphi_{H}^{(0)}(\xi,\phi)=-\sqrt{\frac{p_{x}(x)+ip_{y}}{p_{x}(x)}}+\frac{i}{2}\sin^{2}\phi\,\tan^{-1}\left[-i\,\frac{p_{x}(x)}{p_{y}}\right]\,.$
(67)

We note that equation (67) is valid up to a normalization constant and up to a
phase factor. Thus, our result for gapless graphene with $\phi=0$ differs by a
complex phase $p_{x}(x)\pm ip_{y}$ compared to Ref. [Zalipaev et al., 2015].

For a dice lattice, we obtain

$\varphi_{H}^{(0)}(\xi)=\left\\{\frac{p^{2}_{x}(x)+p_{y}^{2}}{p^{2}_{x}(x)}\right\\}^{1/4}=\left\\{1+\left[\frac{p_{y}}{p_{x}(x)}\right]^{2}\,\right\\}^{1/4}\,.$
(68)

Now the leading order wave function $\Psi_{0}(\xi\,|\,\phi,\tau)$ in Eq. (48)
is determined completely, including the phase differences and the spatial
dependence of each its component. Apart from that, all the next orders of the
wave function expansion (29) could be now calculated, as we will do below for
the first order $\Psi_{1}(\xi\,|\,\phi,\tau)$.

We also note that $\varphi_{H}^{(0)}(\xi)$ in (68) is divergent
($\varphi_{H}^{(0)}(\xi)\longrightarrow\infty$) for $p_{x}(x)=0$ meaning that
WKB approximation is not valid around the turning points, denoted in Fig. 1 as
the points with zero longitudinal momentum. This situation was also earlier
found to be the case for graphene, Zalipaev et al. (2015) as well as a regular
Schrödinger electron.

###  III.3 General solution for the wave function

General transport equation is

$\frac{1}{\sqrt{2}}\,\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\,\tau)\,\Psi_{\lambda}(x,p_{y}\,|\,\phi,\tau)-\,\hat{\Sigma}^{(3)}_{x}(\phi)\,\frac{\partial}{\partial
x}\,\Psi_{\lambda-1}(x,p_{y}\,|\,\phi,\tau)=0\ ,$ (69)

Each of the wavefunctions $\Psi_{\lambda-1}(x,p_{y}\,|\,\phi,\tau)$ and
$\Psi_{\lambda}(x,p_{y}\,|\,\phi,\tau)$ is expended over an orthonormal set of
$\big{|}\mbox{\boldmath$v$}_{1}\bigr{>}$,
$\big{|}\mbox{\boldmath$v$}_{2}\bigr{>}$ and
$\big{|}\mbox{\boldmath$v$}_{3}\bigr{>}$ functions as:

$\displaystyle\Psi_{\lambda-1}(x,p_{y}\,|\,\phi,\tau)=\varphi_{H}^{(\lambda-1,1)}(x)\,\big{|}\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\bigr{>}+\varphi_{H}^{(\lambda-1,2)}(x)\,\big{|}\mbox{\boldmath$v$}_{2}(x\,|\,\phi,\tau)\bigr{>}+\varphi_{H}^{(\lambda-1,3)}(x)\,\big{|}\mbox{\boldmath$v$}_{3}(x\,|\,\phi,\tau)\bigr{>}\,,$
(70)
$\displaystyle\Psi_{\lambda}(x,p_{y}\,|\,\phi,\tau)=\varphi_{H}^{(\lambda,1)}(x)\,\big{|}\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\bigr{>}+\varphi_{H}^{(\lambda,2)}(x)\,\big{|}\mbox{\boldmath$v$}_{2}(x\,|\,\phi,\tau)\bigr{>}+\varphi_{H}^{(\lambda,3)}(x)\,\big{|}\mbox{\boldmath$v$}_{3}(x\,|\,\phi,\tau)\bigr{>}$

Once this is done, each side of Eq. (69) is multiplied by
$\big{\langle}\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\big{|}$,
$\big{\langle}\mbox{\boldmath$v$}_{2}(x\,|\,\phi,\tau)\big{|}$ and
$\big{\langle}|\mbox{\boldmath$v$}_{3}(x\,|\,\phi,\tau)\big{|}$, leading us to
the following three coupled differential equations

$\displaystyle\left[\left(\Theta+\Theta^{\star}\right)\,\frac{\partial\,\varphi_{H}^{(\lambda-1,1)}}{\partial
x}+\frac{\partial\,\Theta}{\partial
x}\,\left(\varphi_{H}^{(\lambda-1,1)}+\varphi_{H}^{(\lambda-1,2)}\right)+\left(\Theta-\Theta^{\star}\right)\frac{\partial\varphi_{H}^{(\lambda-1,2)}}{\partial
x}\right]\cos^{2}\phi-$ (71) $\displaystyle-$
$\displaystyle\left[\left(\frac{\partial\,\Theta^{\star}}{\partial
x}-\frac{\partial\,\Theta}{\partial
x}\right)\varphi_{H}^{(\lambda-1,3)}+\left(\Theta^{\star}-\Theta\right)\,\frac{\partial\varphi_{H}^{(\lambda-1,3)}}{\partial
x}\right]\,\frac{1}{2}\,\sin(2\phi)+$ $\displaystyle+$
$\displaystyle\left[\frac{\partial\,\Theta^{\star}}{\partial
x}\left(\varphi_{H}^{(\lambda-1,1)}+\varphi_{H}^{(\lambda-1,2)}\right)+\Theta\,\left(\frac{\partial\,\varphi_{H}^{(\lambda-1,1)}}{\partial
x}-\frac{\partial\,\varphi_{H}^{(\lambda-1,2)}}{\partial
x}\right)+\Theta^{\star}\,\left(\frac{\partial\,\varphi_{H}^{(\lambda-1,1)}}{\partial
x}+\frac{\partial\,\varphi_{H}^{(\lambda-1,1)}}{\partial
x}\right)\right]\,\cos^{2}\phi=0\,,$

obtained using $\big{\langle}\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\big{|}$,

$\displaystyle\left[\left(\Theta-\Theta^{\star}\right)\,\frac{\partial\,\varphi_{H}^{(\lambda-1,1)}}{\partial
x}+\frac{\partial\,\Theta}{\partial
x}\,\left(\varphi_{H}^{(\lambda-1,1)}+\varphi_{H}^{(\lambda-1,2)}\right)+\left(\Theta+\Theta^{\star}\right)\frac{\partial\varphi_{H}^{(\lambda-1,2)}}{\partial
x}\right]\cos^{2}\phi-$ (72) $\displaystyle-$
$\displaystyle\left[\left(\frac{\partial\,\Theta^{\star}}{\partial
x}-\frac{\partial\,\Theta}{\partial
x}\right)\varphi_{H}^{(\lambda-1,3)}+\left(\Theta^{\star}-\Theta\right)\,\frac{\partial\varphi_{H}^{(\lambda-1,3)}}{\partial
x}\right]\,\frac{1}{2}\,\sin(2\phi)+$ $\displaystyle+$
$\displaystyle\left[\frac{\partial\,\Theta^{\star}}{\partial
x}\left(\varphi_{H}^{(\lambda-1,1)}+\varphi_{H}^{(\lambda-1,2)}\right)+\Theta\,\left(\frac{\partial\,\varphi_{H}^{(\lambda-1,2)}}{\partial
x}-\frac{\partial\,\varphi_{H}^{(\lambda-1,1)}}{\partial
x}\right)+\Theta^{\star}\,\left(\frac{\partial\,\varphi_{H}^{(\lambda-1,1)}}{\partial
x}+\frac{\partial\,\varphi_{H}^{(\lambda-1,1)}}{\partial
x}\right)\right]\,\sin^{2}\phi=$ $\displaystyle=$ $\displaystyle
4\sqrt{2}\,\nu(x)\,\varphi_{H}^{(\lambda,2)}\,,$

when we multiplied by
$\big{\langle}\mbox{\boldmath$v$}_{2}(x\,|\,\phi,\tau)\big{|}$, and, finally,

$\displaystyle\frac{1}{2}\left(\Theta-\Theta^{\star}\right)\,\left[\frac{\partial\,\varphi_{H}^{(\lambda-1,1)}}{\partial
x}-\frac{\partial\,\varphi_{H}^{(\lambda-1,2)}}{\partial
x}\right]\,\sin(2\phi)=\sqrt{2}\,\nu(x)\,\varphi_{H}^{(\lambda,3)}\,.$ (73)

Here we also took into account the fact that
$\big{\langle}\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\big{|}$ is an
eigenfunction of $\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\,\tau)$ and

$\big{\langle}|\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\big{|}\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\,\tau)\big{|}\mbox{\boldmath$v$}_{2,3}(x\,|\,\phi,\tau)\bigr{>}=\big{\langle}\mbox{\boldmath$v$}_{2.3}(x\,|\,\phi,\tau)\big{|}\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\,\tau)\big{|}\mbox{\boldmath$v$}_{1}(x\,|\,\phi,\tau)\bigr{>}=0$
(74)

simply because transport operator
$\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\,\tau)$ is hermitian, i.e.,
$\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\,\tau)=\hat{\mathbb{O}}_{\,T}^{\dagger}(x,p_{y}\,|\,\phi,\,\tau)$.

We are mostly interested in the results for a dice lattice

$\displaystyle-\left(\Theta+\Theta^{\star}\right)\,\left[\frac{\partial\,\varphi_{H}^{(\lambda-1,+)}}{\partial
x}+\frac{\partial\,\varphi_{H}^{(\lambda-1,-)}}{\partial
x}\right]-\left[\frac{\partial\,\Theta}{\partial
x}^{\star}+\frac{\partial\,\Theta}{\partial
x}\right]\,\varphi_{H}^{(\lambda-1,+)}-$ (75)
$\displaystyle-\left(\Theta-\Theta^{\star}\right)\,\frac{\partial\,\varphi_{H}^{(\lambda-1,3)}}{\partial
x}+\left[\frac{\partial\,\Theta^{\star}}{\partial
x}-\frac{\partial\,\Theta}{\partial
x}\right]\,\varphi_{H}^{(\lambda-1,3)}=0\,,$

$\displaystyle\left(\Theta+\Theta^{\star}\right)\,\left[\frac{\partial\,\varphi_{H}^{(\lambda-1,+)}}{\partial
x}-\frac{\partial\,\varphi_{H}^{(\lambda-1,-)}}{\partial
x}\right]+\left[\frac{\partial\,\Theta}{\partial
x}^{\star}+\frac{\partial\,\Theta}{\partial
x}\right]\,\varphi_{H}^{(\lambda-1,+)}+$ (76)
$\displaystyle+\left(\Theta-\Theta^{\star}\right)\,\frac{\partial\,\varphi_{H}^{(\lambda-1,3)}}{\partial
x}-\left[\frac{\partial\,\Theta^{\star}}{\partial
x}-\frac{\partial\,\Theta}{\partial
x}\right]\,\varphi_{H}^{(\lambda-1,3)}=4\sqrt{2}\,\nu(x)\,\varphi_{H}^{(\lambda,2)}\,,$
(77)

$\left(\Theta-\Theta^{\star}\right)\,\frac{\partial\,\varphi_{H}^{(\lambda-1,-)}}{\partial
x}=\sqrt{2}\,\nu(x)\,\varphi_{H}^{(\lambda,3)}\,.$ (78)

Equations (75)-(78) is employed to obtain the three components of the each
consequent order $\varphi_{H}^{(\lambda,(1,2,3))}(x)$ from the previous terms
$\varphi_{H}^{(\lambda-1,(1,2,3))}(x)$ of expansion (29). Therefore, the
general problem of obtaining all the terms of the semiclassical wave function
is now solved, which is the key result of our work.

###  III.4 First-order wave function and applicability of WKB approximation

We apply first-order transport equation:

$\hat{\mathbb{O}}_{\,T}(\xi,p_{y}\,|\,\phi,\,\tau)\,\Psi_{1}(\xi,p_{y}\,|\,\phi,\tau)=\sqrt{2}\,\hat{\Sigma}^{(3)}_{x}(\phi)\,\frac{\partial}{\partial\xi}\,\Psi_{0}(\xi,p_{y}\,|\,\phi,\tau)\
,$ (79)

where $\Psi_{0}(\xi,p_{y}\,|\,\phi,\tau)=v_{1}\varphi_{H}^{(0,1)}$.

We obtain

$\left\\{\frac{\partial\,\Theta^{\star}(x)}{\partial
x}+\frac{\partial\,\Theta(x)}{\partial
x}\right\\}\,\varphi_{H}^{(0)}(x)+2\,\left[\Theta^{\star}(x)+\Theta(x)\right]\,\frac{\partial\,\varphi_{H}^{(0)}(x)}{\partial
x}=0$ (80)

which is equivalent to Eq. (63) and leads us to the zero order wavefunction
$\varphi_{H}^{(0)}(x)$.

$\nu(x)\,\varphi_{H}^{(1,2)}(x)=\frac{\sqrt{2}}{8}\,\left[\frac{\partial\,\Theta^{\star}(x)}{\partial
x}+\frac{\partial\,\Theta(x)}{\partial x}\right]\,\varphi_{H}^{(0)}(x)\,,$
(81)

as well as

$\nu(x)\,\varphi_{H}^{(1,3)}(x)=\frac{\sqrt{2}}{2}\,\left[\Theta^{\star}(x)-\Theta(x)\right]\frac{\partial\,\varphi_{H}^{(0)}(x)}{\partial
x}\,,$ (82)

and, finally, the only remaining and most important equation to find
$\varphi_{H}^{(1,1)}(x)$ from $\varphi_{H}^{(1,1)}(x)$ and
$\varphi_{H}^{(1,1)}(x)$ obtained from Eqs. (81) and (82), correspondingly:

$\displaystyle 2\Theta^{+}\,\frac{\partial\,\varphi_{H}^{(1,1)}}{\partial
x}+\frac{\partial\,\Theta^{+}}{\partial
x}\,\varphi_{H}^{(1,1)}+\frac{\partial\,\Theta^{+}}{\partial
x}\,\varphi_{H}^{(1,2)}-\frac{\partial\,\Theta^{-}}{\partial
x}\,\varphi_{H}^{(1,3)}-\Theta^{-}\,\frac{\partial\,\varphi_{H}^{(1,3)}}{\partial
x}=0\,,$ (83) $\displaystyle\Theta^{\pm}(x)=\Theta^{\star}(x)\pm\Theta(x)\,.$

Using Lagrange multiplier $\mu(x)=\texttt{exp}\big{[}\int f(x)dx\big{]}$, we
can immediately solve equation (83) for $\varphi_{H}^{(1,1)}(x)$ as

$\displaystyle\varphi_{H}^{(1,1)}(x)=\texttt{exp}\left[-\int
f(x)dx\right]\,\left\\{\text{const}-\int\texttt{exp}\left[\int
f(\xi)d\xi\right]\,g(x)dx\right\\}\,,$ (84) $\displaystyle
f(x)=\frac{1}{2\Theta^{+}(x)}\,\frac{\partial\,\Theta^{+}(x)}{\partial x}\,,$
$\displaystyle
g(x)=\frac{1}{2\Theta^{+}(x)}\,\left\\{\frac{\partial\,\Theta^{+}}{\partial
x}\,\varphi_{H}^{(1,2)}-\frac{\partial\,\Theta^{-}}{\partial
x}\,\varphi_{H}^{(1,3)}-\Theta^{-}\,\frac{\partial\,\varphi_{H}^{(1,3)}}{\partial
x}\right\\}\,.$

Equation (84) leads to the last remaining component of the first order wave
function $\varphi_{H}^{(1,1)}(x)$. Thus, we have developed a technique to
obtain the next order wave function, this concludes our effort to build WKB
theory for $\alpha-\mathcal{T}_{3}$ model.

Now, since all the wave function terms for the first order are obtained, it is
reasonable to verify the applicability of WKB approximation, i.e., the
magnitude of each consequent term of expansion (29) is much less compared to
the previous term, $\Psi_{1}(x)\ll\Psi_{0}(x)$ etc.

Let us use equation (81) connecting $\varphi_{H}^{(1,2)}(x)$ and
$\varphi_{H}^{(0)}(x)$ as an example with the most straightforward estimate
procedure. We also consider particle with a large longitudinal momentum
$p_{x}(x)\gg p_{y}$ and $p_{x}(x)\approx\nu(x)$ which for a Dirac cone
particle corresponds either to a high potential energy $V(x)$ or to a large
initial kinetic energy $E_{\gamma}$ of the incoming electron. In this case, we
can use the following approximations

$\displaystyle\left|\Theta^{\star}(x)+\Theta(x)\right|=\frac{2}{p_{x}(x)}{\nu(x)}\approx
1\,,$ (85) $\displaystyle\frac{\partial}{\partial
x}\,\left|\Theta^{\star}(x)+\Theta(x)\right|\approx\frac{p_{y}^{2}}{p_{x}(x)^{2}}\,\frac{\partial
p_{x}(x)}{\partial x}\,.$ (86)

Finally, the first order term wave function term $\varphi_{H}^{(1,2)}(x)$ is
estimated as

$\displaystyle\varphi_{H}^{(1,2)}(x)=\frac{\sqrt{2}}{8}\,\frac{1}{\nu(x)}\,\left\\{\frac{\partial\,\Theta^{\star}(x)}{\partial
x}+\frac{\partial\,\Theta(x)}{\partial
x}\right\\}\,\varphi_{H}^{(0)}(x)\approx$ (87)
$\displaystyle\approx\frac{1}{p_{x}(x)}\,\frac{p_{y}^{2}}{p_{x}(x)^{2}}\,\frac{\partial
p_{x}(x)}{\partial
x}\,\varphi_{H}^{(0)}(x)\approx\left[\frac{p_{y}}{p_{x}(x)^{2}}\right]^{2}\,\frac{\partial
p_{x}(x)}{\partial
x}\,\varphi_{H}^{(0)}(x)=\left[\frac{p_{y}}{p_{x}(x)}\right]^{2}\,\left|\frac{\partial\lambda(x)}{\partial
x}\right|\,\varphi_{H}^{(0)}(x)\,,$
$\displaystyle\left|\frac{\partial\lambda(x)}{\partial
x}\right|=\hbar\,\left[\frac{p_{y}}{p_{x}(x)}\right]^{2}\,\frac{\partial
p_{x}(x)}{\partial x}\,.$

Thus, we obtained a standard applicability criteria for a semiclassical WKB
approximation in quantum mechanics of a slow-changing wavelength of the
incoming particle and condition $p_{x}(x)\gg p_{y}$ discussed above.

Figure 2: (Color online) Classically forbidden regions (C.F.R.’s) as a
function of particle’s position $x$ and transverse momentum $k_{y}$ for an
electron in gapped $\alpha-\mathcal{T}_{3}$ materials with various bandgap
values moving under a linearly increasing potential $V(x)=V_{0}+ax$, where
$a=1.0$ for all plots and $V_{0}$ was chosen to satisfy
$\nu(x)=V(x)-E_{\gamma}=0$ for $x=0$ in each of the considered cases. The
values of the phase $\phi$ were taken $\phi=0$ (graphene), $\phi=\pi/6$ and
$\phi=\pi/4$ (a dice lattice), as labeled. Figure 3: (Color online)
Classically forbidden regions (C.F.R.’s) as a function of particle’s position
$x$ and phase $\phi$ for an electron in gapped $\alpha-\mathcal{T}_{3}$
materials with various bandgap values moving under a linearly increasing
potential $V(x)=V_{0}+ax$, where $a=1.0$ for all plots and $V_{0}$ was chosen
to satisfy $\nu(x)=V(x)-E_{\gamma}=0$ for $x=0$ in each of the considered
cases. The values of the transverse momentum $k_{y}$ were taken $k_{y}=0.1$
for panels $(a)$-$(c)$ and $k_{y}=0$ for plot $(d)$. Figure 4: (Color online)
Classically forbidden regions (C.F.R.’s) as a function of particle’s position
$x$ and transverse momentum $k_{y}$ for an electron in gapped
$\alpha-\mathcal{T}_{3}$ materials with various bandgap values moving under a
linearly increasing potential $V(x)=V_{0}+a\,\text{sign}(x)\,x^{2}$, where
$a=1.0$ for all plots and $V_{0}$ was chosen to satisfy
$\nu(x)=V(x)-E_{\gamma}=0$ for $x=0$ in each of the considered cases. The
values of the phase $\phi$ were taken $\phi=0$ (graphene), $\phi=\pi/6$ and
$\phi=\pi/4$ (a dice lattice), as labeled. Figure 5: (Color online)
Classically forbidden regions (C.F.R.’s) as a function of particle’s position
$x$ and phase $\phi$ for an electron in gapped $\alpha-\mathcal{T}_{3}$
materials with various bandgap values moving under a linearly increasing
potential $V(x)=V_{0}+a\,\text{sign}(x)\,x^{2}$, where $a=1.0$ for all plots
and $V_{0}$ was chosen to satisfy $\nu(x)=V(x)-E_{\gamma}=0$ for $x=0$ in each
of the considered cases. The values of the transverse momentum $k_{y}$ were
taken $k_{y}=0.1$ for panels $(a)$-$(c)$ and $k_{y}=0$ for plot $(d)$.

##  IV Applications of WKB approximation

In this Section, we will briefly comment on how we can use the obtained
approximated semiclassical wave function. The primary applications of a
semiclassical theory is finding the electron transmission for non-trivial
potential profiles where a direct calculation of the transmission amplitude is
either impossible or too complicated, as well as studying the resonant
tunneling or scattering with electron-hole transition.

If we only know action $S(x)$ which is directly related to the $x-$dependent
longitudinal momentum $p_{x}(x)=\partial S(x)/\partial x$, it is sufficient to
estimate the electron transmission through a barrier as an integral of
$p_{x}(x)$

$T(p_{y}\,|\,a,\Delta_{0},\phi)=\texttt{exp}\left[-\frac{2}{\hbar}\,\int\limits_{\rm
CFR}|p_{x}(x)|\,dx\right]$ (88)

over the so-called classically forbidden regions where $p_{x}(x)^{2}<0$, e.i.,
the longitudinal momentum is purely imaginary. Sonin (2009); Anwar et al.
(2020) Such states of the particle result in a significant drop of the
transmission coefficient, which is reflected in Eq. (88).

Let us consider a special case of an $\alpha-\mathcal{T}_{3}$ material with a
finite energy bandgap $\Delta_{0}$. Such a bandgap could be achieved either by
adding a dielectric substrate or irradiating our sample with an external off-
resonance dressing optical field. Gorbar et al. (2019); Iurov et al. (2019)
The Hamiltonian for a gapped $\alpha-\mathcal{T}_{3}$ lattice is obtained by
adding a $\phi$-dependent $\hat{\Sigma}_{z}(\phi)$ term

$\displaystyle\hat{\mathcal{H}}_{\Delta}(\phi)=\frac{\Delta}{2}\,\hat{\Sigma}_{z}(\phi)\,=\Delta\,\left[\begin{array}[]{ccc}\cos^{2}\phi&0&0\\\
0&-\cos 2\phi&0\\\ 0&0&-\sin^{2}\phi\end{array}\right]\ ,$ (92)
$\displaystyle\hat{\hat{\Sigma}}_{z}(\phi)=-i\,\left[\hat{\Sigma}_{x}(\phi),\hat{\Sigma}_{y}(\phi)\right]$

as it was used in Ref. [Biswas and Ghosh, 2018].

The transport operator $\hat{\mathbb{O}}_{\,T}(x,p_{y}\,|\,\phi,\tau)$ for
Hamiltonian (92) is modified as

$\displaystyle\hat{\mathbb{O}}_{\,T}(x,p_{y},\Delta_{0}\,|\,\phi,\tau)=\,\left[\begin{array}[]{ccc}\nu(x)+\Delta_{0}\cos^{2}\phi&\cos\phi\,\Xi(x)_{\mathbb{S},\tau}^{-}&0\\\
\cos\phi\,\Xi(x)_{\mathbb{S},\tau}^{+}&\nu(x)-\Delta_{0}\cos^{2}\phi&\sin\phi\,\Xi(x)_{\mathbb{S},\tau}^{-}\\\
0&\sin\phi\,\Xi(x)_{\mathbb{S},\tau}^{+}&\nu(x)-\Delta_{0}\sin^{2}\phi\end{array}\right]\
,$ (96)
$\displaystyle\Xi(x)_{\mathbb{S},\tau}^{\pm}=-i\hbar\,\tau\,\frac{\partial\mathbb{S}_{\Delta}(x)}{\partial
x}\pm ip_{y}$

where $\Delta_{0}\longrightarrow\Delta/V_{0}$. The longitudinal momentum
$p_{x}(\xi,\Delta_{0}\,|\,\phi)$ is now given by

$[p_{x}(x,\Delta_{0}\,|\,\phi)]^{2}=\nu^{2}(x)-p_{y}^{2}+\frac{\Delta_{0}^{2}}{8\,\nu(x)}\,\sin(2\phi)\,\sin(4\phi)-\frac{\Delta_{0}^{2}}{8}\left[5+3\cos(4\phi)\right]\,.$
(97)

The classically forbidden regions calculated by Eq. (97) for gapped
$\alpha-\mathcal{T}_{3}$ are presented in Fig. 2 \- 5 for a linear
$V(x)=V_{0}+a\,x$ and non-linear $V(x)=V_{0}+a\,\text{sign}(x)\,x^{2}$
potential profiles, correspondingly; $a$ is determined by the magnitude of the
electric field which created the potential barrier $V(x)$. We see a
substantial difference between the two types of potential which also strongly
affects the electron tunneling. We have also noted that the span of the
classically forbidden regions is generally decreased for a larger phase
$\phi$.

it is also important to realize that for $0<\phi<\pi/4$ there is a pole at
$\nu(x)=V(x)-E_{\gamma}=0$ and the sign of the term
$\Delta_{0}^{2}/[8\,\nu(x)]\,\sin(2\phi)\,\sin(4\phi)$ would change. This
means that $x=0$ must always be a boundary of a classically inaccessible
region, which we observe in all plots of Figs. 2 \- 5, except for graphene or
a dice lattice.

##  V Summary and Remarks

To summarize, this paper embodies our effort to build up a complete
semiclassical description for $\alpha-\mathcal{T}_{3}$ model which describes a
wide class of Dirac cone materials with an additional flat band in their
dispersions, also referred to as pseudospin-1 Dirac materials. Up to the
present time, there has been immense and strong experimental evidence for the
existence of a flat energy band in a large number of lab-fabricated materials.
Any of these materials fits the theoretical $\alpha-\mathcal{T}_{3}$ model
with a certain degree of precision so that our results are expected to be
applicable to a plenty of existing two-dimensional lattices.

Our development of Wentzel–Kramers–Brillouin (WKB) approximation includes
finding the semiclassical action specifically for $\alpha-\mathcal{T}_{3}$
model, as well as the $x$-dependent longitudinal momentum; calculating the
leading ($0^{th}$) order wave function, all of its components, their phase
differences and the spatial dependence. Most importantly, we have derived a
complete set of transport equations which connect any two consequent orders of
the semiclassical wave function ($\Psi_{\lambda}(x,p_{y}\,|\,\phi,\tau)$ and
$\Psi_{\lambda+1}(x,p_{y}\,|\,\phi,\tau)$). Therefore, the sought WKB wave
function could be now obtained up to any required order and precision.
Finally, we have solved these transport equations for the $1^{st}$ order wave
function and demonstrated the conditions for the applicability of
Wentzel–Kramers–Brillouin approximation for our model Hamiltonian, i.e.,
$\Psi_{1}(x,p_{y}\,|\,\phi,\tau)\ll\Psi_{0}(x,p_{y}\,|\,\phi,\tau)$.

Our derivations and the obtained results for $\alpha-\mathcal{T}_{3}$ appeared
to be drastically different from standard WKB approximation for a Schrödinger
particle studied in Quantum Mechanics textbooks, as well as the previously
obtained results for graphene. Zalipaev et al. (2015)

The obtained equations and the semiclassical electronic states could be
further used to investigate most of the crucial electronic properties of
$\alpha-\mathcal{T}_{3}$ materials. We have discussed the possible
applications of the obtained wave function, primarily, studying the tunneling
and transport properties of pseudospin-1 Dirac electrons in non-trivial
potential profiles. Specifically, we have considered a bandgap in
$\alpha-\mathcal{T}_{3}$ \- a situation for which the direct computation of
the wavefunction is complicated and is not even possible if the potential
becomes spatially non-uniform. The calculation of the electron transmission in
such cases is done by integrating the longitudinal electron momentum
$|p_{x}(x)|$ over the so-called classically forbidden regions with imaginary
momentum $p_{x}(x)$ ($p_{x}(x)^{2}<0$). We have calculated the transmission
for various kind of $\alpha-\mathcal{T}_{3}$ lattices and, specifically,
demonstrated that the transmission is mostly reduced for a larger phase
$\phi$.

In general, our developed semiclassical theory for a pseudospin-1 Dirac
electron could be employed to investigate resonant tunneling and scattering on
a potential barrier, as well as to recognize trapped or localized electronic
states in $\alpha-\mathcal{T}_{3}$.

We are confident that our work is crucial for investigating main electronic
properties of a whole class of innovative low-dimensional structures. Building
a Wentzel–Kramers–Brillouin approximation has been viewed as one of the key
results for each newly discovered material or a model Hamiltonian. The
obtained results will undoubtedly find their applications in the follow-up
research into electron tunneling and transport and electronic control in
$\alpha-\mathcal{T}_{3}$ materials, building devices and transistors based on
resonant tunneling or scattering of Dirac electrons through a custom-made
potential barrier profile. Our results are also important in terms of the
valleytronic applications since the low-energy states of a gapped material
directly depend on the valley index, as well as many other application of
these unique and innovative materials.

###### Acknowledgements.

A.I. would like to acknowledge the funding received from TRADA-52-113, PSC-
CUNY Award # 64076-00 52. D.H. was supported by the Air Force Office of
Scientific Research (AFOSR). G.G. would like to acknowledge Grant No.
FA9453-21-1-0046 from the Air Force Research Laboratory (AFRL).

## References

  * Vogl et al. (2017) M. Vogl, O. Pankratov, and S. Shallcross, Physical Review B 96, 035442 (2017).
  * Vandecasteele et al. (2010) N. Vandecasteele, A. Barreiro, M. Lazzeri, A. Bachtold, and F. Mauri, Physical Review B 82, 045416 (2010).
  * Zhang et al. (2012) Y. Zhang, Y. Barlas, and K. Yang, Physical Review B 85, 165423 (2012).
  * Weekes et al. (2021) N. Weekes, A. Iurov, L. Zhemchuzhna, G. Gumbs, and D. Huang, Physical Review B 103, 165429 (2021).
  * Zalipaev et al. (2015) V. Zalipaev, C. Linton, M. Croitoru, and A. Vagov, Physical Review B 91, 085405 (2015).
  * Zalipaev (2011) V. Zalipaev, Graphene - synthesis, characterization, properties and applications p. 81 (2011).
  * Berk et al. (1982) H. Berk, W. M. Nevins, and K. Roberts, Journal of Mathematical Physics 23, 988 (1982).
  * Sonin (2009) E. Sonin, Physical Review B 79, 195438 (2009).
  * Anwar et al. (2020) F. Anwar, A. Iurov, D. Huang, G. Gumbs, and A. Sharma, Physical Review B 101, 115424 (2020).
  * Ye et al. (2020) X. Ye, S.-S. Ke, X.-W. Du, Y. Guo, and H.-F. Lü, Journal of Low Temperature Physics pp. 1–12 (2020).
  * Xu and Lai (2019) H.-Y. Xu and Y.-C. Lai, Physical Review B 99, 235403 (2019).
  * Gangaraj et al. (2020) S. A. H. Gangaraj, C. Valagiannopoulos, and F. Monticone, Physical Review Research 2, 023180 (2020).
  * Bercioux et al. (2009) D. Bercioux, D. Urban, H. Grabert, and W. Häusler, Physical Review A 80, 063603 (2009).
  * Dóra et al. (2011) B. Dóra, J. Kailasvuori, and R. Moessner, Physical Review B 84, 195422 (2011).
  * Vidal et al. (1998) J. Vidal, R. Mosseri, and B. Douçot, Physical review letters 81, 5888 (1998).
  * Vidal et al. (2001) J. Vidal, P. Butaud, B. Douçot, and R. Mosseri, Physical Review B 64, 155306 (2001).
  * Illes (2017) E. Illes, Ph.D. thesis (2017).
  * Iurov et al. (2019) A. Iurov, G. Gumbs, and D. Huang, Physical Review B 99, 205135 (2019).
  * Dey and Ghosh (2018) B. Dey and T. K. Ghosh, Physical Review B 98, 075422 (2018).
  * Dey and Ghosh (2019) B. Dey and T. K. Ghosh, Physical Review B 99, 205429 (2019).
  * Goldman and Dalibard (2014) N. Goldman and J. Dalibard, Physical review X 4, 031027 (2014).
  * Iurov et al. (2013) A. Iurov, G. Gumbs, O. Roslyak, and D. Huang, Journal of Physics: Condensed Matter 25, 135502 (2013).
  * Kristinsson et al. (2016) K. Kristinsson, O. Kibis, S. Morina, and I. Shelykh, Scientific reports 6, 20082 (2016).
  * Kibis (2010) O. Kibis, Physical Review B 81, 165433 (2010).
  * Iurov et al. (2017) A. Iurov, L. Zhemchuzhna, G. Gumbs, and D. Huang, Journal of Applied Physics 122, 124301 (2017).
  * Sandoval-Santana et al. (2020) J. Sandoval-Santana, V. Ibarra-Sierra, A. Kunold, and G. G. Naumis, Journal of Applied Physics 127, 234301 (2020).
  * Kibis et al. (2017) O. Kibis, K. Dini, I. Iorsh, and I. Shelykh, Physical Review B 95, 125401 (2017).
  * Kibis (2022) O. Kibis, Physical Review A 105, 043106 (2022).
  * Cunha et al. (2022) S. Cunha, D. da Costa, J. M. Pereira Jr, R. Costa Filho, B. Van Duppen, and F. Peeters, Physical Review B 105, 165402 (2022).
  * Malcolm and Nicol (2016) J. Malcolm and E. Nicol, Physical Review B 93, 165433 (2016).
  * Iurov et al. (2020a) A. Iurov, G. Gumbs, and D. Huang, Journal of Physics: Condensed Matter 32, 415303 (2020a).
  * Abranyos et al. (2020) Y. Abranyos, O. L. Berman, and G. Gumbs, Physical Review B 102, 155408 (2020).
  * Islam and Dutta (2017) S. F. Islam and P. Dutta, Phys. Rev. B 96, 045418 (2017).
  * Huang et al. (2019) D. Huang, A. Iurov, H.-Y. Xu, Y.-C. Lai, and G. Gumbs, Physical Review B 99, 245412 (2019).
  * Cunha et al. (2021) S. Cunha, D. da Costa, J. M. Pereira Jr, R. Costa Filho, B. Van Duppen, and F. Peeters, Physical Review B 104, 115409 (2021).
  * Iurov et al. (2021) A. Iurov, L. Zhemchuzhna, G. Gumbs, D. Huang, P. Fekete, F. Anwar, D. Dahal, and N. Weekes, Scientific Reports 11, 1 (2021).
  * Urban et al. (2011) D. F. Urban, D. Bercioux, M. Wimmer, and W. Häusler, Physical Review B 84, 115136 (2011).
  * Illes and Nicol (2017) E. Illes and E. Nicol, Physical Review B 95, 235432 (2017).
  * Wang et al. (2020) J. Wang, J. Liu, and C. Ting, Physical Review B 101, 205420 (2020).
  * Louvet et al. (2015) T. Louvet, P. Delplace, A. A. Fedorenko, and D. Carpentier, Physical Review B 92, 155116 (2015).
  * Han and Lai (2022) C.-D. Han and Y.-C. Lai, Physical Review B 105, 155405 (2022).
  * Carbotte et al. (2019) J. Carbotte, K. Bryenton, and E. Nicol, Physical Review B 99, 115406 (2019).
  * Iurov et al. (2020b) A. Iurov, L. Zhemchuzhna, D. Dahal, G. Gumbs, and D. Huang, Physical Review B 101, 035129 (2020b).
  * Kovács et al. (2017) Á. D. Kovács, G. Dávid, B. Dóra, and J. Cserti, Physical Review B 95, 035414 (2017).
  * Illes et al. (2015) E. Illes, J. P. Carbotte, and E. J. Nicol, Phys. Rev. B 92, 245410 (2015).
  * Illes and Nicol (2016) E. Illes and E. Nicol, Physical Review B 94, 125435 (2016).
  * Raoux et al. (2014) A. Raoux, M. Morigi, J.-N. Fuchs, F. Piéchon, and G. Montambaux, Physical review letters 112, 026402 (2014).
  * Piéchon et al. (2015) F. Piéchon, J. Fuchs, A. Raoux, and G. Montambaux, in _Journal of Physics: Conference Series_ (IOP Publishing, 2015), vol. 603, p. 012001.
  * Biswas and Ghosh (2016) T. Biswas and T. K. Ghosh, Journal of Physics: Condensed Matter 28, 495302 (2016).
  * Biswas and Ghosh (2018) T. Biswas and T. K. Ghosh, Journal of Physics: Condensed Matter 30, 075301 (2018).
  * Balassis et al. (2020) A. Balassis, D. Dahal, G. Gumbs, A. Iurov, D. Huang, and O. Roslyak, Journal of Physics: Condensed Matter 32, 485301 (2020).
  * Neto et al. (2009) A. C. Neto, F. Guinea, N. M. Peres, K. S. Novoselov, and A. K. Geim, Reviews of modern physics 81, 109 (2009).
  * Tan et al. (2020) C.-Y. Tan, C.-X. Yan, Y.-H. Zhao, H. Guo, and H.-R. Chang, arXiv preprint arXiv:2012.11885 (2020).
  * Verma et al. (2017) S. Verma, A. Mawrie, and T. K. Ghosh, Physical Review B 96, 155418 (2017).
  * Champo and Naumis (2019) A. E. Champo and G. G. Naumis, Physical Review B 99, 035415 (2019).
  * Islam and Zyuzin (2019) S. F. Islam and A. Zyuzin, Physical Review B 100, 165302 (2019).
  * Islam and Saha (2018) S. F. Islam and A. Saha, Physical Review B 98, 235424 (2018).
  * Qiu et al. (2016) W.-X. Qiu, S. Li, J.-H. Gao, Y. Zhou, and F.-C. Zhang, Physical Review B 94, 241409 (2016).
  * Wang and Ran (2011) F. Wang and Y. Ran, Physical Review B 84, 241103 (2011).
  * Orlita et al. (2014) M. Orlita, D. Basko, M. Zholudev, F. Teppe, W. Knap, V. Gavrilenko, N. Mikhailov, S. Dvoretskii, P. Neugebauer, C. Faugeras, et al., Nature Physics 10, 233 (2014).
  * Ahmadkhani and Hosseini (2020) S. Ahmadkhani and M. V. Hosseini, Journal of Physics: Condensed Matter 32, 315504 (2020).
  * Jo et al. (2012) G.-B. Jo, J. Guzman, C. K. Thomas, P. Hosur, A. Vishwanath, and D. M. Stamper-Kurn, Physical review letters 108, 045305 (2012).
  * Ruostekoski (2009) J. Ruostekoski, Physical review letters 103, 080406 (2009).
  * Santos et al. (2004) L. Santos, M. Baranov, J. I. Cirac, H.-U. Everts, H. Fehrmann, and M. Lewenstein, Physical review letters 93, 030601 (2004).
  * Huang et al. (2011) X. Huang, Y. Lai, Z. H. Hang, H. Zheng, and C. Chan, Nature materials 10, 582 (2011).
  * Mukherjee et al. (2015) S. Mukherjee, A. Spracklen, D. Choudhury, N. Goldman, P. Öhberg, E. Andersson, and R. R. Thomson, Physical review letters 114, 245504 (2015).
  * Li et al. (2015) Y. Li, S. Kita, P. Muñoz, O. Reshef, D. I. Vulis, M. Yin, M. Lončar, and E. Mazur, Nature Photonics 9, 738 (2015).
  * Vicencio et al. (2015) R. A. Vicencio, C. Cantillano, L. Morales-Inostroza, B. Real, C. Mejía-Cortés, S. Weimann, A. Szameit, and M. I. Molina, Physical review letters 114, 245503 (2015).
  * Baba (2008) T. Baba, Nature photonics 2, 465 (2008).
  * Romhányi et al. (2015) J. Romhányi, K. Penc, and R. Ganesh, Nature communications 6, 6805 (2015).
  * Franchina Vergel et al. (2020) N. A. Franchina Vergel, L. C. Post, D. Sciacca, M. Berthe, F. Vaurette, Y. Lambert, D. Yarekha, D. Troadec, C. Coinon, G. Fleury, et al., Nano Letters (2020).
  * Leykam et al. (2018) D. Leykam, A. Andreanov, and S. Flach, Advances in Physics: X 3, 1473052 (2018).
  * Iurov et al. (2020c) A. Iurov, L. Zhemchuzhna, P. Fekete, G. Gumbs, and D. Huang, Physical Review Research 2, 043245 (2020c).
  * Li et al. (2017) Z. Li, T. Cao, M. Wu, and S. G. Louie, Nano letters 17, 2280 (2017).
  * Gorbar et al. (2019) E. Gorbar, V. Gusynin, and D. Oriekhov, Physical Review B 99, 155124 (2019).
  * Gorbar et al. (2021) E. V. Gorbar, V. P. Gusynin, and D. O. Oriekhov, Phys. Rev. B 103, 155155 (2021).