## 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(x)\), including the special case with a square barrier \(V(x)=\Theta(x)\,\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}{ccc}V(x )&-\hbar\,\partial/\partial x-p_{y}&\Delta_{0}(-\hbar\,\partial/\partial x-p_ {y})&0\\ -\hbar\,\partial/\partial x-k_{y}&V(x)&0&\Delta_{0}(-\hbar\,\partial/\partial x -p_{y})\\ \hline\Delta_{0}(-\hbar\,\partial/\partial x+p_{y})&0&V(x)&-\hbar\,\partial/ \partial x+p_{y}\\ 0&\Delta_{0}(-\hbar\,\partial/\partial x+p_{y})&-\hbar\,\partial/\partial x-p_ {y}&V(x)\end{array} \right]\,, \tag{21}\]  
which is independent of \(y\), where \(p_{y}=\hbar k_{y}=\) const remains conserved. The key idea of the WKB approximation originates from a series expansion in powers of a small parameter \(\hbar\). Since we will consider a semi-classical particle with a large kinetic energy \(\varepsilon(k)\), 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\longrightarrow\frac{x}{W_{B}}\equiv\xi\quad\text{ and }\quad\varepsilon(k) \longrightarrow\frac{\varepsilon(k)}{V_{0}}\,, \tag{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(x)\longrightarrow\frac{V(x)}{V_{0}}\quad\text{ and }\quad p_{y} \longrightarrow\frac{v_{F}\,p_{y}}{V_{0}}\,. \tag{23}\]  
In this way, the spatial derivative now becomes:  
\[\partial/\partial x\to\partial/(W_{B}\,\partial\xi)\,. \tag{24}\]  
Such a scheme allows us to modify the Planck constant \(\hbar\), which is our chosen small parameter for a series expansion, leading to:  
\[\hbar\to\left(\frac{v_{F}}{W_{B}V_{0}}\right)\,\hbar\, \tag{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:  
\[\Psi_{Y}(s_{1},s_{2}\mid x,k_{y}) = \exp\left[\frac{i}{\hbar_{0}}\,\mathbb{S}(x)\right]\,\sum_{\lambda =0}^{\infty}\,(-i\hbar_{0})^{\lambda}\,\Psi_{\lambda}(x,k_{y})\] \[= \exp\left[\frac{i}{\hbar_{0}}\,\mathbb{S}(x)\right]\,\left[\Psi_ {0}(x,k_{y})-i\hbar_{0}\,\Psi_{1}(x,k_{y})-\hbar_{0}^{2}\,\Psi_{2}(x,k_{y})+ \cdots\right]\,,\]where \(\mathbb{S}(x)\) 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:  
\[\hat{\mathbb{D}}_{x}^{(4)}(\Delta_{0})\left[\,\frac{\partial}{ \partial x}\,\Psi_{\lambda}(s_{1},s_{2}\,|\,x,k_{y})\,\right]-\hat{\mathbb{O}} _{\,T}(s_{1},s_{2}\,|\,x,p_{y},\Delta_{0})\,\Psi_{\lambda+1}(s_{1},s_{2}\,|\,x,k_{y})=0\, \tag{27}\]  
where  
\[\hat{\mathbb{D}}_{x}^{(4)}(\Delta_{0}) = \hat{\mathcal{T}}_{0}^{(2)}\otimes\hat{\Sigma}_{x}^{(2)}+\Delta_{ 0}\,\hat{\mathcal{T}}_{x}^{(2)}\otimes\hat{\Sigma}_{0}^{(2)}\] \[= \left[ \begin{array}{c|c|c}0&1&\Delta_{0}&0\\ 1&0&0&\Delta_{0}\\ \hline\Delta_{0}&0&0&1\\ 0&\Delta_{0}&1&0\end{array}