## 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}(\mathbf{k}\,|\,\Delta_{0})=\hbar v_{F}\hat{\mathcal{T}}_{0}^{(2)} \otimes\left(\mathbf{k}\cdot\hat{\mathbf{\Sigma}}^{(2)}\right)+\hbar v_{F}\,\Delta_{0} \left(\mathbf{k}\cdot\hat{\mathbf{\mathcal{T}}}^{(2)}\right)\,\otimes\hat{\mathbf{\Sigma}} _{0}^{(2)}\, \tag{1}\]  
where \(\mathbf{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 \(\mathbf{K}\) and \(\mathbf{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}\)\(\sim\)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:  
\[\hat{\mathbf{\Sigma}}_{x}^{(2)} = \hat{\mathcal{T}}_{x}^{(2)}=\left[ \begin{array}{cc}0&1\\ 1&0\end{array} \right]\,, \tag{2}\]  
\[\hat{\mathbf{\Sigma}}_{y}^{(2)} = \hat{\mathcal{T}}_{y}^{(2)}=\left[ \begin{array}{cc}0&-i\\ i&0\end{array} \right]\,, \tag{3}\]  
which act in the spin- and pseudo-spin subspaces, respectively. Correspondingly, \(\hat{\mathbf{\Sigma}}_{0}^{(2)}\) and \(\hat{\mathcal{T}}_{0}^{(2)}\) in Equation (1) are two \(2\times 2\) unit matrices in the spin- and pseudo-spin subspaces, written as:  
\[\hat{\mathbf{\Sigma}}_{0}^{(2)} = \hat{\mathcal{T}}_{0}^{(2)}=\left[ \begin{array}{cc}1&0\\ 0&1\end{array} \right]\,. \tag{4}\]  
Equation (1) is presented by using a Kronecker product of two \(2\times 2\) matrices, defined as:  
\[\overleftrightarrow{A}\otimes\overleftrightarrow{B}= \begin{bmatrix}a_{11} \overleftrightarrow{B}&a_{12}\overleftrightarrow{B}\\ a_{21}\overleftrightarrow{B}&a_{22}\overleftrightarrow{B}\end{bmatrix} \,. \tag{5}\]  
Figure 2: (**a**,**b**) Energy dispersions \(\varepsilon_{Y}(s_{1},s_{2}\,|\,k,\Delta_{0})\) for Kek-Y graphene represented by two nonequivalent Dirac cones with different Fermi velocities \(v_{F,1}\) and \(v_{F,2}\) corresponding to \(v_{F}(1\pm\Delta_{0})\). Here, \(s_{1}=\pm 1\) is a band index, while \(s_{2}=\pm 1\) corresponds to two different Fermi velocities due to strain coupling \(\Delta_{0}\). More specifically, for the case of Pauli matrices, we have:  
\[\hat{\Sigma}_{2}^{(2)}\otimes\hat{\mathcal{T}}_{2}^{(2)}=\left[ \begin{array}{ cccc}1&0&0&0\\ 0&-1&0&0\\ 0&0&-1&0\\ 0&0&0&1\end{array} \right], \tag{6}\]  
\[\hat{\Sigma}_{0}^{(2)}\otimes\hat{\mathcal{T}}_{2}^{(2)}=\left[ \begin{array}[] {cccc}1&0&0&0\\ 0&-1&0&0\\ 0&0&1&0\\ 0&0&0&-1\end{array} \right]. \tag{7}\]  
Therefore, any block-diagonal (b-d) matrix can be built by following this procedure, namely:  
\[\hat{\mathcal{M}}_{\text{b-d}}=\hat{\Sigma}_{0}^{(2)}\,\otimes\,\hat{ \mathcal{A}}=\left[ \begin{array}{c|cc}a_{11}&a_{12}&0&0\\ a_{21}&a_{22}&0&0\\ \hline 0&0&a_{11}&a_{12}\\ 0&0&a_{21}&a_{22}\end{array} \right]. \tag{8}\]  
As a result, Equation (1) can take a simple matrix form:  
\[\hat{\mathcal{H}}_{Y}(\mathbf{k}\,|\,\Delta_{0})=\hbar v_{F}\left[ \begin{array}{ c|cc}\mathbf{k}\cdot\hat{\Sigma}^{(2)}&\Delta_{0}\left(k_{x}-ik_{y}\right)\hat{ \Sigma}_{0}^{(2)}\\ \Delta_{0}\left(k_{x}+ik_{y}\right)\hat{\Sigma}_{0}^{(2)}&\mathbf{k}\cdot\hat{ \Sigma}^{(2)}\end{array}