\section{Analytic GBZ results for generic 1D systems with two terms}\label{sec1}
As a foundation for the analysis of the effective BZ of higher-dimensional lattices, we develop in this section a general analytic derivation for the generalized Brillouin zone (GBZ) results for 1D systems with 2 hopping terms, and compare them with numerics. In this supplement, we shall refer to ``GBZ'' and ``effective BZ'' interchangeably. Note that these results may no longer hold when there are more than 2 hopping terms, which is the more interesting scenario which this work is concerned with. \\

\subsection{Preliminaries}

Consider the following 1D single-band model with dissimilar left/right hoppings:
\begin{equation}\label{eqS1}
\begin{split}
H_\text{1D}=\sum_{n} A|n\rangle\langle n+\alpha|+ B|n\rangle\langle n-\beta|
\ ,
\end{split}
\end{equation}
whose eigenvalues are $E(z)=Az^{\alpha}+Bz^{-\beta}$, where $z = \text{e}^{ik}$ describes the momentum component normal to the boundary of interest. This effective 1D single-band system has only 2 hoppings,
one is to move $\alpha$ sites to the left by  amplitude $A$, the other with amplitude $B$,  $\beta$ sites to the right. In higher-dimensional contexts, $A$ and $B$ can depend on the momenta from the other directions.

Corresponding to this eigenvalue $E(z)=Az^{\alpha}+Bz^{-\beta}$ is the eigenstate $|\psi(z)\rangle$. In deference to Bloch's theorem, we assume that it takes the form of a \emph{generalized} Bloch state given by the ansatz $|\psi(z)\rangle=(...,\psi_n(z),...)^{\text{T}}$, $\psi_n(z)=z^{n}$, position index $n$, which satisfies the bulk equation
\begin{equation}\label{eqS2}
A\psi_{n+\alpha}(z)+B\psi_{n-\beta}(z)=E(z)\psi_n(z)\ .
\end{equation}
At a particular energy $E(z)$,  there exists other wavefunction solutions $|\psi(z')\rangle$ which satisfies the same eigenvalue $E(z)$, i.e.
\begin{equation}\label{z'}
\begin{split}
&Az'^{\alpha}+Bz'^{-\beta}-E(z)=0\ .
\end{split}
\end{equation}
This equation above is a polynomial relation of order $\alpha+\beta$ in $z$, and it's easy to get its solutions $z_1,z_2,...,z_{\alpha+\beta}$, which can all be expressed in terms of $z$. These solutions shall provide information about how $z$ is controlled by hopping amplitudes $A$ and $B$.

Hereupon, an eigensolution $ |\Psi \rangle $ with eigenenergy $E(z)$ can be written as
 \begin{equation}
|\Psi \rangle=\sum^{\alpha+\beta}_{i=1}c_i|\psi(z_i)\rangle\  ,
\end{equation}
 with coefficients $c_i$ $(i=1,2,...,\alpha+\beta)$. To determine how they are constrained, we apply open boundary conditions onto $|\Psi \rangle$, arriving at
   \begin{equation}
\Psi_{n'}=\sum^{\alpha+\beta}_{i=1}c_i z_i^{n'}=0\ ,\quad\quad n'=0,-1,-2,...1-\beta\quad \text{or}\quad L+1,L+2,...L+\alpha\ ,
\end{equation}
that is,  there have $\alpha+\beta$ constraints from the open boundary conditions (OBCs), which together combine to form the GBZ characteristic equation
  \begin{equation}\label{eqS2}
\text{det}\ M=\sum^{\alpha+\beta}_{n_1,n_2,...n_{\alpha+\beta}=1}\varepsilon_{n_1,n_2,...,n_{\alpha+\beta}} z_{n_1}^{1-\beta}\times z_{n_2}^{2-\beta}\times...\times z_{n_{\beta}}^{0}\times z_{n_{\beta+1}}^{L+1} \times z_{n_{\beta+2}}^{L+2}\times...\times z_{n_{\alpha+\beta}}^{L+\alpha}=0\ ,
\end{equation}
where the matrix $M$ is $(\alpha+\beta)\times(\alpha+\beta)$ square array , $\varepsilon_{n_1,n_2,...,n_{\alpha+\beta}}$ is $\alpha+\beta$-order arrangement, and $\varepsilon_{n_1,n_2,...,n_{\alpha+\beta}}=-1$, when  $n_1,n_2,...,n_{\alpha+\beta}$ is odd arrangement,   $\varepsilon_{n_1,n_2,...,n_{\alpha+\beta}}=1$ if  $n_1,n_2,...,n_{\alpha+\beta}$ is even arrangement. Since $z_1,z_2,...,z_{\alpha+\beta}$ are all functions of $z$,  the GBZ characteristic equation Eq.\eqref{eqS2} is in actuality only a function of $z$.

 \begin{figure}[t]
    \begin{centering}
    \includegraphics[width=.65\linewidth]{FIGS1-1.pdf}
    \par\end{centering}
    \protect\caption{
	(a1--c1) Solutions $z$ to the characteristic equation (Eq.\eqref{eqS2}) for a few illustrative values of $\alpha,\beta$, with $A=B=1$ and system size of $L=10$ sites. Blue dots represent the full set of solutions, but only those lying on the black constraint curves constitute the actual GBZ solutions for OBCs. (a2--c2)  Their corresponding eigenenergies $E(z)=Az^{\alpha}+bz^{-\beta}$ (black) are the set of $z$ points belonging to the GBZ, which agree excellently with $E_{\text{OBC}}$ (red) of Eq.\eqref{eqS1}. (a3--c3) shows the excellent agreement between numerical eigenstate solutions (red) and their corresponding GBZ solutions (black).}
	\label{fig:S1}
\end{figure}

Taking the model Eq.\eqref{eqS1} with $\alpha=2,\beta=1$ as a sample,  the eigenvalues $E(z)=Az^{2}+B/z$ can be obtained by bulk equation Eq.\eqref{eqS2}. Considering the same eigenenergy $E(z)$,  the eigensolution $ |\Psi \rangle $ can be written as $ |\Psi \rangle =\sum_{i=1}^{3}c_i|\psi(z_i)\rangle $, $\psi_n(z_i)=z_i^n$, and  $z_i$  which satisfy $Az_i^2+B/z_i-E(z)=0$, are constrained by $z$ and hopping amplitude $A,B$
 \begin{equation}\label{eqS7}
z_1=z\ ,\  z_2=-\frac{z}{2}+\sqrt{\frac{B}{Az}+\frac{z^2}{4}}\ ,\ z_3=-\frac{z}{2}-\sqrt{\frac{B}{Az}+\frac{z^2}{4}}\ .
\end{equation}
 Applying open boundary conditions into $|\Psi \rangle$, we can get GBZ characteristic equation
 \begin{equation}\label{EqS8}
 \begin{split}
\text{det}\ M= \text{det}\begin{pmatrix} 1&1&1 \\ z_1^{L+1}&z_2^{L+1}&z_3^{L+1}\\ z_1^{L+2}&z_2^{L+2}&z_3^{L+2}\end{pmatrix}=
(z_1z_2)^{L+1}(z_2-z_1)+(z_1z_3)^{L+1}(z_1-z_3)+(z_2z_3)^{L+1}(z_3-z_2)=0\ ,
\end{split}
\end{equation} which is only a function of $z$. The determinant must vanish so that we have a nontrivial eigenstate solution $(c_1,c_2,c_3)^{\text{T}}$ of the matrix $M$ with eigenvalue $0$.
It is worth noting that not all the  solution $z$ of  GBZ characteristic equation Eq.\eqref{EqS8}  is  GBZ results of the system, as show in FIG.~\ref{fig:S1}(a1---a3).
To put it simply, for fixed $z$, $z_1, z_2, z_3$ which are from Eq.\eqref{eqS7}
can be rearranged by $|z_1|\leq|z_2|\leq|z_3|$. If the absolute value $|z_3|$ is not equal to $|z_2|$, we have $\lim\limits_{L\gg1} c_3=0,  (c_1=-c_2)$, that is, the solution $z_3$ does not belong to GBZ. And in other cases, we  can also select GBZ solutions from the set of solution $z$ by the values of the coefficients  $c_1,c_2,c_3$.

\subsection{Determination of the GBZ for OBCs}

However, not all solutions $z$ of the characteristic equation Eq.\eqref{eqS2} contribute to the actual OBC solutions. Those that do define the GBZ. In Figs.~\ref{fig:S1}(a1---c1), we see that of all solutions $z$ (blue dots), only those that lie on the black curves correspond to values of $E(z)$ that coincide with actual numerically obtained eigenvalues (Figs.~\ref{fig:S1}(a2---c2)). \\

From known results on GBZ construction~\cite{Lee2019anatomy,yang2019auxiliary,lee2020unraveling}, solutions $z$ that belong to the GBZ (i.e. black curves) are those $z$ such that there exists another $z'\neq z$ such that $|z|=|z'|$ and $E(z)=E(z')$. If more than one such pair of solutions exist, the GBZ is defined by the pair with $|z|=|z'|$ closest to unity. This is because $|z|$ determines the spatial decay rate of the eigensolution, and the GBZ pair is given by the pair of solutions that are mostly slowly decaying, and yet able to satisfy OBCs at both boundaries (which can be arbitrarily far) simultaneously. In other words, the GBZ is obtained through the pair among
$z_1,z_2,...,z_{\alpha+\beta}$ with the same absolute value.
\\

Let's now solve for the GBZ of our 1D lattice $H_\text{1D}$, with $z$ parametrized by a wavenumber $k$. First, we write the two degenerate solutions of $z$ as $z_{\pm}=z_0\text{e}^{ik'\pm k }$, with real $k', k$ and possessing the same energy $E(z_\pm)$. Here $z_0$ is a parameter that is unspecified for now. Taking $z_{\pm}$ into the energy equation,
\begin{equation}
\begin{split}
E=E(z_{+})=E(z_{-})=Az_0^{\alpha}\text{e}^{i\alpha (k'+k)}+Bz_0^{-\beta}\text{e}^{i\beta (k'+k)}=Az_0^{\alpha}\text{e}^{i\alpha (k'-k)}+Bz_0^{-\beta}\text{e}^{-i\beta (k'-k)}\ .
\label{Epm}
\end{split}
\end{equation}
As such, we can solve for $z_{\pm}$:
\begin{equation}\label{aaa}
\begin{split}
z_{\pm}&=\left(\frac{B\sin  \beta k}{A\sin \alpha k}\right)^{\frac{1}{\alpha+\beta}}\text{e}^{\pm i k+i\frac{2\pi v}{\alpha+\beta}} \ ,\\
\end{split}
\end{equation}
with $v=1,2,...,\alpha+\beta$. Note that $|z_\pm|\neq 1$ in general, and only when $k=n\pi/(\alpha+\beta)$ do we have $|\sin  \beta k|/|\sin  \alpha k|=1$.
Due to the periodicity of $2\pi/(\alpha+\beta)$, without loss of generality, we shall from now on set the range of parameter $k$ to be $[n\pi-\pi/(\alpha+\beta), n\pi+\pi/(\alpha+\beta))$, with integer $n$. For definiteness, we take $n=0$.\\

In all, the GBZ of our model $H_\text{1D}$ (Eq.\eqref{eqS1}) is given by
\begin{equation}\label{eqS3}
\text{GBZ}_\text{1D}=\left \{z=\left(\left.\frac{B\sin  \beta k}{A\sin \alpha k}\right)^{\frac{1}{\alpha+\beta}}\text{e}^{ i k+i\frac{2\pi v}{\alpha+\beta}}\,\right| k=\left(-\frac{\pi}{\alpha+\beta},\frac{\pi}{\alpha+\beta}\right],v=1,2,..\alpha+\beta \right\}\ ,
\end{equation}
giving energies
\begin{equation}\label{eqS4}
\begin{split}
\bar E(k) =E(z)&=Az^{\alpha}+Bz^{-\beta} =\left(\frac{A^{\beta}B^{\alpha}}{(\sin (\alpha k))^{\alpha}(\sin (\beta k))^{\beta}}\right)^{\frac{1}{\alpha+\beta}} \sin \left((\alpha+\beta)k\right) \text{e}^{\frac{2i\pi \alpha v}{\alpha+\beta}}\ ,
\end{split}
\end{equation}
which agrees closely with numerical OBC eigenenergies in Figs.~\ref{fig:S1}(a2---c2). Here, the bar above $\bar E(k)$ indicates that the energy function is evaluated on the GBZ, i.e. it depends on $k$ not on the ordinary BZ, but through the GBZ. Previously, it was already conceived that complex momentum can be used to represent state decay~\cite{kohn1959analytic,he2001exponential,lee2015free}, but in this GBZ formalism, the imaginary part of the momentum is specifically solved to give the profile of OBC eigenstates.
Eq.~\ref{eqS3} is a main result that will be used in the GBZ computation of higher-dimensional cases later on.
As consistent with the fact that OBC spectra cannot undergo any further NHSE (in the same direction), the spectral winding of $E(z)$ above is always zero, that is, the distribution of $E(z)$ in the complex plane cannot form closed curves.\\

In generic 1D lattices with multiple hopping terms, Eq.~\ref{Epm} will have to be replaced by a simultaneously polynomial relations which has to be solved numerically. The resultant GBZ is still well-defined, although it is likely obtainable only numerically.

\section{2D non-Hermitian lattices with 2D GBZs}\label{sec2}
This section introduces how the GBZ can be obtained for 2D lattices through a few concrete examples with different hopping terms.
\subsection{2D lattice model with 2 hopping terms}
To connect with our previous treatment of 1D lattices with 2 hopping terms, we first consider the simplest case of 2D lattices with 2 hopping terms:
\begin{equation}\label{eqS21}
H_\text{2D,2}= \sum_{m,n}t|m,n\rangle\langle m+\alpha,n+a|+t'|m,n\rangle\langle m-\beta,n-b|\ ,
\end{equation}
with hoppings of amplitudes $t$ (or $t'$) corresponding to transitions of $(-\alpha,-a)$ sites ( or $(\beta,b)$ sites) on the lattice.\\

We consider the ansatz $|\psi\rangle\propto\sum_{m,n}\psi_{mn}|m,n\rangle$ with $\psi_{mn}=z_x^mz_y^n$.
In the bulk, it gives eigenenergies $E(z_x,z_y)$
  \begin{equation}
E(z_x,z_y)=t z_x^{\alpha}z_y^a+t' z_x^{-\beta}z_y^{-b}.
\label{Ezxzy}
\end{equation}
For a fixed $z_x$, there are $a+b$ solutions $z_{y,i}$, ($i=1,2,...,a+b$)  corresponding to the same energy $E(z_x,z_y)$. Similarly, we can also find $\alpha+\beta$ solutions $z_{x,j}$, ($j=1,2,...,\alpha+\beta$) corresponding to fixed $E(z_x,z_y)$ and $z_y$. Both sets of solutions $z_{y,i}$, ($i=1,2,...,a+b$) and $z_{x,j}$, ($j=1,2,...,\alpha+\beta$) can be written entirely in terms of $z_x$ and $z_y$.\\

We next show a detailed for derivation of the GBZ of $H_\text{2D,2}$. For a quick heuristic approach, the reader may directly skip to Eqs.~\ref{eqS22} and~\ref{eqS23}.

Implementing open boundary conditions on the wavefunction $|\Psi\rangle$ for a fixed eigenenergy $E(z_x,z_y)$ gives
    \begin{equation}\label{eqS21.5}
    \begin{split}
|\Psi\rangle=|\Psi^1\rangle=|\Psi^2\rangle\ \propto \sum_{m,n}\Psi_{mn}|m,n\rangle=\sum_{m,n}\Psi^1_{mn}|m,n\rangle=\sum_{m,n}\Psi^2_{mn}|m,n\rangle\ ,\\
\quad \Psi^1_{mn}=\sum_{z_y} \left( \left(\sum_{j}^{\alpha+\beta}f_{z_{x,j}}z_{x,j}^m\right)g_{z_y} z_y^n\right)\ , \quad \Psi^2_{mn}=\sum_{z_x}\left( f_{z_x} z_x^m \left(\sum_{i}^{a+b}{g}_{z_{y,i}}z_{y,i}^n\right)\right)\ ,\\
\end{split}
\end{equation}
Here, we have written $|\Psi\rangle$ either as $|\Psi^1\rangle$ or $|\Psi^2\rangle$, depending on the order of expansion in terms of the $x$ and $y$ coefficients.
These two ways to expand are of course equivalent, and note the relation $f_{z_x}g_{z_{y,i}}=g_{z_{y,i}}f_{z_x}$. The $\sum_i,\sum_j$ summations in $\Psi^{1}_{mn},\Psi^{2}_{mn}$ refers to sums over all sets $z_x,z_y$ which have same energy $E$. To treat boundary conditions along different direction (i.e.,${\bm x, \bm y}$ directions), we can choose the appropriate expansion and proceed like in the 1D case (see Sec.~\ref{sec1}), paying careful attention to the indices i.e. for OBCs along the ${\bm x}$ direction, we have  $\Psi_{m'{n}}=\Psi^1_{m'{n}}=0 $ for any $n$ and $m'=-\beta+1,-\beta+2,...0,L+1,L+2,..L+\alpha $ (assuming $\alpha,\beta>0$).
Thus we obtai
\begin{equation}\label{eqS22}
\begin{split}
z_x^{\alpha+\beta}=\frac{t'}{t }\frac{z_1^{\alpha+\beta}}{z_y^{a+b}}, \quad z_1^{\alpha+\beta}=\frac{\sin \beta k_1}{\sin \alpha k_1}\text{e}^{i(\alpha+\beta)k_1}\ ,
\end{split}
\end{equation}
with momentum parameter $k_1\in  \left(-\frac{\pi}{\alpha+\beta},\ \frac{\pi}{\alpha+\beta}\right]$ as in a quasi-1D case (see Sec.~\ref{sec1}. for details). Likewise, considering the expansion via $|\Psi\rangle=|\Psi^2\rangle$  with boundary conditions along the ${\bm y}$ direction, we obtain the alternative expressions
\begin{equation}\label{eqS23}
\begin{split}
z_y^{a+b}=\frac{t'}{t }\frac{z_2^{a+b}}{z_x^{\alpha+\beta}}, \quad z_2^{a+b}=\frac{\sin b k_2}{\sin a k_2} \text{e}^{i(a+b)k_2}\ ,
\end{split}
\end{equation}
with $k_2\in  \left(-\frac{\pi}{a+b},\ \frac{\pi}{a+b}\right]$. To solve for the GBZs of simultaneous x-OBCs and y-OBCs in 2D, we simultaneously solve Eqs.~\ref{eqS22} and~\ref{eqS23}. Note that heuristically, Eqs.~\ref{eqS22} and~\ref{eqS23} can be simply written down by treating $H_\text{2D,2}$ as a 1D chain with hoppings dependent on the transverse momentum. For instance, we can write Eq.~\ref{Ezxzy} either as $E(z_x)=(tz_y^a)z_x^\alpha+(tz_y^{-b})z_x^{-\beta}$ or $E(z_y)=(tz_x^\alpha)z_y^a+(tz_x^{-\beta})z_y^{-b}$, and apply Eq.~\ref{eqS3} to obtain Eqs.~\ref{eqS22} and~\ref{eqS23} respectively. Importantly, we shall see that this approach fails in general 2D lattices, even though it is valid in this case where there are only two hopping terms from each site.
\begin{figure}[h]
\centering
\includegraphics[width=.75 \linewidth]{FIGS1-2.pdf}
\caption{Spectra and corner-localized eigenstates for $H_\text{2D,2}$ in the $\alpha/\beta\neq a/b$ case. (a--c) Exactly zero (flatband) OBC eigenenergies $E_{\text{OBC}}$ due to non-Bloch collapse, which coincide with $E(z_x,z_y)$ from Eq.\eqref{eqS25}. (a1---c1) Corresponding illustrative eigenstates are perfectly localized at the boundary dictated by the direction of localization. The model parameters are $L_x=L_y=15$, $t=2,\quad t'=1$  and the specific values of $\alpha,\beta,a,b$ are indicated in the figure.}\label{fig:S22}
\end{figure}

In the following, we shall specialize to the case where the two hoppings are pointing in the same direction i.e. $\alpha /\beta=a/b$, for reasons that will be become evident soon. In this case, the GBZ and energy satisfy
\begin{equation}\label{eqS24}
\begin{split}
&\text{GBZ}_\text{2D,2}=\left\{\left.z_x^{\alpha+\beta}z_y^{a+b}=\frac{t'\sin b k}{t\sin a k} \text{e}^{i(a+b)k}\right|k=\left(-\frac{\pi}{a+b},\frac{\pi}{a+b}\right]\right\}\ ,\\
&\bar E(k)=E(z_x,z_y)=\left(\frac{t^{b}t'^{a}}{(\sin (a k))^{a}(\sin (b k))^{b}}\right)^{\frac{1}{a+b}} \sin \left((a+b)k\right) \text{e}^{\frac{2i\pi  av}{a+b}}\ ,
\end{split}
\end{equation}
with $k=\left(-\frac{\pi}{a+b},\frac{\pi}{a+b}\right]$ and $v=1,2,...,a+b$, as plotted in FIG.~\ref{fig:S21} (a---c). As this case only contains the combination $z_x^{\alpha+\beta}z_y^{a+b}$,  rather than $z_x^{\alpha+\beta}$ and $z_y^{a+b}$, it is essentially a 1D model along the $\alpha\bm x+a\bm y$ direction, which is consistent with the results of numerical diagonalization as shown in FIG.~\ref{fig:S21}(a1---a3). \\
\begin{figure}[h]
\centering
\includegraphics[width=.75 \linewidth]{FIGS1-3.pdf}
\caption{Spectra and illustrative eigenstates for $H_\text{2D,2}$ in the $\alpha/\beta=a/b$ case. (a--c) Excellent agreement between the GBZ eigenenergies $E(z_x,z_y)$ from Eq.\eqref{eqS24} and OBC eigenenergies $E_{\text{OBC}}$ of Eq.\eqref{eqS21}. (d--f) Spatial profiles for the $E_{\text{OBC}}=0$ eigenstate of a few other illustrative cases, clearly showing that the skin states are aligned along the $\alpha \hat x +a \hat y$ (or $\beta \hat x +b \hat y$) direction. The model parameters are $L_x=L_y=15$, $t=1,\  t'=2$  and the specific values of $\alpha,\beta,a,b$ are indicated in the figure.}\label{fig:S21}
\end{figure}

We next discuss the other case with $\alpha /\beta\neq a/b$, where $|z_1|$ and $|z_2|$ can only coincide at $z_1^{\alpha+\beta}=z_2^{a+b}=-1$. Hereby, its GBZ and energy are simply given by
\begin{equation}\label{eqS25}
\begin{split}
&\text{GBZ}_\text{2D,2}=\left\{z_x^{\alpha+\beta}z_y^{a+b}=-\frac{t'}{t}\right\}\ ,\\
&\qquad\qquad \quad  E(z_x,z_y)=0\ ,
\end{split}
\end{equation}
As confirmed in FIG.~\ref{fig:S22}(a--c), the eigenenergies are indeed zero, and the eigenstates (FIG.~\ref{fig:S22}(a1--c1)) are perfectly localized at a corner determined by $\alpha/\beta$, $a/b$ and $t/t'$. This is because of uncompensated unbalanced hoppings along at least one direction, which leads to non-Bloch collapse~\cite{longhi2020non,crespi2013dynamic}. Since there is no nontrivial dynamics to speak of in this case, we shall not discuss it further.

\begin{figure}[h]
\centering
\includegraphics[width=.65 \linewidth]{FIGS1-4.pdf}
\caption{The three directed hopping configurations in each of our type I,II and III models. For type I, two hoppings are orthogonal; for type II, two hoppings have common $\beta$ displacements in the direction antiparallel to the third hopping; for type III, there exists two common $\beta$ component displacements orthogonal to the two common $a$ displacements. }\label{fig:S221}
\end{figure}

\subsection{2D lattice model with 3 hopping terms}

We next discuss 2D lattices with 3 unbalanced hopping terms, such that their combined effect is no longer either trivial or just that along a 1D subspace. We shall study two types of hopping configurations here, and reserve the third type, which turns out to interestingly requires a GBZ of lower dimensionality, to the next section on its own.

  \subsubsection{Type I: 2D lattice model with 3 hopping terms, 2 perpendicular to each other}
We first consider type I models, the simplest of ``entangled'' 2D lattice models. They contain 2 perpendicular terms $t_1,t_2$ in directions $(\alpha,0)$ and $(0,a)$, and a third one in an oblique direction $-(\beta,b)$. We have
  \begin{equation}\label{eqS210}
H_\text{2D,I}= \sum_{m,n}t_1|m,n\rangle\langle m+\alpha,n|+ t_2|m,n\rangle\langle m,n+a|+t_3|m,n\rangle\langle m-\beta,n-b|\ .
\end{equation}
Similarly as before, by substituting the ansatz $|\psi\rangle\propto\sum_{m,n}\psi_{mn}|m,n\rangle$ with $\psi_{m,n}=z_x^mz_y^n$ into the bulk equations,  we obtain the energy $E(z_x,z_y)$
  \begin{equation}\label{EqS211}
E(z_x,z_y)=t_1 z_x^{\alpha}+t_2z_y^a+t_3z_x^{-\beta}z_y^{-b}\ .
\end{equation}\\
By heuristically regarding $E(z_x,z_y)$ as 1D models $E(z_x)$ or $E(z_y)$ with nonconstant hoppings and using Eqs.~\ref{eqS3} and~\ref{eqS4}, or by considering boundary conditions on the wave function $|\Psi\rangle$ (Eq.\eqref{eqS21.5}), we obtain the following relationships between $z_x$ and $z_y$:
 \begin{equation}\label{EqS212}
 \begin{split}
z_x^{\alpha+\beta}z_y^b&=\frac{t_3}{t_1} z_1^{\alpha+\beta},\quad z_1^{\alpha+\beta}=\frac{\sin \beta k_1}{\sin \alpha k_1} \text{e}^{i(\alpha+\beta)k_1},\  k_1\in  \left(-\frac{\pi}{\alpha+\beta},\ \frac{\pi}{\alpha+\beta}\right]\ ,\\
z_x^{\beta}z_y^{a+b}&=\frac{t_3}{t_2} z_2^{a+b},\quad\  z_2^{a+b}=\frac{\sin b k_2}{\sin a k_2} \text{e}^{i(a+b)k_2},\quad k_2\in  \left(-\frac{\pi}{a+b},\ \frac{\pi}{a+b}\right]\ .
 \end{split}
\end{equation}
Substituting these $z_x$ and $z_y$ into the energy expression Eq.\eqref{EqS211}, we obtain the GBZ energies (Note the slight abuse of notation - below, we write $E(z_x,z_y)=E(z_1,z_2)$ to convey that the set of $E$ depends on $z_1,z_2$ through $z_x,z_y$, even though $E$ of course takes different functional dependencies in either case. )
 \begin{equation}\label{eqS213}
 \begin{split}
\bar E(k_1,k_2)=E(z_x,z_y)&=\left( \frac{t_1^{\beta a}t_2^{\alpha b}t_3^{\alpha a}}{z_1^{(\alpha+\beta)a\beta}z_2^{(a+b)b\alpha}} \right)^{\frac{1}{(\alpha+\beta)a+\alpha b}}\times\left(1+z_1^{\alpha+\beta}+z_2^{a+b}\right)\ , \\
 \end{split}
\end{equation}
and
 \begin{equation}\label{eqS214}
 \begin{split}
 \text{GBZ}_\text{2D,I}=\left\{ z_x,z_y \left|
z_x^{(\alpha+\beta)a+\alpha b}=\frac{t_2^{b}t_3^a}{t_1^{a+b}}\frac{z_1^{(\alpha+\beta)(a+b)}}{z_2^{(a+b)b}} \ ,\,
z_y^{(\alpha+\beta)a+\alpha b}=\frac{t_1^{\beta}t_3^{\alpha}}{t_2^{\alpha+\beta}}\frac{z_2^{(\alpha+\beta)(a+b)}}{z_1^{(\alpha+\beta)\beta}}\right.\right\}\ ,
 \end{split}
\end{equation}
with $z_1,z_2$ functions of $k_1,k_2$, as defined in Eq.\eqref{EqS212}. Some comments on the relationship between $z_x,z_y$ and $z_1,z_2$ are in order. As defined by the ansatz $|\psi(x,y)\rangle\propto \sum_{m,n}z_x^mz_y^n|m,n\rangle$, $z_x,z_y$ control the spatial growth and decay of the wavefunction. However, unlike in an ``unentangled'' case like $H_\text{2D,2}$, both $z_x$ and $z_y$ depend on the GBZ spanning parameters $k_1,k_2$ in complicated manners given by Eq.~\ref{EqS212}. In other words, Eq.~\ref{EqS212} dictates how the ``non-Bloch'' scaling factors $z_x,z_y$ are related to the GBZ coordinates $k_1,k_2$ through intermediate quantities $z_1,z_2$, which are related to the effective 1D chain projections of the Hamiltonian.

From this explicit expression Eq.~\ref{eqS214}, we can already deduce that the locus of $E$ is a star-like set of branches, since $E$ factorizes into a product of terms involving the phase factors of $k_1$ and $k_2$ separately. By considering these phase factors, we find that the number of branches for a type-I hopping lattice is given by
\begin{equation}
\text{LCM}\left[\displaystyle\frac{(\alpha+\beta)a+\alpha b}{\text{GCD}((\alpha+\beta)a+\alpha b,(\alpha+\beta)a\beta)},\displaystyle\frac{(\alpha+\beta)a+\alpha b}{\text{GCD}((\alpha+\beta)a+\alpha b,(a+b)\alpha b)}\right].
\end{equation}
Their agreement with numerical OBC results is given in FIG.~\ref{fig:S23}. The phase factor is given by $\exp(\frac{2\pi i( (\alpha+\beta)a\beta k_1+(a+b)\alpha bk_2}{(\alpha+\beta)a+\alpha b})$, which evaluates to $\exp(\frac{2\pi i(  2k_1+2k_2)}{3}) \Leftrightarrow 3$ branches for (a4); $\exp(\frac{2\pi i (  3k_1+4k_2)}{5}) \Leftrightarrow 5$ branches for (b4); and $\exp(\frac{2\pi( 3 k_1+12 k_2)}{7}) \Leftrightarrow 7$ branches for (c4).

\begin{figure}[h]
\centering
\includegraphics[width=.99\linewidth]{FIGS2-1.pdf}
\caption{Plots of the GBZ torus (Upper row a1-c1) , GBZs $z_x$ (Middle row a2-c2) and $z_y$ (Middle row a3-c3), as well as the energy spectra (Lower row a4-c4) of illustrative type-I lattices ($H_\text{2D,I}$ from Eq.~\eqref{eqS210}) with hoppings given by parameters $\alpha,\beta,a,b$. The upper row plots (a1-c1) are parametrized such that a torus is traced out in the trivial case without NHSE ($|z_x|=|z_y|=1$ for all $k_1,k_2$); departures from a toroidal shape depict the extent of 2D NHSE. Belonging to a 2D model, each of $z_x,z_y$ (a1-c1,a2-c2,a3-c3,) traces out a 2D region parametrized by $k_1,k_2$ (Eq.~\ref{eqS214}) as its GBZ, even though it trivially collapses into 1D loops for case (a). Perfect agreement of GBZ spectra $E(z_x,z_y)$ from Eq.~\eqref{eqS213} with OBC spectra $E_{\text{OBC}}$ is demonstrated for all cases, which for this model fills the interior of a $[(\alpha+\beta)a+\alpha b]$-sided figure (a4-c4).
Parameters are $L_x=L_y=15$, $t_1=t_1=t_3=1$, and the GBZ predictions are generated with a mesh defined by $k_1=-\frac{\pi}{\alpha+\beta}:\frac{\pi}{30}:\frac{\pi}{\alpha+\beta}$,  $k_2=-\frac{\pi}{a+b}:\frac{\pi}{30}:\frac{\pi}{a+b}$.}\label{fig:S23}
\end{figure}

\clearpage
\subsubsection{Type II: 2D lattice model with 3 hopping terms, two with common $\beta$ displacement}
\begin{figure}[h]
\centering
\includegraphics[width=.99\linewidth]{FIGS2-2.pdf}
\caption{
Plots of the  GBZ torus (Upper row a1-c1) , GBZs $z_x$ (Middle row a2-c2) and $z_y$ (Middle row a3-c3), as well as the energy spectra (Lower row a3-c3)of illustrative type-II lattices ($H_\text{2D,II}$ from Eq.~\eqref{eqS215}) with hoppings given by parameters $\alpha,\beta,a,b$. (a1-c1,a2-c2,a3-c3) While $z_x$ traces out a 2D region parametrized by the two ``momenta'' $k_1$ and $k_2$, $z_y$ only depends on one such momentum parameter, as given by Eq.~\eqref{eqS216}), thereby tracing out a 1D GBZ loop (a3-c3). (a4-c4) Perfect agreement of GBZ spectra $E(z_x,z_y)$ from Eq.~\eqref{eqS216e} with OBC spectra $E_{\text{OBC}}$ is demonstrated for all cases, with their star-like spectra consistent with the 1D nature of their effective GBZ description. The $|z_x|=0$ central dot corresponds to $k_2=\pm\pi/(a+b)$ where the factor $\left(z_2^a+z_2^{-b}\right)$ in $z_x$ is equal to 0. Parameters are $L_x=L_y=15$, $t_1=t_1=t_3=1$, and the GBZ predictions are generated with a mesh defined by $k_1=-\frac{\pi}{\alpha+\beta}:\frac{\pi}{30}:\frac{\pi}{\alpha+\beta}$, $k_2=-\frac{\pi}{a+b}:\frac{\pi}{30}:\frac{\pi}{a+b}$.
}
\label{fig:S24}
\end{figure}

We next consider a slightly more complicated 2D lattice model (type II). None of the hoppings are orthogonal to each other, so all 3 hoppings are ``entangled''. However, it has the simplifying property that the $t_2$ and $t_3$ hoppings are equidistant in the direction parallel to the $t_1$ hopping, such that if we set the $t_1$ hopping  normal against the horizontal $x$ open boundary, there are only two unique hopping distances in this direction. The Hamiltonian is given by
  \begin{equation}
H_\text{2D,II}= \sum_{m,n}t_1|m,n\rangle\langle m+\alpha,n|+ t_2|m,n\rangle\langle m-\beta,n+a|+t_3|m,n\rangle\langle m-\beta,n-b|\ .
\label{eqS215}
\end{equation}
Taking the ansatz $|\psi\rangle\propto\sum_{m,n}\psi_{mn}|m,n\rangle$ with $\psi_{m,n}=z_x^mz_y^n$ into the bulk equations, we obtain the energy $E(z_x,z_y)$
  \begin{equation}\label{eqS216e}
E(z_x,z_y)=t_1 z_x^{\alpha}+t_2z_x^{-\beta}z_y^a+t_3z_x^{-\beta}z_y^{-b}\ .
\end{equation}
Considering the boundary conditions in the same way as before, we have two relations between $z_x$ and $z_y$ from which the GBZ and OBC energy $E(z_x,z_y)$ can be obtained:
 \begin{equation}\label{eqS216}
 \begin{split}
 \text{GBZ}_\text{2D,II}=\left\{z_x,z_y\left|z_x^{\alpha+\beta}=t_1^{-1}t_2^{\frac{b}{a+b}}t_3^{\frac{a}{a+b}}\left(z_2^a+z_2^{-b}\right)z_1^{\alpha+\beta}, z_y^{a+b}=\frac{t_3}{t_2} z_2^{a+b}\right.\right\}\ ,
 \end{split}
\end{equation}
\begin{equation}\label{eqS217}
\bar E(k_1,k_2)=E(z_x,z_y)=t_1^{\frac{\beta}{\alpha+\beta}} \left(t_2^{\frac{b}{a+b}}t_3^{\frac{a}{a+b}}\left(z_2^a+z_2^{-b}\right)\right)^{\frac{\alpha}{\alpha+\beta}}\left(z_1^{\alpha}+z_1^{-\beta}\right),
\end{equation}
the forms of $z_1,z_2$ given in Eq.\eqref{EqS212} as before.

From the explicit expression Eq.~\ref{eqS217}, we can similarly deduce that the locus of $E$ is a star-like set of branches. We find that the number of branches for type-II lattice hoppings is given by
\begin{equation}
\text{LCM}\left[\displaystyle\frac{(a+b)(\alpha+\beta)}{\text{GCD}((a+b)(\alpha+\beta),a\alpha)},\frac{\alpha+\beta}{\text{GCD}(\alpha+\beta,\alpha)}\right].
\end{equation}
We verify this formula with the three examples in FIG.~\ref{fig:S24}, for which there exist excellent agreement between the numerical OBC results and the above GBZ expression for $E$. The phase factor is given by$\exp(\frac{2\pi( a\alpha k_1+\alpha(a+b))k_2}{(a+b)(\alpha+\beta)}$, which evaluates to $\exp(\frac{2\pi i(  k_1+2k_2)}{4}) \Leftrightarrow 4$ branches for (a4); $\exp(\frac{2\pi i (  k_1+k_2)}{6/2}) \Leftrightarrow 3$ branches for (b4); and $\exp(\frac{2\pi( 2 k_1+6k_2)}{9}) \Leftrightarrow 9$ branches for (c4). It is not surprising that these type II 2D lattice models considered above have star-like spectra that resemble that of 1D NHSE models, since they after quickly reduce to a simple 1D effective model with x-OBCs with 2 effective hoppings.

\begin{comment}
\subsection{Failure of successively taking GBZs in different directions}
\textcolor{magenta}{CH: Will be good to have another figure showing additional examples where our approach gives excellent agreement with the numerical spectrum, but which taking x and y OBCs sequentially doesnt. Maybe plot something analogous to the spectral plots in the above figures, but with points of another color representing the sequential approach that doesnt agree. No need to describe these models in the text (I will intro them somewhere) We can put this subsection either here, or in Sect III.}
\begin{figure}[h]
\centering
\includegraphics[width=5.0in]{gFIGS25.pdf}
\caption{\textcolor[rgb]{0,0.4,0.8}{JH: I think the type I model with $\alpha=\beta=a=b=1$ is the best example to show that taking x and y OBCs sequentially doesn't give excellent agreement with the numerical spectrum. I also take the type I model with $\alpha=a=b=1,\beta=2$ in (b) as another example, it not as good as (a) because the size of system. And for type II model , taking x and y OBCs sequentially have the same results with our  our approach. If we consider the model without type I II,  it seems  that the connection with the context is not so strong. And there have same results with different x and y OBC sequence. }
}
\label{fig:S25}
\end{figure}
\end{comment}

\newpage
\section{2D lattices with dimensionally-reduced GBZs}\label{sec3}

The previous section gives the approach to obtain the GBZ of 2D lattice models. However, as explained in the main text, there are some classes of models where the correct GBZ is not even of the same dimensionality as the lattice. In this section, we shall provide a detailed account of the additional steps and analysis required to dimensionally reduce the GBZ to the correct one.

\subsection{Model description and setup}

We consider 2D lattice models of the form (type III in Fig.~\ref{fig:S221})
  \begin{equation}\label{eqS31}
H_\text{2D,III}=\sum_{m,n}t_1|m,n\rangle\langle m+\alpha,n+a|+t_2|m,n\rangle\langle m-\beta,n+a|+t_3|m,n\rangle\langle m-\beta,n-b|\ ,
\end{equation}
with oblique hoppings of amplitudes $t_1,t_2,t_3$. They constitute minimal models that require dimensional reduction of the GBZ, and are still simple enough such that they can be analyzed completely analytically.  Using the ansatz $|\psi\rangle\propto\sum_{m,n}\psi_{mn}|m,n\rangle$ with $\psi_{m,n}=z_x^mz_y^n$ as before, the bulk relations give $t_1\psi_{m+\alpha,n+a}+t_2\psi_{m-\beta,n+a}+t_3\psi_{m-\beta,n-b}=E(z_x,z_y)\psi_{m,n}$, which yields
\begin{equation}\label{eqS31.5}
\begin{split}
E(z_x,z_y)=\,t_1z_x^{\alpha}z_y^{a}+t_2z_x^{-\beta}z_y^{a}+t_3z_x^{-\beta}z_y^{-b}\
\end{split}
\end{equation}

\begin{comment}
Fixed $z_x$, we have $a+b$ solutions $z_{y,i}$, ($i=1,2,...,a+b$)  corresponding to the same energy $E(z_x,z_y)$. Similarly, there also find $\alpha+\beta$ solutions $z_{x,j}$, ($j=1,2,...,\alpha+\beta$) by fixed $E(z_x,z_y)$ and $z_y$. Both solutions $z_{y,i}$, ($i=1,2,...,a+b$) and $z_{x,j}$, ($j=1,2,...,\alpha+\beta$)  are the functions of $z_x$ and $z_y$.  Similarly, there have lots of sets $z_x,z_y$ corresponding to same energy $E$.\\

Considering the boundary conditions, the wave function $|\Psi\rangle$ with same eigenvalue $E(z_x,z_y)$ can be despicte
by
    \begin{equation}\label{eqS29}
    \begin{split}
&\quad \quad |\Psi\rangle=|\Psi^1\rangle=|\Psi^2\rangle\ \propto \sum_{m,n}\Psi_{mn}|m,n\rangle=\sum_{m,n}\Psi^1_{mn}|m,n\rangle=\sum_{m,n}\Psi^2_{mn}|m,n\rangle\ ,\\
& \Psi^1_{mn}=\sum_{z_y} \left( \left(\sum_{j}^{\alpha+\beta}f_{z_{x,j}}z_{x,j}^m\right)g_{z_y} z_y^n\right)\ , \quad \Psi^2_{mn}=\sum_{z_x}\left( f_{z_x} z_x^m \left(\sum_{i}^{a+b}{g}_{z_{y,i}}z_{y,i}^n\right)\right)\ ,\\
\end{split}
\end{equation}
$|\Psi^1\rangle,\ |\Psi^2\rangle$ are the different manifestations of the wave function $|\Psi\rangle$,  which are equivalent with $f_{z_x}g_{z_{y,i}}=g_{z_{y,i}}f_{z_x}$, and the $\sum$ in equation $\Psi^{1,2}_{mn}$ mean we must consider all the sets $z_x,z_y$ which have same energy $E$.\\

Taking the  different manifestations of the wave function $|\Psi\rangle$ into boundary conditions along different direction (i.e.,${\bm x, \bm y}$ direction), like boundary conditions along ${\bm x}$ direction, we have  $\Psi_{m'{n}}=\Psi^1_{m'{n}}=0 $ with special value $m'$( like $m'=-\beta+1,-\beta+2,...0,L+1,L+2,..L+\alpha $ if $\alpha,\beta>0$ )and any values $n$, like `1D' system with paramters $z_y$ in Sec.~\ref{sec1}.
Thus,  for the same energy, we have
\end{comment}

By treating the system as a quasi-1D system in $x$ or $y$ as before i.e. by expressing it in the form of Eq.~\ref{eqS1} with GBZ and energies given by Eqs.~\ref{eqS3} and~\ref{eqS4}, we obtain
\begin{equation}\label{eqS32}
\begin{split}
z_x^{\alpha+\beta}=\frac{t_2+t_3z_y^{-a-b}} {t_1}z_1^{\alpha+\beta}, \quad z_1^{\alpha+\beta}=\frac{\sin \beta k_1}{\sin \alpha k_1}\text{e}^{i(\alpha+\beta)k_1},\ k_1\in  \left(-\frac{\pi}{\alpha+\beta},\ \frac{\pi}{\alpha+\beta}\right]\ ,
\end{split}
\end{equation}
and
\begin{equation}\label{eqS33}
\begin{split}
z_y^{a+b}=\frac{t_3}{t_1z_x^{\alpha+\beta}+t_2} z_2^{a+b}, \quad z_2^{a+b}=\frac{\sin bk_2}{\sin ak_2}\text{e}^{i(a+b)k_2},\ \  k_2\in  \left(-\frac{\pi}{a+b},\ \frac{\pi}{a+b}\right].
\end{split}
\end{equation}
That is, the boundary conditions along the two different direction give us two relations on $z_x$ and $z_y$ respectively. In order to get $z_x,z_y$ which satisfy all the boundary conditions, we solve the above to obtain
\begin{equation}\label{eqS34}
\begin{split}
z_x^{\alpha+\beta}=\frac{t_2}{t_1}\frac{1+z_2^{a+b}}{z_2^{a+b}-z_1^{\alpha+\beta}}z_1^{\alpha+\beta},\quad z_y^{a+b}=\frac{t_3}{t_2}\frac{z_2^{a+b}-z_1^{\alpha+\beta}}{1+z_1^{\alpha+\beta}}\ .
\end{split}
\end{equation}
Note that now, $z_x$ and $z_y$ are not even proportional to a single phase factor, and thus no longer take any conventional ``non-Bloch'' form. Substituting $z_x,z_y $ Eq.\eqref{eqS34} into the energy $E(z_x,z_y)$ Eq.\eqref{eqS31.5}, we furthermore obtain
\begin{equation}\label{eqS35}
\begin{split}
\bar E(k_1,k_2)=E(z_x,z_y)&=\frac{\left(t_1^{\beta}(1+z_2^{a+b})^{\alpha}\right)^{\frac{1}{\alpha+\beta}}\cdot\left(t_3^{a}(1+z_1^{\alpha+\beta})^{b}\right)^{\frac{1}{a+b}}}{((z_2^{a+b}-z_1^{\alpha+\beta})/t_2)^{\frac{\alpha b-\beta a}{(\alpha+\beta)(a+b)}}}z_1^{-\beta} \
\end{split}
\end{equation}
with $v=1,2,...,a+b$, the forms of $z_1,z_2$ given by Eq.(\ref{eqS32} and~\ref{eqS33}).

Importantly, this 2D GBZ $(z_x,z_y)$ as it is currently defined does \emph{not} form a valid GBZ because there exists certain paths on it where the spectral winding
 \begin{equation}\label{eqS36}
   \begin{split}
   \omega_{i}(E)&=\oint\text{d} z_i \partial_{z_i} \log(E(\bm k)-E_0)\ ,\\
   \end{split}
\end{equation}
$i=1,2 $ is nonzero for some arbitrary reference energies $E_b$. In other words, there exist some closed paths in $(k_1,k_2)$ space (defined in Eqs.~\ref{eqS32},~\ref{eqS33}) such that the energy $E(k_1,k_2)$ loop encloses a nonzero area as we cycle over $k_1$ or $k_2$. An illustrative example is shown in Fig.~\ref{fig:S31}.

\begin{figure}[htb]
\centering
\includegraphics[width=5.5in]{FIGS3-1.pdf}
\caption{
(a,b) Closed paths from the unconstrained GBZ from Eq.~\ref{eqS32} to~\ref{eqS34}  that result in nontrivial spectral winding, in contradiction to the known fact that OBC spectra should not enclose any nonzero area (nontrivial winding). The unconstrained GBZ $z_x$ (a) and $z_y$ (b) of $H_\text{2D,III}$ with $t_1=t_2=t_3=1$ are plotted in black for $k_1,k_2$ ranging over all values in steps of $\pi/100$, with illustrative green and blue closed paths enclosing nonzero area given by $k_1=(-\pi/3,\pi/3), k_2=\pi/6$ and $k_2=(-\pi/3,\pi/3), k_1=\pi/6$ respectively. (c) The corresponding spectrum (Eq.~\ref{eqS35}) from the unconstrained GBZ (black) and the corresponding green and blue spectral loops that enclose non-vanishing areas, which they are not supposed to do. In particular, they do not agree with the correct, numerically obtained OBC spectrum (red), which lies on a subset of $E(z_x,z_y)$ points, but which does not exhibit any nontrivial winding.
 }\label{fig:S31}
\end{figure}

Since it is known\cite{okuma2020topological,zhang2019correspondence} that OBC spectra should never possess nonzero spectral winding, which presupposes incomplete NHSE equilibration, the above results Eqs.~\ref{eqS34} and~\ref{eqS35} cannot possibly give the correct OBC spectrum of $H_\text{2D, III}$. Indeed, this is shown in the lack of agreement between the numerically obtained $E_\text{OBC}$ spectrum and the (currently unconstrained) GBZ spectrum $E(z_x,z_y)$ in Fig.~\ref{fig:S31}. Yet, the above results can be rigorously traced to satisfy all the relations pertaining to the real-space hoppings, and cannot be incorrect. Hence the natural conclusion is that the correct GBZ must be a subset of the GBZ, as currently defined by Eqs.~\ref{eqS32} to~\ref{eqS34} with the range of $k_1,k_2$ prescribed above. In the below, we shall carefully derive the constraints that extracts the subset of $k_1,k_2$ that generates the correct GBZ for this model; at the end of this supplement, we shall generalize this correct GBZ construction to arbitrary non-Hermitian lattice models.

\subsection{Detailed analysis of constraints leading to dimensional reduction}\label{sec32}

Here, we elaborate on the possible constraints that can make the energy spectrum (Eq.~\ref{eqS35}) exhibit zero winding. Since the unconstrained GBZ is a 2D torus that is spanned by 2 homotopy generators, a sure way to remove all possible nonzero spectral windings is to remove the homotopy generator/s that lead to particular nontrivial spectral winding paths. In principle, there can be many ways to remove a homotopy generator, since a combination of generators form another generator. However, for our model $H_\text{2D,III}$, it turns out that the requisite constraints can simply take the linear form of $k_1=\gamma k_2$, with $\gamma$ taking particular values that we shall elaborate on. (Note that we cannot possibly remove both homotopy generators, since there will be no GBZ left then.)\\

\subsubsection{Analysis of different possible parametrizations}
We now justify why it suffices to consider $k_1=\gamma k_2$, such that the energy takes the 1-parameter form $\bar E(k_2)$ (or equivalently $\bar E(k_1)$), with an slight abuse of notation. First, we establish the range and offset of $k_1,k_2$. given that $z_1,z_2$ (Eqs.~\ref{eqS32},\ref{eqS33}) satisfy periodicity conditions $z_{1}(k_{1})=z_{1}(k_{1}+2m'\pi/(\alpha+\beta))$, $z_{2}(k_{2})=z_{2}(k_{2}+2n'\pi/(a+b))$ with integer $n',m'$ ($n',m' \in \mathbb{Z}$). We can hence write down a valid linear reparametrization as $k_1=\gamma k_2+{2 n' \pi \gamma}/{(a+b)}+{2 m' \pi}/{(\alpha+\beta)}$ with integer $n',m'$, such that the energy Eq.\eqref{eqS35} is transformed into $\bar E'$ as given by
\begin{equation}
\begin{split}
   &\bar E'= E(z_x,z_y) \Big |_{k_1=\gamma k_2+\frac{2 n' \pi \gamma}{a+b}+\frac{2 m' \pi}{\alpha+\beta}}= E\left(z_1(k_1),z_2(k_2)\right) \Big |_{k_1=\gamma k_2+\frac{2 n' \pi \gamma}{a+b}+\frac{2 m' \pi}{\alpha+\beta}}\\
&\quad=E(z_1(\gamma k_2+2n'\pi\gamma/(a+b)),z_2(k_2))\\
&\quad= \bar E\left(k_2+\frac{2n'\pi}{a+b}\right)\ ,
\end{split}
\end{equation}
such as $E'\in \{E(k_2)\}$. Since this hence just a translation of a single coordinate, we can hence just consider constraints of the form $k_1=\gamma k_2$ without any constant offset. With that,
\begin{equation}
\bar E(k_2)= E(z_x,z_y) \Big |_{k_1=\gamma k_2}= E\left(z_1(k_1),z_2(k_2)\right) \Big |_{k_1=\gamma k_2}=E(z_1(\gamma k_2),z_2(k_2))\ .
\end{equation}\\

Next, we shall show that there exists $3$ possible values of $\gamma$ that guarantee zero spectral winding number while keeping $z_x,z_y$ periodic, namely $\gamma=b/\beta,\ a/\alpha,\  (a+b)/(\alpha+\beta)$.
Below are the detailed justifications for the possible choices of $\gamma$ (we call them $\gamma_1,\gamma_2$ and $\gamma_3$)
\begin{itemize}
\item The case with $k_1=\gamma_1 k_2$, $\gamma_1=b/\beta$: Here $z_x,z_y$~\eqref{eqS34} and energy~\eqref{eqS35} with $k_1=b k, k_2=\beta k$ are given by
\begin{equation}\label{eqS37}
\begin{split}
 &z_{x,1}^{\alpha+\beta} =\left. z_x^{\alpha+\beta}\right|_{k_1=\gamma_1 k_2}=\frac{t_2\sin((a+b) \beta k)}{t_1\sin((\alpha b-a\beta )k)}\text{e}^{i(\alpha b+\beta b)k}\ , \\
& \left.  z_{y,1}^{a+b} = z_y^{a+b}\right|_{k_1=\gamma_1 k_2}=\frac{t_3\sin(\beta b k)\sin((\alpha b- a \beta)k)}{t_2\sin( a \beta k)\sin((\alpha b+\beta b)k)}\ ,
 \end{split}
\end{equation}
\begin{equation}\label{eqS38}\begin{split}
\bar E_1(k)&=\left. E(z_x,z_y)\right|_{k_1=\gamma_1 k_2=bk}\ ,\\
 &=t_1^{\frac{\beta}{\alpha+\beta}}t_2^{\frac{\alpha b-\beta a}{(\alpha+\beta)(a+b)}}t_3^{\frac{a}{a+b}} \text{e}^{\frac{i2\pi av}{a+b}+\frac{2\pi\alpha v'}{\alpha+\beta}} \left(\frac{\sin(\alpha b+\beta b) k}{\sin \alpha bk}\right)^{\frac{b}{a+b}}\\
&\qquad \times\left(\frac{\sin(a\beta+\beta b)k}{\sin a\beta k}\right)^{\frac{\alpha}{\alpha+\beta}}\left(\frac{\sin\beta  bk}{\sin \alpha bk}\right)^{-\frac{\beta}{\alpha+\beta}} \left(\frac{\sin(a+b)\beta k}{\sin a\beta k}-\frac{\sin(\alpha b+\beta b)k}{\sin \alpha b k}\right)^{\frac{\beta a-\alpha b}{(\alpha+\beta)(a+b)}}\ ,\\
\end{split}
\end{equation}
where $v=1,2,..,a+b$ ,  $v'=1,2,..,\alpha+\beta$. The range of $k$ has to be dependent on the ratio between $\alpha/\beta$ and $a/b$ as follows:

\begin{itemize}
    \item[-]  $\alpha/\beta>a/b$: $k$ can take values in $\left(-\frac{\pi}{(\alpha+\beta) b},\frac{\pi}{(\alpha+\beta) b}\right]$.  Both $z_{x,1}$ and $z_{y,1}$ are connected end to end with `momentum' $k$, and the path of $z_{x,1}$ in the complex plane forms a closed loop.  Hence, the results of $z_{x,1},z_{y,1}$ (Eq.\eqref{eqS37}) are legitimate, and zero winding in the energy $\bar E_1$ (Eq.\eqref{eqS38}) is respected.
    \item[-] $\alpha/\beta=a/b$:  $k$ can take values in $\left(-\frac{\pi}{(a+b) \beta},\frac{\pi}{(a+b) \beta}\right]$. And $z_{x,1}^{\alpha+\beta}\rightarrow \infty$, \ $z_{y,1}^{a+b}\rightarrow 0$, but the path of $ z_{x,1}^{\alpha+\beta}z_{y,1}^{a+b}=\frac{\sin(\beta bk)}{\sin(a\beta k)}\text{e}^{i(a+ b)\beta k}$ in the complex plane  forms a closed loop with `momentum'  $k$.  Hence, the results of $z_{x,1},z_{y,1}$ Eq.\eqref{eqS37} are legitimate,  and zero winding in the  energy $\bar E_1$ Eq.\eqref{eqS38} is respected.
   \item[-] $\alpha/\beta<a/b$:  $k$ can  take values in  $\left(-\frac{\pi}{(a+b) \beta},\frac{\pi}{(a+b) \beta}\right]$. None of   $z_{x,1}^{\alpha+\beta}$, $z_{y,1}^{a+b}$, $z_{x,1}^{\alpha+\beta}z_{y,1}^{a+b}$ can form  closed loop parametrized by  `momentum' $k$. The results of $z_{x,1},z_{y,1}$ and the zero winding requirement of the energy $\bar E_1$ Eq.\eqref{eqS38} are inconsistent. It is worth noting that when $z^{\alpha+\beta}_1=z_2^{a+b}$ in Eq.\eqref{eqS34}, the $z_{x}^{\alpha+\beta}\rightarrow \infty,z_{y}^{a+b}\rightarrow 0$ which can ignore the end-to-end condition. And in this case, the energy $\bar E=0$.

\end{itemize}
\item The case with $k_1=\gamma_2 k_2$, $\gamma_2=a/\alpha$: $z_x,z_y$~\eqref{eqS34} and the energy~\eqref{eqS35} with $k_1=a k', k_2=\alpha k'$ are given by
\begin{equation}\label{eqS39}
\begin{split}
& z_{x,2}^{\alpha+\beta} =\left. z_x^{\alpha+\beta}\right|_{k_1=\gamma_2 k_2}=\frac{t_2\sin(\beta a k')\sin((\alpha a+ \alpha b)k')}{t_1\sin( a \alpha k')\sin((\alpha b-\beta a)k')}\ , \\
 &\left.  z_{y,2}^{a+b} = z_y^{a+b}\right|_{k_1=\gamma_2 k_2}=\frac{t_3\sin((\alpha b-a\beta )k'))}{t_2\sin((\alpha a+a\beta )k')}\text{e}^{i(\alpha a+\alpha b)k'}\ ,
 \end{split}
\end{equation}
\begin{equation}\label{eqS310}\begin{split}
\bar E_2(k')&=\left. E(z_x,z_y)\right|_{k_1=\gamma_2 k_2=ak'}\ ,\\
 &=t_1^{\frac{\beta}{\alpha+\beta}}t_2^{\frac{\alpha b-\beta a}{(\alpha+\beta)(a+b)}}t_3^{\frac{a}{a+b}} \text{e}^{\frac{i2\pi av}{a+b}+\frac{2\pi\alpha v'}{\alpha+\beta}} \left(\frac{\sin(\alpha a+\beta a) k'}{\sin \alpha ak'}\right)^{\frac{b}{a+b}}\\
&\qquad \times\left(\frac{\sin(a\alpha+\alpha b)k'}{\sin a\alpha k'}\right)^{\frac{\alpha}{\alpha+\beta}}\left(\frac{\sin\beta ak'}{\sin \alpha ak'}\right)^{-\frac{\beta}{\alpha+\beta}} \left(\frac{\sin(\alpha b-a\beta ) k'}{\sin a\alpha k'}\right)^{\frac{\beta a-\alpha b}{(\alpha+\beta)(a+b)}}\ ,\\
\end{split}
\end{equation}
where $v=1,2,..,a+b$ , $v'=1,2,..,\alpha+\beta$. The range of $k'$ similarly has to depend on the ratio between $\alpha/\beta$ and $a/b$:
\begin{itemize}
    \item[-]  $\alpha/\beta>a/b$,  $k'$ can take values in $\left(-\frac{\pi}{(a+b) \alpha},\frac{\pi}{(a+b) \alpha}\right]$.  Both $z_{x,2}$ and $z_{y,2}$ are connected end to end with `momentum' $k$, and the path of $z_{y,2}$ in complex plane forms a closed loop. Hence, the results of $z_{x,2},z_{y,2}$ (Eq.\eqref{eqS39}) are legitimate, and  zero winding in the energy $\bar E_2$ Eq.\eqref{eqS310} is repected.
    \item[-] $\alpha/\beta=a/b$,  $k'$ can take values in  $\left(-\frac{\pi}{(a+b) \alpha},\frac{\pi}{(a+b) \alpha}\right]$. And  $z_{x,2}^{\alpha+\beta}\rightarrow \infty$, \ $z_{y,2}^{a+b}\rightarrow 0$, but the path of $ z_{x,2}^{\alpha+\beta}z_{y,2}^{a+b}=\frac{\sin(\alpha b k')}{\sin(\alpha a k')}\text{e}^{i( a+ b)\alpha k'}$  in the complex plane  forms a closed loop with `momentum' $k$. Hence, the results of $z_{x,2},z_{y,2}$ Eq.\eqref{eqS39} are legitimate, and  zero winding in the energy $\bar E_2$ Eq.\eqref{eqS310} is repected.
   \item[-] $\alpha/\beta<a/b$,  $k'$ can take values in  $\left(-\frac{\pi}{(\alpha+\beta) a},\frac{\pi}{(\alpha+\beta) a}\right]$. None of   $z_{x,2}^{\alpha+\beta}$, $z_{y,2}^{a+b}$, $z_{x,2}^{\alpha+\beta}z_{y,2}^{a+b}$ can form closed loop parametrized by `momentum' $k'$. The results of $z_{x,2},z_{y,2}$ and the zero winding requirement of the energyy $\bar E_2$ Eq.\eqref{eqS310} seem unreasonable. It is worth noting that when $z^{\alpha+\beta}_1=z_2^{a+b}$ in Eq.\eqref{eqS34}, the $z_{x}^{\alpha+\beta}\rightarrow \infty,z_{y}^{a+b}\rightarrow 0$ which can ignore the end-to-end condition. And in this case, the energy $\bar E=0$.

\end{itemize}

\item The case with $k_1=\gamma_3 k_2$, $\gamma_3=(a+b)/(\alpha+\beta)$: $z_x,z_y$~\eqref{eqS34} and energy~\eqref{eqS35}  with  $k_1=(a+b) k'', k_2=(\alpha+\beta) k''$  take the forms
\begin{equation}\label{eqS311}
\begin{split}
& z_{x,3}^{\alpha+\beta} =\left. z_x^{\alpha+\beta}\right|_{k_1=\gamma_3 k_2}=\frac{t_2\sin (\beta (a+b)k'')}{t_1\sin(\alpha b-\beta a)k'')}\text{e}^{ib(\alpha+\beta)k''}\ , \\
 &\left.  z_{y,3}^{a+b} = z_y^{a+b}\right|_{k_1=\gamma_3 k_2}=\frac{t_3\sin(\alpha b-\beta a)k'')}{t_2\sin(a(\alpha+\beta)k'')}\text{e}^{i\alpha (a+b)k''}\ ,
 \end{split}
\end{equation}
\begin{equation}\label{eqS312}\begin{split}
\bar E_3(k'')&=\left. E(z_x,z_y)\right|_{k_1=\gamma_3 k_2=(a+b)k''}\ ,\\
 &=t_1^{\frac{\beta}{\alpha+\beta}}t_2^{\frac{\alpha b-\beta a}{(\alpha+\beta)(a+b)}}t_3^{\frac{a}{a+b}} \text{e}^{\frac{i3\pi av}{a+b}+\frac{3\pi\alpha v'}{\alpha+\beta}}  \left(\frac{\sin(\alpha+\beta)(a+b)k''}{\sin \alpha (a+b)k''}\right)^{\frac{b}{a+b}}\left(\frac{\sin(a+b)(\alpha+\beta)k''}{\sin a(\alpha+\beta)k''}\right)^{\frac{\alpha}{\alpha+\beta}}\\
&\qquad \times\left(\frac{\sin \beta (a+b)k''}{\sin \alpha (a+b)k''}\right)^{-\frac{\beta}{\alpha+\beta}}\left(\frac{\sin b (\alpha+\beta)k''}{\sin a(\alpha+\beta)k''}-\frac{\sin \beta (a+b)k''}{\sin \alpha (a+b)k''}\right)^{\frac{\beta a-\alpha b}{(\alpha+\beta)(a+b)}}\ ,\\
\end{split}
\end{equation}
where $v=1,2,..,a+b$ , $v'=1,2,..,\alpha+\beta$. Here, the range of $k''$ is always $\left(-\frac{\pi}{(\alpha+\beta)(a+b)},\frac{\pi}{(\alpha+\beta)(a+b)}\right]$.
 \begin{itemize}
    \item[-]  $\alpha/\beta\neq a/b$: none of $z_{x,3}^{\alpha+\beta}$, $z_{y,3}^{\alpha+\beta}$, $z_{x,3}^{\alpha+\beta}z_{y,3}^{\alpha+\beta}$ can form closed loop parametrized by momentum $k''$. Hence the results of $z_{x,3},z_{y,3}$ (Eq.\eqref{eqS311}) are inadmissible, since they are inconsistent with the zero winding requirement of the energy $E_3$ Eq.\eqref{eqS312}.
        \item[-] $\alpha/\beta=a/b$: this case is equivalent to the case $k_1=\gamma_{1} k_2$ or $k_1=\gamma_{2} k_2$ with $\alpha/\beta=a/b$.
\end{itemize}
Hence, the case with $k_1=\gamma_3 k_2$, $\gamma_3=(a+b)/(\alpha+\beta)$ can  incorporate the cases with $k_1=\gamma_{1,2} k_2$ ($\gamma_1=b/\beta,\gamma_2=a/\alpha$).
\end{itemize}

\subsubsection{Summary}
\begin{table*}[h!]
\renewcommand\arraystretch{1.5}
  \centering
  \caption{Brief summary of admissible GBZs for our type III model (Eq.~\ref{eqS31})}
  {(Shaded boxes represent valid contributing GBZ sectors)}\label{tab:sum}
  \begin{tabular}{c|c|c|c}
\hline
\hline
& \bm{ $k_1=\gamma_1k_2$} {\small$(\gamma_1=\frac{b}{\beta})$} & \bm{$k_1=\gamma_2k_2$} {\small $(\gamma_2=\frac{a}{\alpha})$}& \bm{ $k_1=\gamma_3k_2$}{\small$(\gamma_3=\frac{a+b}{\alpha+\beta})$}\\
&$z_{x(y),1}=\left. z_{x(y)}\right|_{k_1=\gamma_1 k_2}$ &$z_{x(y),2}=\left. z_{x(y)}\right|_{k_1=\gamma_2 k_2}$&$z_{x(y),3}=\left. z_{x(y)}\right|_{k_1=\gamma_3 k_2}$\\
\hline
\multirow {3}{*}{${\alpha/\beta>a/b}$}&\multicolumn{1}{>{\columncolor{mygray}}c|}{\bm{$z_{x,1}$}{\bf: closed loop}}  & \multicolumn{1}{>{\columncolor{mygray}}c|}{{$z_{x,2}$}{: real}}  &  \multirow {2}{*}{$z_{x,3}$, $z_{y,3}$, $z_{x,3}^{\alpha+\beta}z_{y,3}^{\alpha+\beta}$:  open loop}\\
&\multicolumn{1}{>{\columncolor{mygray}}c|}{{$z_{y,1}$}{: real}} &\multicolumn{1}{>{\columncolor{mygray}}c|}{\bm{$z_{y,2}$}{\bf: closed loop}} &\\
\cline{2-4}
& \multicolumn{3}{c}{ Hence {GBZ=$\text{GBZ}_1\cup\text{GBZ}_2$}\quad(with $\text{GBZ}_1=\left\{z_{x,1},z_{y,1}\right\}$, $\text{GBZ}_2=\left\{z_{x,2},z_{y,2}\right\}$)}\\
\hline
\multirow {3}{*}{${\alpha/\beta=a/b}$}&\multicolumn{3}{>{\columncolor{mygray}}c}{$z_{i,1}=z_{i,2}=z_{i,3}$\quad ($i=x,y$) and $z_{x,1}\rightarrow \infty,\ z_{y,1}\rightarrow 0$}  \\
&\multicolumn{3}{>{\columncolor{mygray}}c}{\bm{$z^{\alpha+\beta}_{x,1} z^{a+b}_{y,1}$}{\bf: closed loop}} \\
\cline{2-4}
& \multicolumn{3}{c}{Hence $\text{GBZ}=\left\{z^{\alpha+\beta}_{x,1}z^{a+b}_{y,1}\right\}$ }\\
\hline
\multirow {2}{*}{${\alpha/\beta<a/b}$}&\multicolumn{3}{c}{$z_{x,j}\rightarrow \infty, \ z_{y,j}\rightarrow 0 $ \quad ($j=1,2,3$)}  \\
&\multicolumn{3}{c}{ $\bar E=0$}  \\
\hline
\hline
\end{tabular}
\end{table*}

To summarize, the (correct) constrained GBZ of $H_\text{2D, III}$ of Eq.\eqref{eqS31} is not 2D, but is made up of one or more disconnected 1D loops as described above, depending on the hopping lengths $\alpha,\beta, a,b$. We call the effective Hamiltonian constrained to this reduced GBZ $H_\text{2D-red}$:
 \begin{itemize}
    \item[$\blacktriangleright$]   When  $\alpha/\beta>a/b$, the GBZ and spectrum consists of the union of the two sectors
    \begin{equation}\label{eqS313}
\begin{split}
&\ \quad\qquad\qquad\text{GBZ}=\text{GBZ}_1\cup\text{GBZ}_2,\\
&\text{GBZ}_1=\left\{z_{x,1},z_{y,1}\left|k\in\left(-\frac{\pi}{(
\alpha+\beta) b},\frac{\pi}{(\alpha+\beta) b}\right]\right.\right\}, \\
&\text{GBZ}_2=\left\{z_{x,2},z_{y,2}\left|k'\in\left(-\frac{\pi}{(a+b) \alpha},\frac{\pi}{(a+b) \alpha}\right]\right.\right\},
 \end{split}
\end{equation}
\begin{equation}\label{eqS314}
\begin{split}
\{\bar E\}=\{\bar E_1\}\cup\{\bar E_2\},\\
 \end{split}
\end{equation}
where the forms of $z_{x(y),1(2)}$ can be found in Eq.(\ref{eqS37},\ref{eqS39}) and that of $E_{1,2}$ can be found in Eq.(\ref{eqS38},\ref{eqS310}). Unlike with the unconstrained GBZ, they agree excellently with numerics, as shown in the example in FIG.~\ref{fig:S32}. In other words, the effective surrogate OBC Hamiltonian of the case $\alpha/\beta>a/b$ takes the form
\begin{equation}\label{eqS314.5}
\bar H_\text{2D-red}=\bar H_\text{2D,III}=\bar H_1\oplus\bar H_2 \ ,
\end{equation}
where $\bar H_1$ and $\bar H_2$ are the Hamiltonian operators with corresponding spectra $\bar E_1$ and $\bar E_2$ (Eqs.~\ref{eqS38} and~\ref{eqS310}).

\begin{figure}[htb]
\centering
\includegraphics[width=5in]{FIGS3-2.pdf}
\caption{The constrained GBZs and corresponding spectra $\bar E$ of two illustrative cases of $H_\text{2D, III}$ with hopping lengths $\alpha/\beta > a/b$, as indicated on top of the panels. Top) Coincident spectral contributions $\bar E_1,\bar E_2$ from 1D GBZ sectors $1$ and $2$ of these 2D lattices, as defined in Eq.(\ref{eqS38},\ref{eqS310}). They agree well with the numerical $E_\text{OBC}$ which is here computed with $10\times 10$ lattice sites, with density of states comparison presented later.
Bottom) 1D GBZs $z_x,z_y$ Eq.(\ref{eqS37},\ref{eqS39}) for sectors $1$ and $2$, which traces loops or flattened loops in the complex plane, consistent with the conclusions from Table.~\ref{tab:sum}. The two GBZ sectors may look totally different from each other, but they combine to describe the system as a coherent whole. Parameters are  $t_1=t_1=t_3=1$ and the GBZs are plotted with $k,k'$ points at intervals of $\pi/300$.
}\label{fig:S32}
\end{figure}

 \item[$\blacktriangleright$] When  $\alpha/\beta=a/b$, the GBZ and spectrum are given by
 \begin{equation}\label{eqS315}
\begin{split}
\text{GBZ}=\left\{z_{x,1}^{\alpha+\beta}z_{y,1}^{a+b}=\frac{\sin(\beta bk)}{\sin(a\beta k)}\text{e}^{i(a\beta+\beta b)k}  \left|k\in\left(-\frac{\pi}{(a+b) \beta},\frac{\pi}{(a+b) \beta}\right]\right.\right\},
 \end{split}
\end{equation}
\begin{equation}\label{eqS316}
\begin{split}
\bar E(k)=t_1^{\frac{b}{a+b}}t_3^{\frac{a}{a+b}}\left( \frac{\sin (a+b)\beta k}{\sin a \beta k}\right)\left( \frac{\sin \beta b  k}{\sin a \beta k}\right)^{-\frac{b}{a+b}}\text{e}^{\frac{2i\pi a v'} {a+b}}\ ,\\
 \end{split}
\end{equation}
with excellent numerical agreement as shown in FIG.~\ref{fig:S33}.

\begin{figure}[htb]
\centering
\includegraphics[width=3.4in]{FIGS3-3.pdf}
\caption{
The constrained GBZs and corresponding spectra $\bar E$ of $H_\text{2D, III}$ of two illustrative cases with hopping lengths $\alpha/\beta=a/b$ as indicated in the panels. Top) Spectrum according to the GBZ (Eq.\eqref{eqS316}), which agree well with the numerical $E_\text{OBC}$ spectrum which is here computed with $20\times 20$ lattice sites.
Bottom) Corresponding constrained 1D GBZs $z_x,z_y$ (Eq.~\eqref{eqS315}) of these 2D lattices, with cusps that are not present in the GBZs of 1D non-Hermitian lattices. Parameters are  $t_1=t_1=t_3=1$ and the GBZ is plotted with $k$ points at intervals of $\pi/600$.
}\label{fig:S33}
\end{figure}

\item[$\blacktriangleright$]   When  $\alpha/\beta<a/b$, the energy is always zero and $z_{x}^{\alpha+\beta}\rightarrow \infty$, \ $z_{y}^{a+b}\rightarrow 0$, which means that the eigenstates always corner-localize, as illustrated in FIG.~\ref{fig:S34}. This is due to non-Bloch collapse, since hoppings $t_1$ and $t_3$ do not contain any net component that in the inverse direction of hopping $t_2$.
\end{itemize}

\begin{figure}[h]
\centering
\includegraphics[width=4.3in]{FIGS3-4.pdf}
\caption{(Top) Identically zero energies for the case of $\alpha\beta<a/b$ of the type III system Eq.\eqref{eqS31} ($ \bar E=0$ are from theoretical analysis , $E_{\text{OBC}}=0$ are numerical OBC eigenvalues). (Bottom) corresponding corner-localized distribution of eigenstates. The system parameters are $t_1=t_1=t_3=1$, with $\alpha,\beta,a,b$ indicated in the figures.}
\label{fig:S34}
\end{figure}

\subsubsection{Number of OBC spectral branches for type III hopping lattices}

To find the number of OBC spectral branches for $H_\text{2D,III}$ satisfies the symmetry
\begin{equation}\label{mastermodelsymmetry}
\begin{split}
&U_{v,v'}^{-1}H_\text{2D,III}U_{v,v'}=\exp\left(i\frac{2 \alpha \pi}{\alpha+\beta}v+i\frac{2a\pi}{a+b}v'\right)  H_\text{2D,III}\ , \\
&U_{v,v'}=\exp\left(i\frac{2 m\pi}{\alpha+\beta}v+i\frac{2n\pi}{a+b}v'\right)|m,n\rangle\langle m,n|\ ,
\end{split}
\end{equation}
with integer indices $v,v'$ ($v=1,2,..,\alpha+\beta,\,v'=1,2,...,a+b$), where $|m,n\rangle$ are real space basis orbitals. $U_{v,v'}$ is a unitary operator which satisfies
$U_{v,v'}U_{v,v'}^{\dagger}=\text{I}$ , (`I' represents the unit operator).  Assuming that the wave function $|\Psi^{n}\rangle $ is an eigenstate of the model Eq.\eqref{eqS31}   with eigenvalue $E_n$, the other eigenstates $U_{v,v'}|\Psi^{n}\rangle$ can be found via symmetry Eq.\eqref{mastermodelsymmetry} with eigenvalue $\exp\left(i\frac{2 \alpha \pi}{\alpha+\beta}v+i\frac{2a\pi}{a+b}v'\right) E_n$, that is,
  \begin{equation}
\begin{split}
H_\text{2D,III}|\Psi^{n}\rangle&=E_n|\Psi^{n}\rangle\ , \\
\exp\left(i\frac{2 \alpha \pi}{\alpha+\beta}v+i\frac{2a\pi}{a+b}v'\right)H_\text{2D,III}|\Psi^{n}\rangle&=\exp\left(i\frac{2 \alpha \pi}{\alpha+\beta}v+i\frac{2a\pi}{a+b}v'\right)E_n|\Psi^{n}\rangle\ ,\\
U_{v,v'}^{-1}H_\text{2D,III}U_{v,v'}|\Psi^{n}\rangle&=\exp\left(i\frac{2 \alpha \pi}{\alpha+\beta}v+i\frac{2a\pi}{a+b}v'\right)E_n|\Psi^{n}\rangle\ ,\\
H_\text{2D,III}U_{v,v'}|\Psi^{n}\rangle&=\exp\left(i\frac{2 \alpha \pi}{\alpha+\beta}v+i\frac{2a\pi}{a+b}v'\right)E_nU_{v,v'}|\Psi^{n}\rangle\ .\\
\end{split}
\end{equation}
Hence we see that the energy spectrum of $H_\text{2D,III}$ is parametrized by a 1-parameter family of states, with indices $v,v'$ taking $\alpha+\beta$,$a+b$ possible values. As such, there are
\begin{equation}
\text{LCM}\left[\displaystyle\frac{a+b}{\text{GCD}(a+b,a)},\frac{\alpha+\beta}{\text{GCD}(\alpha+\beta,\alpha)}\right]\ ,
\end{equation}
OBC spectral branches for type III lattices.

\clearpage
\subsubsection{Flowchart of our approach}
All in all, the procedure of obtaining the effective BZ and spectrum is illustrated in the following flowchart for generic 2D models. Higher-dimensional models can be dealt with analogously, and will be discussed later.

\begin{figure}[h]
\centering
\includegraphics[width=3in]{FIGS3-5.pdf}
\caption{Flowchart of our approach for generic 2D models. Starting with the original model which can be expressed as $E(z_x,z_y)$ by identifying $z_x=e^{ik_x}$, $z_y=e^{ik_y}$, we find how $z_x,z_y$ depend on $k_1,k_2$ by implementing x and y-OBCs separately. If the resulting spectrum does not contain nontrivial spectral winding, the effective BZ is still 2D and we are done. If not, we find 1D parametrizations $k$ of a submanifold of the $(k_1,k_2)$ torus such that the spectral windings are trivial. The effective 1D BZ comprises all such 1D paths. }
\label{fig:S35}
\end{figure}

\clearpage
\subsection{Density of states in the dimensionally-reduced effective system}
\label{sec33}

In the subsection above, we have shown that in order for the GBZ spectrum of $H_\text{2D,III}$ (Eq.~\ref{eqS31}) to exhibit nonzero winding for any closed path, there is a need to constrain the valid GBZ to a 1D subspace of the original naively obtained 2D GBZ. Here, we show that the numerically obtained OBC spectrum indeed correspond to equally spaced states in the 1D constrained BZ, just like how the momentum eigenstates are distributed in 1D BZ corresponding to a physically 1D system.

To verify that the spectrum of $H_\text{2D-red}$ indeed comprises equally-spaced ``momentum'' eigenstates, we examine the density of states (DOS) of the OBC spectrum. The DOS is given by
\begin{equation}\label{DOS}
\begin{split}
\rho (E)=\frac1{N}\frac{\delta n}{\delta E}\ ,
\end{split}
\end{equation}
which means that there are $\delta n$  eigenstates in the energy range $[E,E+\delta E)$ and  normalization constant $N$ is equal to the length of the energy to be measured. What we would like to check is the quantity $\sigma(k)$, which is the DOS in momentum space.
That is,  the momentum range $[k_i,k_i+\delta k_i)$ $(i=1,2) $ which corresponds to energy range $[E,E+\delta E)$ need to be considered, where $i=1,2$ refers to the possible sectors for the GBZs.
 If the states are indeed uniformly labeled by momentum, $\sigma_i(k)\propto \rho_i(E){\delta E_i}/{\delta k}$ should be uniform in each sector $i$, and $\sigma(k)=\sum_i\sigma_i(k)$ should only be proportional to the number of sectors corresponding to a particular value of $k$. For OBC states predicted by the GBZ sector energies $\bar E_i$, we have

\begin{equation}\label{eqS322}
\begin{split}
\sigma (k) &=\sum _i\sigma_i (k)=\frac{2\pi}{\mathcal{N}}\sum _i\frac{\delta n_i(k)}{\delta E_i(k)}\frac{\delta  E_i(k)}{\delta k} =\sum_i\delta n_i(k)\ ,\\
\end{split}
\end{equation}
with $\delta k=2\pi/ \mathcal{N}$.
For each energy interval $\delta E(k)= \bar E_i(k+\delta k)-\bar E_i(k)$,
one can obtain the number $\delta n(k)$ of occupied states within the interval.  With the knowledge of the
dependence between $\bar E_i$ of the GBZ and momentum $k$ , it's easy to get the DOS in momentum space $\sigma(k)$. This is shown in FIG.\ref{fig:S35}(a,b) for two illustrative sets of parameters.

\begin{comment}
\textcolor[rgb]{0,0.4,0.8}{For new function of $\sigma (k)$, (1) $\sigma_i(k)$ is  an integer,
$2\pi$ is equal to the length of the momentum to be measured $( \delta k=2\pi/\mathcal{N})$. $\sum_k \sigma (k) = \sum_k \sum _i\sigma_i (k)= \sum_k \sum _i n_i (k)$ is several times of  total no. of states, like FIG.\ref{fig:S35}(a), $E_1=E_2$,$\sum_k \sigma (k) =2\times$ total no. of states.\\
(2) we can not tell which energy $E_{\text{OBC}}$ belongs to $E_1$ or $E_{2}$.  $\sigma_1$,$\sigma_2$ can be seen as 2  independent channels.
For same momentum space $k-k+\delta k$,  there have energy $E_{\text{OBC}}$ both belong to $\delta E_1, \delta E_2$. We just add $\delta n_1$ and $\delta n_2$ to get $\sigma(k)$.\\
(3)the system have LCM$(a + b, \alpha+\beta)$ spectral branches, $ \sigma (k)$ must be LCM$(a + b, \alpha+\beta)\times n$ with integer $n$.}
\end{comment}

\subsubsection{Model with the GBZ sectors possessing the same energies}

For definiteness, we first consider $H_\text{2D,III}$ Eq.\eqref{eqS31} with $\alpha=2,\beta=1,a=1,b=2$, such that the energy takes the form
\begin{equation}\label{eqS318}
\begin{split}
E(z_x,z_y)=t_1z_x^{2}z_y+t_2z_x^{-1}z_y+t_3z_x^{-1}z_y^{-2},
\end{split}\end{equation}
where each of $z_x,z_y$ belongs to the 1D GBZ sectors $\text{GBZ}_1$ and $\text{GBZ}_2$, which from Sec.~\ref{sec32} Eq.(\ref{eqS313},~\ref{eqS314}) simplify to
\begin{equation}\label{eqS319}
\begin{split}
\text{GBZ}_1&=\left\{\left.z_{x,1}^3=\frac{t_2}{t_1}\text{e}^{ik}, z_{y,1}^3=\frac{t_3}{t_2}\frac{1}{2\cos \left(k/3\right)-1}\right|k\in[-\pi,\pi)\right\}\ ,\\
\text{GBZ}_2=&\left\{\left.z_{x,2}^3=\frac{t_2}{t_1}\left(2\cos \left(k'/3\right)-1\right),z_{y,2}^3=\frac{t_3}{t_2}\text{e}^{-ik'},\right|k'\in[-\pi,\pi)\right\}\ ,\\
&\qquad \qquad \text{GBZ}=\text{GBZ}_1\cup \text{GBZ}_2\ .
\end{split}
    \end{equation}
In other words, the set of GBZ eigenstates consists of eigenstates from both GBZ sectors. The OBC eigenenergies are approximated by GBZ energies, which are given either by $E(z_{x,1},z_{y,1})=\bar{E}_1(k)$ or   $E(z_{x,2},z_{y,2})=\bar{E}_2(k')$, where $k$ and $k'$ parametrize their respective GBZs. The GBZ Hamiltonian is hence given by
\begin{equation}\label{eqS320}
\begin{split}
H_\text{2D-red}&=\bar H_\text{2D,III}={H}_1\oplus {H}_2\ , \\
\bar{E}_1(k)=\bar{E}_2(k)&= t_1^{1/3}t_2^{1/3}t_3^{1/3} (2\cos (k/3)+1)^{1/3}(2\cos (2k/3)+1)^{2/3}\text{e}^{2i\pi v/3}\ ,
\end{split}
\end{equation}
with $k,k'\in(-\pi,\pi]$, $v=1,2,3$, see FIG.~\ref{fig:S35}(a). $ \bar E_1$ and $\bar E_2$ are the corresponding eigenenergies of $ H_1$ and $ H_2$; in this case, after substituting the sector-dependent forms of $z_{x,\nu}$ and $z_{y,\nu}$, they happen to take identical functional forms.

To justify the integrity of this GBZ construction, we first note that the spectral winding number is zero, since  Eq.\eqref{eqS320} is obviously real. Also, both the GBZ$_1$ and GBZ$_2$ are period in $k,k'$. The correctness of Eq.(\ref{eqS319},\ref{eqS320}) is demonstrated in FIG.~\ref{fig:S32}(b), which not only shows that the results of the GBZ and GBZ Hamiltonian can be trusted, but also that the relation between OBC energies and momentum $k$ are as expected.
\begin{figure}[htb!]
\centering
\includegraphics[width=4.5in]{FIGS4-1.pdf}
\caption{Density of states (DOS) in momentum space. (a,b) The solid curves denote the dependence of effective 1D GBZ sector energies $E_1,E_2$ with momentum $k$ for two possible models Eq.\eqref{eqS31.5}. Numerically, we can obtain $k\approx k+\delta k$ from each $\delta E(k)$ point from $E_\text{OBC}$ data (red circles, shown on the central horizontal), and obtain the momentum space DOS $\sigma(k)$. \\
(a1,b1) The momentum space DOS profile $\sigma(k)$ in our model Eq.\eqref{eqS31.5}, with $\delta k=\pi/8000$. For both models, the DOS (black dots) can be approximately fitted to a sinusoidal curve (black curve). Parameters are  $t_1=t_1=t_3=1$, with a OBC lattice of $L^2=10^2=100$ sites.
}\label{fig:S35}
\end{figure}

\subsubsection{Model with the GBZ sectors possessing different energies}

\indent  In the above, we considered $H_\text{2D,III}$ (Eq.\eqref{eqS31}) with  $\alpha=b$ and $\beta =a$, such that the two 1D GBZ sectors correspond to eigenenergies $\bar E_1$ and $\bar E_2$ Eq.\eqref{eqS320} with equivalent forms. That, however, is not the case for more generic models. Here, we consider an illustrative model with  $\alpha=2 $,  $\beta =1, a=1,b=1$ with dissimilar GBZ eigenenergies for its two sectors. We first write it in terms of $z_x,z_y$:
\begin{equation}\label{eqS3181}
\begin{split}
E(z_x,z_y)=t_1z_x^{2}z_y+t_2z_x^{-1}z_y+t_3z_x^{-1}z_y^{-1}\ ,
\end{split}\end{equation}
where each of $z_x,z_y$ belongs to either of the two GBZ sectors
\begin{equation}\label{eqS3191}
\begin{split}
\text{GBZ}_1&=\left\{\left.z_{x,1}^3=\frac{t_2}{t_1}2\cos(k/3)\text{e}^{ik}, z_{y,1}^3=\frac{t_3}{t_2}\frac{1}{2\cos \left(2k/3\right)+1}\right|k\in[-\pi,\pi)\right\}\ ,\\
\text{GBZ}_2=&\left\{\left.z_{x,2}^3=\frac{t_2}{t_1}\frac{\sin(3k'/4)}{\sin(k'/2)},z_{y,2}^3=\frac{t_3}{t_2}\frac{1}{2\cos \left(k/2\right)+1}\text{e}^{-ik'},\right|k'\in[-\pi,\pi)\right\}\ ,\\
&\qquad \qquad \text{GBZ}=\text{GBZ}_1\cup \text{GBZ}_2\ ,
\end{split}
    \end{equation}
as introduced in Sec.~\ref{sec32} Eq.(\ref{eqS313},~\ref{eqS314}).
GBZ eigenenergies $\bar E$ take the form of either $E(z_{x,1},z_{y,1})=\bar E_1(k)$ or $E(z_{x,2},z_{y,2})=\bar E_2(k')$, with their GBZ Hamiltonian and explicit eigenenergies
\begin{equation}\label{eqS3201}
\begin{split}
&\bar H_\text{2D-red}=\bar{H}_\text{2D,III}= \bar H_1\oplus \bar H_2 ,\\
\bar{E}_1(k)=\pm t_1^{1/3}t_2^{1/6}t_3^{1/2} &(2\cos (2k/3)+1)^{1/2}(2\cos (k/3))^{2/3}\text{e}^{2i\pi v/3}\ ,\\
\bar {E}_2(k')=\pm t_1^{1/3}t_2^{1/6}t_3^{1/2} &(2\cos (k'/2)+1)^{1/2}(2\cos (k'/2))^{2/3}(\sin(k'/4))^{2/3}\text{e}^{2i\pi v/3}\ ,
\end{split}
\end{equation}

with $k,k'\in(-\pi,\pi]$, $v=1,2,3,$ for the model Eq.\eqref{eqS3181}. As in the previous example, the spectral winding in Eq.\eqref{eqS3201} are both zero, and both GBZ$_1$ and GBZ$_2$ are periodic. However, here $\bar E_1(k)$ and $\bar E_2(k)$ manifestly take different functional forms. As presented in  FIG.~\ref{fig:S32}(a), the numerical results demonstrate the correctness of Eq.(\ref{eqS3191},\ref{eqS3201}), with both $\bar E_1(k)$ and $\bar E_2(k)$ agreeing with the distribution of numerically obtained $E_\text{OBC}$. Note that strictly, the above GBZ results are exact only in the continuum limit ($L\rightarrow\infty$) where the set of momenta $k$ tends towards an equally spaced set of points with $\sim L^{-1}$ separation.

\subsection{Circuit simulation of non-Hermitian lattices }

\indent Any electrical circuit network can be represented by a graph whose nodes and edges correspond to the circuit junctions and connecting wires/elements. The circuit
 behavior is fundamentally described by Kirchhoff’s law. Non-Hermiticity can be introduced in a RLC electrical circuit by means of negative impedance converters with current inversion (INICs)\cite{ningyuan2015time,lee2018topolectrical,hofmann2019chiral,stegmaier2021topological,circuitliu2021non,hofmann2020reciprocal,helbig2020generalized,ezawa2019non} (Fig.~\ref{fig:circuit3}(c1)). By means of Kirchhoff's law, it can be shown that~\cite{ningyuan2015time,lee2018topolectrical,hofmann2019chiral,stegmaier2021topological,imhof2018topolectrical,helbig2019band,circuitPhysRevB.100.184202,helbig2020generalized,lee2019imaging,circuitPhysRevLett.124.046401,hofmann2020reciprocal,lenggenhager2021electric,circuitliu2021non,zou2021observation,circuitliu202001,hohmann2022observation,zhang2022anomalous,bergholtz2021exceptional,ezawa2019non} each INIC possess the reduced 2-node Laplacian

\begin{equation}
\begin{split}
  \left(
    \begin{array}{c}
      I'_{\text{in}} \\
      I'_{\text{out}} \\
    \end{array}
  \right)
  =\frac{1}{2i\omega L}
  \left(
            \begin{array}{cc}
              1& -1 \\
              1 & -1 \\
            \end{array}
          \right)
          \left(
    \begin{array}{c}
      V_{\text{in}} \\
      V_{\text{out}} \\
    \end{array}
  \right).
  \end{split}
\end{equation}
If we further connect an INIC with an inductor of inductance $2L$ in parallel, we obtain
\begin{equation}
\begin{split}
  \left(
    \begin{array}{c}
      I_{\text{in}} \\
      I_{\text{out}} \\
    \end{array}
  \right)
  =\frac{1}{i\omega L}
  \left(
            \begin{array}{cc}
              1& -1 \\
              0 & 0 \\
            \end{array}
          \right)
          \left(
    \begin{array}{c}
      V_{\text{in}} \\
      V_{\text{out}} \\
    \end{array}
  \right).
  \end{split}
\end{equation}
for this parallel configuration pair. This INIC-inductor pair is very versatile, and can be used as the building block of arbitrary non-Hermitian circuit Laplacians. In general, circuit Laplacians $J$ are obtained via Kirchhoff rule $\bm{I} = J\bm{V}$ , where $I$ denotes the current input and the voltage $\bm{V}$ measures against ground at each node. As an initial step towards identifying circuits with tight-binding lattice models, we can get $\bm{I} = J\bm{V}$ in compact matrix form (Fig.~\ref{fig:circuit3}(c))
\begin{equation}
\begin{split}
  \left(
    \begin{array}{c}
      I_{1}^{m,n} \\
     I_{2}^{m,n} \\
     I_{3}^{m,n} \\
     I_{4}^{m,n} \\
    \end{array}
  \right)
  =
  \left(
            \begin{array}{cccc}
              -{1}/{(i\omega L_1)}& 0&0& {1}/{(i\omega L_1)}\\
              0 & -{1}/{(i\omega L_2)}& 0&{1}/{(i\omega L_2)} \\
            0&0  & -{1}/{(i\omega L_3)}&{1}/{(i\omega L_3)}\\
            0&0&0&i\omega C
            \end{array}
          \right)
          \left(
    \begin{array}{c}
      V_{m+2,n+1}\\
      V_{m-1,n+1} \\
      V_{m-1,n-2}\\
      V_{m,n}
    \end{array}
  \right).
  \end{split}
\end{equation}
For instance, the circuit in Fig.~\ref{fig:circuit3} (c), shown for just one unit cell, is mathematically described by a Laplacian matrix, which in momentum space under PBCs, takes the form
\begin{equation}\label{cceq5}
\begin{split}
&I(\bm{k})=J(\bm{k})V(\bm{k})\ ,\\
J(\bm{k})=\mu+&t_1\text{e}^{2ik_x+ik_y}+t_2\text{e}^{-ik_x+ik_y}+t_3\text{e}^{-ik_x-2ik_y}\ ,
\end{split}
\end{equation}
with $t_j=-1/i\omega L_j$, $j=1,2,3,$ and $\mu=i\omega C-\sum_{j=1,2,3}t_j$.
This takes the same form as a tight-binding Hamiltonian of the form (FIG.~\ref{fig:circuit3}(a))
\begin{equation}\label{cceq2}
\begin{split}
&H=\sum_{m,n}t_1|m,n\rangle\langle m+2,n+1|+t_2|m,n\rangle\langle m-1,n+1+t_3|m,n\rangle\langle m-1,n-2|\ ,
\end{split}
\end{equation}
with $H=t_1\text{e}^{2ik_x+ik_y}+t_2\text{e}^{-ik_x+ik_y}+t_3\text{e}^{-ik_x-2ik_y}=J(\bm{k})-\mu$ when expressed in momentum space.\\
\\
Notably, the system is sensitive to the boundary orientation, in that the effective lattice takes different forms under different boundary orientations, as shown in FIG.~\ref{fig:circuit3} (b1-b3), which correspond to the type I to III models discussed in this work. They correspond to their respective circuits in FIG.~\ref{fig:circuit3} (d1-d3).

 \begin{figure*}[htbp]
    \begin{centering}
    \includegraphics[width=0.8\linewidth]{FIGS4-2.pdf}
    \par\end{centering}
    \protect\caption{\label{fig:circuit3}
(a) The three hoppings (black,blue,red) that define the model of Eq.\eqref{cceq2}.
(b1-b3) Different double OBC boundary orientations on this model lattice gives rise to effective models with hoppings that are oriented differently (labeled as Parts I, II and III).
(c) To realize this model with an electrical circuit, each asymmetric hopping can be realized with an INIC $L_i$, $i=1,2,3$. (c1) gives the internal make-up of an INIC, with the operation amplifier (triangle) giving rise to the hopping asymmetry. (d1-d3) Circuit realizations of the model variations (b1) to (b3).
}
\end{figure*}

\section{2D model with topological zero modes from 1D topological invariant}\label{sec4}

Here we present a nontrivial implication of the effectively 1D GBZ of our 2D model of the class $H_\text{2D,III}$. Since the GBZ construction pertains to not just the spectrum, but represents a complex analytic deformation of the Hamiltonian \emph{operator} itself (to a so-called \emph{surrogate} Hamiltonian, see Ref.~\cite{lee2020unraveling}), it means that the OBC properties of such Hamiltonians are also nontrivially modified. In particular, an OBC feature of particular interest is the presence of topological zero modes. Below, we shall see how our model, which is defined in 2D, nevertheless host topological zero modes defined by a 1D topological invariant.

Noting that the GBZ construction is completely determined by the characteristic polynomial (bulk energy dispersion), we can construct a 2-band topological Hamiltonian with a 1D constrained GBZ by writing down a model with exactly the same eigenenergy dispersion as Eq.\eqref{eqS31.5} up to a constant offset $c$ that gaps the system, i.e. \begin{equation}\label{eqS41}
 H=\left(
      \begin{array}{cc}
        0 &H_{12}\\
        H_{21}& 0 \\
      \end{array}
    \right),
    \end{equation}
with eigenenergies $\mathcal{E}_{1,2}$, which are stipulated to satisfy $\left(\mathcal{E}_{1,2}\right)^2=H_{12}H_{21}=E+c=t_1z_x^{\alpha}z_y^{a}+t_2z_x^{-\beta}z_y^{a}+t_3z_x^{-\beta}z_y^{-b}+c$, and
\begin{equation}\label{eqS42}
\begin{split}
    H_{12}=\left(E+c\right)/Z&=\left(t_1z_x^{\alpha}z_y^{a}+t_2z_x^{-\beta}z_y^{a}+t_3z_x^{-\beta}z_y^{-b}+c\right)/Z\ ,\\
    &\quad\qquad H_{21}=Z\ ,
    \end{split}
    \end{equation}
where we allow $Z$ to assume either of the two forms for illustrative purposes: $Z=z^{a}_y$ or  $Z=z^{-\beta}_x$. $E$ is the energy function of Eq.\ref{eqS31.5}.
The gap induced by $c$ does not affect the GBZ solutions since energy degeneracies are not sensitive to a constant offset. \\

First, for the model Eq.\eqref{eqS31} with OBCs, either along one or both boundary directions, the spectrum of the GBZ hamiltonian must have a net zero winding number, even though the winding number of the off-diagonal term $H_{12}$ is dependent on the $Z$ term.  \\

With $\alpha/\beta>a/b$,  its GBZ Hamiltonian takes the form

\begin{equation}\label{eqS44}
\bar H=\bar H_1\oplus \bar H_2=\left(
      \begin{array}{cc}
        0 &\bar H_{12}\\
        \bar H_{21}& 0 \\
      \end{array}
    \right)=\left(
      \begin{array}{cc}
        0 &\bar H^1_{12}\\
        \bar H^1_{21}& 0 \\
      \end{array}
    \right)\oplus\left(
      \begin{array}{cc}
        0 &\bar H^2_{12}\\
        \bar H^2_{21}& 0 \\
      \end{array}
    \right)\ ,
    \end{equation}
		\begin{equation}
\bar H^i_{12}=(\bar E_i+c)/Z_i,\quad \bar H^i_{21}=Z_i\  .
\end{equation}
with $\bar E_{i}=\bar H^i_{12}\bar H^i_{21}$, ($i=1,2$) Eq.(\ref{eqS38},\ref{eqS310}).   Both $\bar H^i_{12}$, $i=1,2$  have the same forms of $z_{x,i},z_{y,i}$:
\begin{equation}
\begin{split}
&z_{x,1}^{\alpha+\beta} =\frac{t_2\sin((a+b) \beta k)}{t_1\sin((\alpha b-a\beta )k)}\text{e}^{i(\alpha b+\beta b)k}, \quad  z_{y,1}^{a+b} =\frac{t_3\sin(\beta b k)\sin((\alpha b- a \beta)k)}{t_2\sin( a \beta k)\sin((\alpha b+\beta b)k)}, \qquad k\in\left[-\frac{\pi}{(\alpha +\beta) b},\frac{\pi}{(\alpha +\beta) b}\right],\\
& z_{x,2}^{\alpha+\beta} =\frac{t_2\sin(\beta a k')\sin((\alpha a+ \alpha b)k')}{t_1\sin( a \alpha k')\sin((\alpha b-\beta a)k')}\ ,
\  z_{y,2}^{a+b} =\frac{t_3\sin((\alpha b-a\beta )k'))}{t_2\sin((\alpha a+a\beta )k')}\text{e}^{i(\alpha a+\alpha b)k'},\  k'\in\left[-\frac{\pi}{\alpha(a +b)},\frac{\pi}{\alpha(a +b)}\right]\ ,
 \end{split}
\end{equation}
For each GBZ sector $i$, we have $Z_i=z^a_{y,i}$ or $Z_i=z^{-\beta}_{x,i}$ depending on the choice of illustrative example.

Since we have designed this model such that it harbors nontrivial topological zero modes under double OBCs, from established results~\cite{Lee2019anatomy}, both $\bar H_{12}$ and $\bar H_{21}$ must have nonzero winding about $E=0$ (we shall suppress the GBZ sector label $i$ unless it is explicitly referred to), and these windings sum to zero.
For the off-diagonal term $\bar H_{12}$, the winding number about $E=0$
 \begin{equation}
   \begin{split}
   \omega(\bar H_{12})&=\oint\text{d}k\, \partial_k \log\bar H_{12}=\omega(\bar H^1_{12})+\omega(\bar H^2_{12})\neq 0\ ,\\
   \end{split}
\end{equation}
is also nonzero. We emphasize that although this nonzero winding criterion for topological modes was originally formulated for a 1D system, we have now used it for the 1D effective GBZ of a \emph{physically 2D} system. Some examples of cases with topological zero modes for double OBCs but not single OBCs are given in Fig.~\ref{fig:S41}.\\

However, this model with a single OBC exhibit very different results.
If we choose $Z=z^{a}_y$,  only when the $\bm x$ PBC, $\bm y$ OBC,  the system has  topological zero modes, the winding number of  $\bar H_{12}$  is non-zero on the single OBC GBZ.  By contrast, the same system with $\bm x$ OBC, $\bm y$ PBC has no topological zero modes (Fig.~\ref{fig:S41}). With double OBCs, the 1D GBZ for the double OBCs give rise to zero winding number for $\bar H^1_{12}$, and nonzero winding number for $\bar H^2_{12}$, both summing to a nonzero total winding of $\bar H_{2}$ from Eq.~\ref{eqS44}.

If we choose $Z=z^{-\beta}_x$ instead, the opposite is the case, the system with $\bm x$ OBC, $\bm y$ PBC have topological zero modes, the system with $\bm x$ PBC, $\bm y$ OBC does not have.
The winding number of  $\bar H^1_{12}$  is non-zero, the winding number of  $\bar H^2_{12}$  is zero in the GBZ Hamiltonian with double OBC. In an analogous way, the 2D model  with double OBC also  has  topological zero modes which is  contributed by $\bar H_{1}(k)$ from Eq.\ref{eqS44}.\\

We demonstrate the above arguments for two illustrative cases of $\alpha,\beta,a,b$ parameters.
We use $Z=z_y$ and examine the energy spectrum of our 2D two-band model Eq.\ref{eqS41} with different boundary conditions.  The observed numerical results prove the correctness of the above theoretically established results, in FIG.~\ref{fig:S41},\ref{fig:S43}. Specifically, it demonstrates that the construction of the effective 1D GBZ Hamiltonian can correctly predict the topological zero mode, despite the system being physically a 2D system.

 \begin{figure}[h!]
\centering
\includegraphics[width=7in]{FIGS4-3.pdf}
\caption{ Energy spectrum of the 2-band model (Eq.(\ref{eqS41})) with (a) single OBCs and (b) double OBCs, showing double OBC topological modes that originate from a 1D topological invariant. We have $Z=z_y$,  $\alpha=2,\beta=1,a=1,b=2$, and the constant is set to $c=5$ to make sure that the 2-band model is gapped under PBC. In (a), there is no topological zero mode, which is corroborated by the fact that $H_{12}$ in (a1) has zero winding. In (b), the presence of a double OBC topological zero mode agrees with the fact that for each of $H^i_{12}$ and $H^i_{21}$, the sum of their windings over $i=1,2$ is nonzero. The winding of $H_{12}^2$ is $\omega(H_{12}^2)=-1=-\omega(H_{21}^2)$  and $\omega(H_{12}^1)=\omega(H_{21}^1)=0$. Note that for the single OBC (xOBC, yPBC) case, the winding is about $k_x$, not $k_y$, so the winding of $H_{12}$ and $H_{21}$ is zero.
The plots are generated with (a) $k_y=-\pi:\pi/30:\pi$, $L=10$,  (b) with $L_x=L_y=10$.
}\label{fig:S41}.
\end{figure}

\begin{comment}\begin{figure}[h!]
\centering
\includegraphics[width=4.5in]{gfigures42.pdf}
\caption{ Energy spectrum of the 2-band model (Eq.(\ref{eqS41})) with different boundary conditions and $Z=z_y$ with a different set of parameters ($\alpha=2,\beta=1,a=1,b=1$), also showing double OBC topological modes that originate from a 1D topological invariant.  The constant is set to $c=5$ to make sure that the 2-band model is gapped under PBC,
(a) with $k_y=-\pi:\pi/30:\pi$, $L=10$, (b)  with $L_x=L_y=10$. \\
}
\label{fig:S42}
\end{figure}
\end{comment}

\begin{figure}[h!]
\centering
\includegraphics[width=3.5in]{FIGS4-4.pdf}
\caption{ Energy spectrum of the 2-band model (Eq.(\ref{eqS41})) with $Z=z_y$, under single and double open boundary conditions. The Left and Right columns denote  different sets of parameters (Marked in the figure). Evidently, under double OBCs, we observe topological zero modes that originate from a 1D topological invariant, indicating that the GBZ is an effective 1D GBZ, not a 2D GBZ. The constant offset in Eq.(\ref{eqS41}) is set to $c=5$ to make sure that the 2-band model is gapped under PBCs. The system is numerically discretized such that (a,b) $k_y=-\pi:\pi/30:\pi$, $L_x=10$, (c,d) $L_x=L_y=10$.}
\label{fig:S43}
\end{figure}

\section{generalizations to higher dimensions }\label{sec5}
Our approach for getting the correct 2D GBZs as well as their 1D constrained GBZs, if necessary, can be generalized to higher dimensions. In generic number of dimensions, the GBZ will first be symmetrically constructed in terms of relations between complex $z_x,z_y,z_z,...$ and real $k_1,k_2,k_3,...$. If there exist nonzero spectral windings, they will also have to be constrained to a lower dimension such that the all spectral windings in the constrained GBZ become zero.

Specializing to 3 dimensions for definiteness (higher dimension cases can be analogously written down), suppose we have
  \begin{equation}\label{eqS51}
H_\text{3D}=\sum_{m,n,l} \sum_{a,b,c}t_{abc}|m,n,l\rangle\langle m+a,n+b,l+c|\ ,
\end{equation}
with hopping amplitudes $t_{abc}$ and $a,b,c$ generic integers which represent hopping directions and distances. In the same way as before,  we insert the ansatz $|\psi\rangle\propto\sum_{m,n.l}\psi_{mnl}|m,n,l\rangle$ with $\psi_{m,n,l}=z_x^mz_y^nz_z^l$ into the bulk equations to obtain
\begin{equation}\label{eqS52}
\begin{split}
E(z_x,z_y,z_z)=\sum_{a,b,c}t_{abc}z_x^az_y^bz_z^c\ .
\end{split}
\end{equation}
Like before, there are generally many sets of $(z_x,z_y,z_z)$ that corresponds to a particular energy $\bar E=E(z_x,z_y,z_z)$.\\

To incorporate the boundary conditions, we can express a wavefunction $|\Psi\rangle$ with eigenenergy $\bar E =E(z_x,z_y, z_z)$ as different 1D chains in the x,y or z-directions, with two other momenta as parameters:
    \begin{equation}\label{eqS53}
    \begin{split}
|\Psi\rangle=|\Psi^1\rangle&=|\Psi^2\rangle=|\Psi^3\rangle=\ \propto \sum_{m,n,l}\Psi_{mnl}|m,n,l\rangle\ ,\\
&\Psi_{mnl}=\Psi^1_{mnl}=\Psi^2_{mnl}=\Psi^3_{mnl}\ ,\\
\Psi^1_{mn}&=\sum_{z_y,z_z} \left(\sum_{j}f_{z_{x,j}}z_{x,j}^m\right)g_{z_y}  z_y^nh_{z_z}z_z^l\ ,\\
\Psi^2_{mnl}&=\sum_{z_x,z_z} f_{z_x} z_x^m \left(\sum_{j'}{g}_{z_{y,j'}}z_{y,j'}^n\right)h_{z_z}z_z^l\ ,\\
\Psi^3_{mnl}&=\sum_{z_x,z_y} f_{z_x} z_x^m g_{z_y} z_y^n \left(\sum_{j''}h_{z_{z,j''}}z_{z,j''}^l\right).\\
\end{split}
\end{equation}
Here $|\Psi^1\rangle,\ |\Psi^2\rangle,\ |\Psi^3\rangle$ are different direction decompositions of the same wave function $|\Psi\rangle$,  and the first summation in the equation $\Psi^{j}_{mnl}$, $j=1,2,3$ means we must consider all the sets $z_x,z_y,z_z$ which have same energy $\bar E$.\\

Taking the  different manifestations of the wave function $|\Psi\rangle$ into the boundary conditions along different directions (i.e.,${\bm x, \bm y,\bm z}$ direction), like boundary conditions along ${\bm x}$ direction, we have  $\Psi_{m'{n}}=\Psi^1_{m'{n}}=0 $ with special value $m'$ and any values $n$, like `1D' system with paramters $z_y,z_z$.
Thus, for a particular fixed energy, we can in general write down a relation between $z_x$ and $z_y,z_z$, i.e.
\begin{equation}\label{eqS54}
\begin{split}
\mathscr{F}_1(z_x,z_y,z_z,k_1)=0\ ,
\end{split}
\end{equation}
with real $k_1$ ``momentum'' parametrizing one direction in the GBZ.
In the same way, we can get other  relations along ${\bm y},\bm z$ directions
\begin{equation}\label{eqS55}
\begin{split}
\mathscr{F}_2(z_x,z_y,z_z,k_2)=0\ ,
\end{split}
\end{equation}
\begin{equation}\label{eqS56}
\begin{split}
\mathscr{F}_3(z_x,z_y,z_z,k_3)=0\ ,
\end{split}
\end{equation}
 with real $k_2,k_3$. By construction, the $z_x,z_y,z_y$ satisfy all the boundary conditions. Hence we have 3 unknowns $z_x,z_y,z_z$ with 3 independent relational equation (Eq.~\ref{eqS54},~\ref{eqS55},\ref{eqS56}), which can in principle be simultaneously solved to yield relations
 \begin{equation}\label{eqS57}
\begin{split}
z_i=\mathcal{F}_i(k_1,k_2,k_3)\ ,
\end{split}
\end{equation}
with $i=x,y,z$. Taking $z_x,z_y,z_z$ Eq.\eqref{eqS57} into the energy equation Eq.\eqref{eqS52}, we can equivalently get the GBZ-predicted OBC energy
 \begin{equation}\label{eqS58}
\begin{split}
E(z_x,z_y,z_z)=\bar E(k_1,k_2,k_3)\ .
\end{split}
\end{equation}
If the winding number of $\bar E(k_1,k_2,k_3)$  Eq.\eqref{eqS52} is zero for all closed loops, the resultant GBZ of the 3D system will just be given by
\begin{equation}\label{eqS59}
\begin{split}
\text{GBZ}_\text{3D}=\{z_x,z_y,z_z|z_i=\mathcal{F}_i(k_1,k_2,k_3),i=x,y,z\}\ .
\end{split}
\end{equation}

And if the spectral winding number of $E(k_1,k_2,k_3)$  Eq.\eqref{eqS58} is non-zero, we must determine what constraints we need on the parameters $k_1,k_2$ and $k_3$ such that the spectral winding of all possible paths is zero, and that the GBZs defined by $z_x,z_y,z_z$ are periodic as we vary any of the remaining $k_i$. The explicit construction of such constraints to obtain the correct GBZ depends on the model, and will be the subject of future work - in the next subsection, we give a minimal example in 3D. In general, the constraints can either reduce the correct GBZ to a 2D subspace, or even a 1D subspace.

\subsubsection{Illustrative simple 3D model}
For a minimal example of the GBZ construction of a 3D model, we consider 2 hopping terms in 3D space i.e.\begin{equation}\label{eqS510}
H_\text{3D}=\sum_{m,n,l} t_{abc}|m,n,l\rangle\langle m+a,n+b,l+c|+t_{a'b'c'}|m,n,l\rangle\langle m-a',n-b',l-c'|\ ,
\end{equation}
with real $a,b,c,a',b',c'$. Assuming the wave function ansatz $|\psi\rangle\propto\sum_{m,n.l}\psi_{mnl}|m,n,l\rangle$ with $\psi_{m,n,l}=z_x^mz_y^nz_z^l$, we arrive at two possible cases. In the first case, $a/a'=b/b'=c/c'$ and previous derivations can be straightforwardly generalized to give
 \begin{equation}\label{eqS511}
\begin{split}
&\text{GBZ}_\text{3D}=\left\{\left. z_x^{a+a'}z_y^{b+b'}z_z^{c+c'}=\frac{t_{a'b'c'}\sin(a'k)}{t_{abc}\sin (ak)}\text{e}^{i(a+a')k}\right|k\in \left(-\frac{\pi}{a+a'}\ ,\frac{\pi}{a+a'}\right] \right \}\ ,\\
&\qquad E(z_x,z_y,z_z)=E(k)=\left(\frac{t^{a'}_{abc}t_{a'b'c'}^{a}}{(\sin (a k))^{a}(\sin (a' k))^{a'}}\right)^{\frac{1}{a+a'}} \sin \left((a+a')k\right) \text{e}^{\frac{2i\pi  av}{a+a'}}\ ,
\end{split}
\end{equation}
with $E(z_x,z_y,z_z)=t_{abc}z_x^{a}z_y^bz_z^c+t_{a'b'c'}z_x^{-a'}z_y^{-b'}z_z^{-c'}$, as shown in FIG.~\ref{fig:S51} (b,b1,c,c1).  In this model Eq.\eqref{eqS511} with  $a/a'=b/b'=c/c'$, the GBZ is only determined by a condition on the combination $z_x^{a+a'}z_y^{b+b'}z_z^{c+c'}$  rather than $z_x$, $z_y$ and $z_z$ separately, and is akin to a 1D model  along the $a\bm x+b\bm y+c\bm z$ direction, which is  consistent with the results of numerical diagonalization, as in FIG.~\ref{fig:S21}(a1---a3). \\

In the other case where $a/a'=b/b'=c/c'$ does not hold, there is non-Bloch collapse due to uncompensated hoppings in certain directions, and the GBZ and energy are simply given by
 \begin{equation}\label{eqS512}
\begin{split}
\text{GBZ}&=\left\{z_x^{a+a'}z_y^{b+b'}z_z^{c+c'}=-1 \right\}\ ,\\
&\qquad E(z_x,z_y,z_z)=0\ ,
\end{split}
\end{equation}
as shown in FIG.~\ref{fig:S51} (a,a1). In both of these cases, the GBZ spectrum have zero winding number, and further dimensional reduction is not necessary.\\

  \begin{figure}[htb]
\centering
\includegraphics[width=5in]{FIGS5-1.pdf}
\caption{Spectra and illustrative eigenstates for $H_\text{3D}$  with different $a,b,c,a',b',c'$ case. (a--c) Excellent agreement between the GBZ eigenenergies $E(z_x,z_y)$ of GBZ and OBC eigenenergies $E_{\text{OBC}}$ of Eq.\eqref{eqS510}.   (a1---c1) Spatial profiles for the $E_{\text{OBC}}$ eigenstate of a few other illustrative cases, clearly showing that the skin states are aligned along the ($a\bm x+b\bm y+c\bm z$) direction (b1,c1). The model parameters are $L_x=L_y=L_z=10$, $t_{abc}=1,t_{a'b'c'}=2$  and the specific values of $a,b,c,a',b',c'$ and  $E_{\text{OBC}}$ are indicated in the figure.}
\label{fig:S51}
\end{figure}
In general, in arbitrarily high dimensional systems, the number of independent $z_i$ is always the same as the  number of independent relations $\mathcal{F}_i$, with zero spectral winding either automatically satisfied or fulfilled via dimensional reduction to a lower-dimensional GBZ.

  \begin{figure}[htb]
\centering
\includegraphics[width=5in]{FIGS5-2.pdf}
\caption{Schematic illustration of effective BZ for higher dimensional systems.  (Top) In 2D lattices that are Hermitian or ``unentangled'', the BZ is simply a 2D torus; yet for ``entangled'' non-Hermitian cases, the effective BZ becomes the union of one or more 1D tori. (Bottom) Similarly, 3D lattices that are Hermitian or ``unentangled'' have 3D tori as their BZ, but for ``entangled'' non-Hermitian cases, the effective BZ is in general a union of various 1D and 2D tori. }
\label{fig:S61}
\end{figure}