\section{supplemental materials}
\subsection{Polynomial fitting for non-Hermitian Anderson model and U(1) model}
We study the tight-binding model on a three-dimensional cubic lattice (Anderson model, AM),
\begin{align}
{\cal H}=\sum_{i}\varepsilon_i c^\dagger_{i}c_{i}+\sum_{\langle i,j \rangle}V_{i,j}c^\dagger_{i}c_{j}\label{AM}
\end{align}
and U(1) model,
\begin{align}
{\cal H}=\sum_i \varepsilon_i c_i^{\dagger} c_i+\sum_{\langle i,j \rangle} e^{2\pi i\cdot \theta_{i,j}}V_{i,j}c_i^{\dagger}c_j\label{U1}
\end{align}
where $c_i^\dagger$ ($c_i$) is the creation (annihilation) operator for electrons at site $i$ and $\varepsilon_i$ is random onsite potential. 
$\langle i,j\rangle$ means that $i$ and $j$ are the nearest neighbor lattice sites to each other, 
$V_{i,j}=V_{j,i}=1$ are the nearest neighbor hopping term, and $\theta_{i,j}=-\theta_{j,i}$ is the random phase 
that distributes uniformly within $[0,1)$. To study the Anderson transition (AT) in the 
non-Hermitian (NH) system, we set $\varepsilon_{j}=w_{j}^{r}+i w_{j}^i$ with the imaginary unit $i$. 
$w_{j}^{r}$ and $w_{j}^{i}$ are independent random numbers, both of which distribute uniformly within 
$[-W/2,W/2]$ for a given disorder strength $W$.

\begin{table}[b]
	\caption{Polynomial fitting results for (a) Anderson model and (b) U(1) model with non-Hermitian (NH) disorders. 
        The goodness of fit (GOF), critical disorder $W_{c}$, critical exponent $\nu$, 
		the scaling dimension of the least irrelevant scaling variable $-y$, and critical level spacing ratio $\langle r\rangle_c$ are shown for different 
		orders of the Taylor expansion of the scaling function $(m_1,n_1,m_2,n_2)$, for different system size range, and for 
        different disorder range. The square bracket stands for the 95\% confidence interval for each fitting result. Polynomial fittings with 
        various energy windows (`percent'), expansion orders, system size range, and disorder range have been carried out, for the purpose of proving 
        the stability of the fitting results against these changes.}
	\begin{tabular}{cccccccccccc}
		&&&&&&&&&&&\\
		\multicolumn{12}{l}{(a)NH Anderson model} \\
		\hline
percent &	$L$& $W$&	$m_1$& $n_1$ & $m_2$ & $n_2$ & GOF & $W_c$ & $\nu$ & $y$&$\langle r\rangle_c$   \\
\hline
10\%&	8-24& [6, 7.12]&   3&3&0&1&0.11&6.28[6.26, 6.30]&1.046[1.012, 1.086]&1.75[1.65, 1.84]&0.7169[0.7163, 0.7177]\\
10\%&	10-24& [6, 7.19]&  3&3&0&1&0.15&6.32[6.30, 6.34]&0.990[0.945, 1.040]&2.10[1.87, 2.35]&0.7155[0.7146, 0.7164]\\
10\%&	12-24& [5.9, 7.2]& 3&3&0&1&0.13&6.34[6.32, 6.36]&0.942[0.897, 0.989]&2.53[2.14, 2.90]&0.7145[0.7138, 0.7154]\\
\hline
5\%&	8-24& [6.14, 7.3]	&3&3&0&1&0.21&6.38[6.35, 6.40]&0.977[0.938, 1.022]&1.82[1.73, 1.90]&0.7157[0.7149, 0.7165]\\
5\%&	10-24& [6.14, 7.26] &3&3&0&1&0.17&6.37[6.33, 6.41]&0.948[0.878, 1.039]&1.77[1.54, 2.02]&0.7159[0.7143, 0.7179]\\
5\%&	12-24& [5.9, 7.3]   &3&3&0&1&0.12&6.42[6.39, 6.45]&0.825[0.765, 0.893]&2.27[1.89, 2.60]&0.7134[0.7124, 0.7151]\\
\hline
&&&&&&&&&&&\\
\multicolumn{12}{l}{(b)NH U(1) model} \\
\hline
percent  &	$L$& $W$&	$m_1$& $n_1$ & $m_2$ & $n_2$ & GOF & $W_c$ & $\nu$ & $y$ & $\langle r\rangle_c$  \\
\hline
10\% &	8-24& [7, 7.56]& 1&3&0&1&0.32&7.14[7.13, 7.15]&1.065[1.036, 1.100]&2.60[2.31, 2.89]&0.7178[0.7171, 0.7188]\\
10\% &	8-24& [7, 7.56]& 2&3&0&1&0.43&7.15[7.14, 7.16]&1.068[1.034, 1.105]&2.63[2.35, 2.92]&0.7177[0.7169, 0.7186]\\
10\% &	8-24& [7, 7.56]& 3&3&0&1&0.49&7.15[7.14, 7.16]&1.065[1.031, 1.103]&2.64[2.35, 2.92]&0.7177[0.7169, 0.7186]\\
10\% &	10-24& [6.8, 7.6]&3&3&0&1&0.12&7.14[7.12, 7.16]&1.091[1.050, 1.151]&2.50[1.88, 3.16]&0.7187[0.7170, 0.7201]\\
10\% &	12-24& [6.68, 7.64]&3&3&0&1&0.46&7.13[7.06, 7.16]&1.133[1.065, 1.411]&2.29[0.83, 4.62]&0.7187[0.7166, 0.7281]\\
\hline
5\% & 8-24 & [7, 7.8]	&2&3&0&1&0.09&7.23[7.22, 7.24]&1.028[1.000, 1.059]&2.37[2.18, 2.57]&0.7165[0.7155, 0.7176]\\
5\% & 8-24 & [7, 7.8]	&3&3&0&1&0.18&7.23[7.22, 7.24]&1.012[0.979, 1.048]&2.35[2.15, 2.54]&0.7167[0.7155, 0.7178]\\
5\% & 10-24 & [6.8, 8] &3&3&0&1&0.23&7.23[7.22, 7.25]&1.012[0.979, 1.047]&2.40[2.09, 2.73]&0.7164[0.7152, 0.7179]\\
5\% & 12-24 & [6.8, 8]	&3&3&0&1&0.23&7.24[7.22, 7.25]&1.027[0.974, 1.079]&2.68[2.09, 3.38]&0.7158[0.7143, 0.7177]\\
\hline
	\end{tabular}
	\label{table}
\end{table}

Level spacing ratio $r_i$ for each complex-valued eigenvalue $E_i$ as 
\begin{align}
	r_i\equiv  |z_i|,
\end{align}
with
\begin{align}
z_{i}\equiv  \frac{E_i-E_{\rm NN}}{E_i-E_{\rm NNN}}, \label{z_i}
\end{align} 
where $E_{\rm NN}$ and $E_{\rm NNN}$ are the nearest neighbor and next nearest neighbor to $E_i$ in the complex Euler plane. 
$z_i$ is a complex number with modulus and angle. Here we focus on the modulus of $z_i$ first. $r_i$ will be averaged within 
an energy window (see below) and then averaged over $M$ realizations of disordered systems. This gives a precise value of 
$\langle r\rangle$ with a standard deviation $\sigma^2_{\langle r\rangle}=\frac{1}{M-1}(\langle r^2\rangle-\langle r\rangle^2)$.

In order to have enough  energy levels whose critical $W$ for the AT are sufficiently close to that for $E=0$, 
we take only $10\%$ eigenvalues near $E=0$ in the Euler plane as the energy window, and calculate them for each disorder realization. Here, the periodic boundary condition is imposed in the 
three directions for both of the two NH models. We prepare $M$ disorder 
realizations so that the total number of eigenvalues for the statistics ($M\ \times10\% \times L^3$) 
reaches $5\times 10^7$ ($L<24$) or $10^7$ ($L=24$) for the NH AM and $10^7$ ($L\le 24$) for the NH U(1) model; 
$M=10^6$, $5\times 10^5$, $3\times 10^5$, $1.2\times 10^5,$ $6\times 10^4$, $10^4$ for $L=8,\!\ 10, \!\ 12, \!\ 16, \!\ 20, \!\ 24$  
in the AM, and $M=2\times 10^5$, $10^5$, $6\times 10^4$, $2.5\times 10^4$, $1.2\times 10^4$, 
$6.4\times 10^3$ for $L=8,\!\ 10,\!\ 12,\!\ 16,\!\ 20,\!\ 24$ in the U(1) model.

\begin{figure}[t]
	\centering
	\includegraphics[width=0.8\linewidth]{sfig1}
	\caption{Averaged level spacing ratio $\langle r\rangle$ as a function of disorder. The level spacing ratio $r$ calculated from 
		$10\%$ eigenvalues around $E=0$ is averaged over $M$ number of different disorder realizations. 
		We take $M=10^4$, $10^3$ for the system size $L=8$, $10$, respectively, in non-Hermitian (NH) Anderson 
		model (AM) and NH U(1) model. When $W$ is either very small or very large, 
		$\langle r\rangle $ approaches  constant values.}
	\label{r_constant}
\end{figure}

$\langle r\rangle$ goes to almost constant values when $W$ is  far away from the critical point (FIG. \ref{r_constant}).
Consider first the very strong disorder region; $W\gg W_c$. In the thermodynamic limit, energy levels are distributed randomly
without any correlations with others. Thus, $r$ could be any value within $0<r<1$ with  equal probability. Hence 
\begin{align}
	\langle r\rangle_{\rm insulator}=\int_{0}^{1}r\rho(r) dr=\frac{2}{3}, \label{r_insulator}
\end{align}
where the density for $r$ is determined from $\langle \rangle_{\rm insulator} =1$; $\rho(r)=2r$.
Both models with the very strong disorder indeed show $\langle r\rangle\approx 0.66$ (FIG. \ref{r_constant}).

Consider next the metal phase; $W\ll W_c$. The energy levels in the metal phase have 
repulsive interactions with others, because of spatial overlaps between the extended eigenfunctions. Accordingly, the density 
of $r$ near zero will be smaller than in the insulator phase, hence
$\langle r\rangle_{\rm metal}>\langle r\rangle_{\rm insulator}$. 

For the Hermitian case (Gaussian ensemble), a mean value of $r_i\equiv s_i/s_{i-1}$, where $s_i \equiv E_i-E_{i-1}$ on  
the real axis with $\{E_i\}$ ordered, takes a certain constant value in the metal phase. The value depends on the symmetry class in the 
three Wigner-Dyson (WD) classes \cite{Atas13}. For the non-Hermitian case, we found that 
$\langle r\rangle$ reaches a constant value in the weaker disorder region 
($W\le W_c$), and the value increases with the system size for both models (FIG. \ref{r_constant}). For example, 
$\langle r\rangle=0.7182$, $0.7201$ for $L=12$, $20$ at $W=1$ for the NH AM and 
$\langle r\rangle=0.7307$, $0.7329$, $0.7353$, $0.7363$ for $L=10$, $12$, $16$, $20$ at $W=1$ for the NH U(1) model.  
By an extrapolation, we speculate $\langle r\rangle_{\rm metal}$ in the thermodynamic limit as 
$0.720$, and $0.736$ for the NH AM and U(1) models, respectively. Thus, the distributions of $r$
in the metal phases are different in the two NH models (FIG. \ref{Pr_NH_O1_U1_Metal}).




In order to characterize the Anderson transition, the scale-invariant quantity $\langle r\rangle$ is adopted. 
We estimate the critical exponent (CE) $\nu$ by polynomial fitting method \cite{Slevin14}, with 
various system size range, disorder range and expansion orders in the polynomial fitting.
Moreover, we narrow the energy window from the $10\%$ eigenvalues around $E=0$ to $5\%$ to see the stability of the polynomial fitting results. 
The results are shown in TABLE \ref{table} (a) for NH AM 
and (b) for NH U(1) model. 

CEs are consistent with each other for data sets with various system size range, disorder range, and 
expansion orders. This implies that our results are stable and precise. 
Moreover, CEs estimated from the 
$5\%$ eigenvalues energy window are consistent with that from the $10\%$ eigenvalues energy window, 
although they have a tendency to become smaller.
CEs estimated from the $10\%$ eigenvalues energy window is preferable because of the abundant energy levels. 
We conclude $\nu=0.99\pm 0.05$ for NH AM and $\nu=1.09\pm 0.05$ for NH U(1) model.

When the energy window is narrowed from 10\% to 5\% eigenvalues, the critical disorders in the both models 
tend to get larger; delocalized states at the band center are more robust against the disorder.
Besides, the CE for the NH AM changes to smaller values, when the data points for the smaller system 
sizes are excluded. The change of the CE in the NH AM becomes more prominent with 
the 5\% energy window than with the 10\% energy window. This means that the CE for the NH AM 
suffers from a systematic error by the choice of the energy windows. 
For the NH U(1) model, on the one hand, CEs extracted from the 5\% energy window stay nearly constant 
against the exclusion of the data from the smaller system sizes; 
the fitting of the NH U(1) model is much more stable than that of the NH AM.  

The difference of the stability in the fittings between the two models could be explained as follows.  
In the limit of the strongly localized phase, the level spacing ratio reaches the same value for the both models; 
$\langle r\rangle= 2/3$, while $\langle r\rangle$ in the limit of the delocalized phase goes to two 
different constant values in the two models respectively (FIG. \ref{r_constant}); there are two plateau regimes of 
$\langle r\rangle$ as a function of the disorder strength in these models. Note that data points near the plateau regimes 
should {\it not} be included for the scaling analysis, for they are likely outside the critical regime.
FIG. 1 in the main text shows that for the NH AM, the intersection of curves of 
$\langle r\rangle$ is quite close to the plateau value of $\langle r\rangle$ in the limit of the delocalized phase.    
Therefore, it is expected that a valid data range of the polynomial fitting in the NH AM becomes 
quite small in the side of the metal phase. On the other hand, the intersection for the NH U(1) model 
is relatively far away from the plateau value in the delocalized phase; 
the valid data range of the fitting becomes much wider in the NH U(1) model. 




\subsection{Spectral rigidity}
The spectral rigidity is defined by a number variance within an energy window 
\begin{align}
\Sigma_{2}\equiv \langle \delta N^2 \rangle=\langle (N-\langle N\rangle )^2 \rangle,
\end{align}
where $N$ is the number of eigenvalues within the same energy window, and 
$\langle \cdots \rangle$ stands for the average over $M$ disorder realizations. $M$ is 
typically on the order of or larger than $10^4$. The spectral compressibility $\chi$ can be extracted by 
\begin{align}
\chi\equiv \lim_{L \rightarrow \infty} 
\lim_{N\rightarrow\infty}\frac{d \Sigma_2(N)}{d\langle N\rangle}. 
\end{align}
Here, we set a circular energy window around $E=0$ as $\{|E|<E_{\rm bound}\}$ with $E_{\rm bound}>0$. To calculate the spectral compressibility, 
we decrease the energy window (reduce $E_{\rm bound}$) for a fixed (but sufficiently large) system size. 


Let us focus on the spectral compressibility at a critical point for the AT; $W=W_c$. 
From TABLE \ref{table}, We choose $W_c=6.3$ for the NH AM and $W_c=7.16$ for the NH U(1) model. We prepare
$M$ disorder realizations with $M=2\times 10^5$, $10^5$, $3\times 10^5$, $1.2\times 10^5$, $6\times10^4$, 
$10^4$ for $L=8$, $10$, $12$, $16$, $20$, $24$ for the NH AM  
and $N=2\times 10^5$, $10^5$, $6\times10^4$, $2.5\times 10^4$, $1.2\times10^4$, $6.4\times10^3$ for 
$L=8$, $10$, $12$, $16$, $20$, $24$ for the NH U(1) model.
In order to extract $\chi$, we change the energy window, namely $E_{\rm bound}$, within $10\%$ eigenvalues for a fixed system size $L$, and 
obtain various $\langle N\rangle$ and $\Sigma_2$ for each system size $L$. Then we carry out the linear fitting for 
$\Sigma_{2}$ vs. $\langle N\rangle$ for each $L$ (TABLE \ref{table_chi}).
From TABLE \ref{table_chi}, we found that $\chi$ thus obtained is stable against changing the 
system size $L$. We thus choose $\chi$ as of the largest system size in TABLE \ref{table_chi}; 
$\chi\approx0.46$ for the NH AM and $\chi\approx0.55$ for the NH U(1) model.

For the 3D Hermitian AM, the critical spectral compressibility $\chi$ has been already studied by others 
\cite{Ndawana02,Bogomolny11,Ghosh17}. For the comparison with the NH cases, we recalculate the same quantity for 
the Hermitian case for much smaller size system ($L=10$). The 10$\%$ eigenvalues around $E=0$ are calculated at the 
critical disorder strength of the AT ($W=W_c$), and the average number and number variance are taken over $10^5$ disorder 
realizations for the AM ($W_c=16.5$, see \cite{Slevin18}), and over $6400$ disorder realizations for the U(1) model ($W_c=18.8$, see \cite{Slevin16}).
We narrow the energy window, obtain the various $\langle N\rangle$ and $\Sigma_2$, and carry out the linear fitting 
for $\Sigma_{2}$ vs. $\langle N\rangle$. This gives $\chi\approx 0.28$ \cite{Bogomolny11} for the Hermitian AM 
and $\chi\approx 0.31$ for the Hermitian U(1) model at the critical point of the AT.

\begin{table}[b]
	\centering
	\setlength{\tabcolsep}{6.5mm}
	\caption{Linear fitting results of critical spectral compressibility 
    $\chi$ of various system size $L$ for non-Hermitian (NH) Anderson model and U(1) model. }
	\begin{tabular}{c|cccccc}
		\hline
		\diagbox{model}{$L$}& 8 & 10 & 12 & 16 & 20 & 24 \\
		\hline
		NH Anderson model & 0.4658 & 0.461 & 0.4593 & 0.4585 & 0.457 & 0.4592 \\
		NH U(1) model & 0.5439 & 0.5425 & 0.5456 & 0.5469 & 0.5397 & 0.5533 \\
		\hline
	\end{tabular}\label{table_chi}
\end{table}


\subsection{Level spacing distribution for Gaussian random matrix}
For the Hermitian case, 
the level spacing distribution $P(s)$ in a metal phase 
is well described by the Wigner-Dyson (WD) surmise \cite{Wigner51} in random matrix theory; 
\begin{align}
P_{\rm GOE}(s)&=\frac{\pi}{2}s e^{-\frac{\pi}{4}s^2},  \  \ \beta=1 \nonumber \\
P_{\rm GUE}(s)&=\frac{32}{\pi^2}s^2e^{-\frac{4}{\pi}s^2},  \  \ \beta=2 \nonumber \\
P_{\rm GSE}(s)&=\frac{2^{18}}{3^6\pi^3}s^4e^{-\frac{64}{9\pi}s^2},  \  \ \beta=4 
\end{align}
where $\beta$ is the power exponent for smaller $s$ region (the Dyson index), and $\beta=1,2,4$ for Gaussian 
Orthogonal ensemble (GOE),Gaussian Unitary ensemble (GUE), and Gaussian Symplectic ensemble (GSE), respectively.
$P(s)\propto s^{\beta}$ for the smaller $s$ region is caused by repulsive interactions between the energy levels in the metal phase. 

\begin{figure}[b]
	\centering
	\includegraphics[width=0.65\linewidth]{sfig2}
	\caption{Level spacing distribution $P(s)$ for the Hermitian case. 
		$W=5, 200$ for Hermitian Anderson model (AM), and $W=2, 200$ 
         for U(1) models, corresponding to metal and insulator phase 
        ($W_c=16.5$ for the AM, and $W_c=18.8$ for the U(1) model). 
         The $10\%$ eigenvalues around $E=0$ 
         are calculated over $10^5$ samples for Hermitian AM and U(1) models ($L=10$). All the eigenvalues of 
         the random matrix are calculated over $6\times 10^6$ samples with the matrix dimension
        $D=10$ for the GOE, GUE, and over $64000$ samples with the matrix dimension
        $D=1000$ for GSE. $P(s)$ with $s$ that is not unfolded are also tried for comparison and unfolding makes no difference here, 
    because density of states is almost constant in the region calculated.}
	\label{Ps_HermitianSystem}
\end{figure}

In insulator phase, $P(s)$ obeys the Poisson distribution,
\begin{align}
 P(s)=e^{-s},
\end{align}
since energy levels are uncorrelated with others. At the critical point of the AT, 
\begin{align}
P(s) \propto s^{\beta_c},
\end{align}
for smaller $s$ region and 
\begin{align}
P(s) \propto e^{-\alpha s}, 
\end{align}
for larger $s$ region.


For the purpose of the comparison, we recalculated the level spacing distribution $P(s)$ of the Hermitian WD random matrix.
We prepare the Hermitian random matrix $H$ as,
\begin{align}
	H=(A+A^{\dagger})/2, \label{dagger}
\end{align}
where $A$ is a random matrix with a restriction according to the symmetry.
For GOE, $H^*=H$, so $A$ is a real random matrix. 
For GUE, there is no restriction for $H$, so that $A$ is a complex random matrix.  
For GSE, $\Sigma_y H^* \Sigma_y=H$, with
\begin{align}
	\Sigma_y=
	\begin{pmatrix}
	 0& -i\\
	 i& 0\\
	\end{pmatrix},
\end{align}
 so that $A$ has the following structure, 
\begin{align}
    A=\begin{pmatrix}
    X & Y \\
    -Y^* & X^* \\
    \end{pmatrix}, \label{symplectic}
\end{align} 
where $X$, $Y$ are complex random matrices. Each element of the real random matrix 
is produced by the same Gaussian distribution, and  is independent. Real and imaginary parts of each element 
of the complex random matrix are produced by the same Gaussian distribution, and are independent. 
Eigenvalues in the GSE are doubly degenerate by the symmetry. We excluded this trivial degeneracy in the energy level 
statistics of $P(s)$
We calculate all the eigenvalues of the random matrix with the matrix dimension $D=10$
and take energy level statistics over $6\times 10^6$ 
different realizations of the random matrix for GOE and GUE, while energy level statistics of the random matrix
with the matrix dimension $D=1000$ is taken over $64000$ different realizations of the random matrix for GSE. For the unfolded energy levels in FIG.~\ref{Ps_HermitianSystem}, 
we first calculate an averaged density of state out of many samples. In terms of the 
averaged density of states $\rho(E)$, we map an energy level in each sample into an 
integrated density of states (IDOS); 
\begin{align}
{\rm IDOS}(E_i)=\int_{}^{E_i}\rho(E)dE. 
\end{align} 
The level spacing corrected by the averaged density of states, $s_i$, 
is given by a difference between the neighboring IDOSs,
\begin{align}
	s_i \equiv {\rm IDOS}(E_{i+1})-{\rm IDOS}(E_i).
\end{align}   


Our recalculation in FIG. \ref{Ps_HermitianSystem} reproduces the Dyson index;   
$\beta\approx 1$ for GOE, $\beta\approx 2$ for GUE, $\beta\approx 4$ for GSE. 



\subsection{Level spacing distribution for Anderson model and U(1) model with Hermitian disorder}
For comparison, we also (re)calculate the level spacing distribution for Hermitian 
AM [Eq. (\ref{AM})] and U(1) models [Eq. (\ref{U1})] with real random onsite potentials. 
We calculate $10\%$ of the whole eigenvalues near $E=0$ with the system size $L=10$ for 
every disorder realization, and take the statistics over $10^5$ different disorder realizations.
FIG. \ref{Ps_HermitianSystem} shows that $P(s)$ in metal phase thus obtained  
for AM and U(1) model are consistent with $P(s)$ for GOE and GUE, respectively. 
Moreover, $P(s)$ with and without unfolding are almost identical to each other in metal phase. 
$P(s)$ with and without unfolding are so close to each other,  because the density of states 
within the $10\%$ energy windows is nearly constant in energy in metal phase. 
The unfolding process produces an error in $P(s)$ for small $s$ region, breaking a linear relationship; $\ln P(s) \propto \ln s$. 
We do not use the unfolding process when focusing on the behaviors of $P(s)$ at small-$s$ region. 

At the critical point, $P(s) \propto s^{\beta_c}$ for small $s$ region, where critical power-law exponent $\beta_c$ takes almost 
the same exponent as the corresponding Dyson index $\beta$ for every WD classes; 
$\beta_c \approx \beta$ \cite{Kawarabayashi96}. $P(s) \propto e^{-\alpha s}$ for large $s$ region, where $\alpha$ takes  
almost the same value for the three WD classes; $\alpha = 1.8\pm 0.1$~\cite{Batsch96,Hofstetter96,Zharekeshev95}.
In our calculation (FIG. \ref{Ps_O1_U1_H_large_s}), $\beta_c\approx0.98$, and $\alpha\approx 1.79$ at the critical point 
of the Hermitian AM ($W_c=16.5$), and $\beta_c\approx2.01$ and $\alpha\approx 1.71$ at critical point of the Hermitian U(1) model 
($W_c=18.8$).

\begin{figure}[t]
	\centering
	\subfigure[Hermitian Anderson model, $W_c$=16.5;  $\beta_c \approx 0.98$]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig3_1}
		
		\end{minipage}%
	}%
	\subfigure[Hermitian U(1) model, at $W_c$=18.8;  $\beta_c \approx 2.01$]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig3_2}
		
		\end{minipage}%
	}%
	
	\subfigure[Hermitian Anderson model, at $W_c$=16.5; $\alpha\approx 1.79$]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig3_3}
		
		\end{minipage}%
	}%
	\subfigure[Hermitian U(1) model, $W_c$=18.8;  $\alpha\approx 1.71$]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig3_4}
		
		\end{minipage}%
	}%
	\caption{Level spacing distribution $P(s)$ for small $s$ ((a), (b); without unfolding), and for large $s$ ((c), (d); with unfolding). 
		The data comes from $10^5$ disorder realizations with $L=10$ for both models. The data are fitted by  
		$P(s)\propto s^{\beta_c}$ at small s and $P(s) \propto e^{-\alpha s}$ at large s. Red points are data excluded 
		for the linear fitting.}
	\label{Ps_O1_U1_H_large_s}
\end{figure}

\subsection{Level spacing distribution for Ginibre ensemble}
Ginibre ensembles are classes of ensembles for non-Hermitian random matrices~\cite{Ginibre65}, and therefore they might be useful for understanding the 
energy level statistics and the AT in a non-Hermitian disorder system.
According to the classification, there exist three kinds of the Ginibre ensembles;  
Ginibre Orthogonal ensemble (GinOE) ($H^*=H$), Ginibre Unitary ensemble (GinUE) (no restriction on $H$), and Ginibre Symplectic ensemble (GinSE) ($\Sigma_y H^* \Sigma_y=H$).  
The matrix $A$ for GOE, GUE, and GSE without Eq.~(\ref{dagger})  
corresponds to random matrix in GinOE, GinUE, and GinSE, respectively. 
GinOE, GinUE, GinSE correspond to the symmetry class AI, A, AII in the 
classification for the NH system. Because of the symmetry, 
eigenvalues in GinOE and GinSE come in pairs; 
$\{E_i, E_i^*\}$, and eigenvalues in the upper-half Euler 
plane are sufficient for the energy level statistics. The double degeneracy 
on the real axis needs to be excluded. We thus use only those eigenvalues 
whose imaginary parts are greater than 1, to determine the energy level 
statistics. Now the eigenvalues are complex number and the level 
spacing $s$ is defined by
\begin{align}
	s_i=|E_i-E_{\rm NN}|
\end{align}
where $E_{\rm NN}$ is the nearest neighbor for $E_i$.
Here the density of states in complex Euler plane is almost constant in the 
region calculated, so that we omit the unfolding process.
For small $s$, $P(s)$ of all these three kinds of the random matrix, GinOE, GinUE, and GinSE,  
obeys the same distribution (FIG. \ref{Ps_NH_insulator}(b)) with a cubic repulsion; $P(s) \propto s^3$~\cite{Hamazaki20}

\subsection{Level spacing distribution for Anderson model and U(1) model with non-Hermitian disorder}
Level spacing distribution are calculated for the NH AM (Eq. \ref{AM}) and U(1) models (Eq. \ref{U1}). 
In insulator, $P(s)$ takes the 2D Poisson distribution \cite{Grobe88},
\begin{align}
	P_{P}^{2D}(s)=\frac{\pi}{2}se^{-\pi s^2/4}.
\end{align}
To test this formula, we calculates $10\%$ eigenvalues around $E=0$ in the complex Euler plane 
for the NH AM and U(1) model with $L=12$ at $W=100$, where the eigenstates in both models 
are in insulator phase. We determine $P(s)$ out of the $10\%$ eigenvalues calculated over $6\times 10^4$ 
different disorder realizations. $P(s)$ thus determined takes the same 2D Poisson distribution in the 
both models (FIG. \ref{Ps_NH_insulator}(a)).

To calculate $P(s)$ in metal phase, we calculate the $10\%$ eigenvalues around $E=0$  
for NH AM and U(1) model with $L=16$ at $W=3$, where it is guaranteed that 
all the eigenstates within the $10\%$ circular energy window are in the metal phase. We find that 
$P(s)$ for NH U(1) model is consistent with that for GinOE, GinUE, GinSE (FIG. \ref{Ps_NH_insulator}(b)), 
but $P(s)$ for NH AM deviates from that for GinOE, GinUE, GinSE.
In small $s$ region, $P(s)$ behaves as $s^{\beta}$, where $\beta\approx 2.74$ 
for NH AM and $\beta\approx 2.93$ for NH U(1) model 
(FIG. \ref{Ps_O1_U1_NH_c}). This is consistent with ref. \cite{Hamazaki20} where the 
class AI$^{\dagger}$ shows unique $P(s)$, different from $P(s)$ for the 
class A (FIG. \ref{Ps_NH_insulator}(b)).

To calculate $P(s)$ at the critical point,  we calculate the $10\%$ eigenvalues around $E=0$  
for the NH AM and U(1) model at $W=W_c$ ($W_c \approx 6.3$ for NH AM and $W_c\approx 7.16$ for NH U(1) model). 
The system size $L$ and the number of disorder realizations $M$ are set as 
$M=2\times 10^5$, $10^5$, $3\times 10^5$, $1.2\times 10^5$, $6\times10^4$, $10^4$ for $L=8$, $10$, $12$, 
$16$, $20$, $24$ for the NH AM, and $M=2\times 10^5$, $10^5$, $6\times10^4$, $2.5\times 10^4$, $1.2\times10^4$, 
$6.4\times10^3$ for $L=8$, $10$, $12$, $16$, $20$, $24$ for the NH U(1) model. For each $L$, $P(s)$ is determined  
from the $10\%$ eigenvalues calculated over the $M$ different disorder realizations. We fit $P(s)$ thus obtained 
as $s^{\beta_c}$ for small $s$ and as $e^{-\alpha s}$ for large $s$. $\beta_c$ and $\alpha$ are estimated for 
each system size $L$ (FIG. \ref{Ps_O1_U1_NH_c} and TABLE \ref{table_ps_c}). We find that the fitted values of 
$\beta_c$ and $\alpha$ are robust against the change of system size; $\beta_c= 2.6\pm 0.05$ , 
$\alpha=5.0\pm0.1$ for NH AM and $\beta_c= 2.9\pm 0.05$ , $\alpha= 4.5\pm 0.1$ for NH U(1) model.

\begin{table}[h]
	\centering
	\setlength{\tabcolsep}{5.5mm}
	\caption{Linear fitting result of critical level spacing distribution $P(s)$ for small and 
large $s$ regions. The fitting by $e^{-\alpha s}$ for the large $s$ region gives $\alpha$ and 
the fitting by $s^{\beta_c}$ for the small $s$ region gives $\beta_c$.}
	\begin{tabular}{c|c|cccccc}
		\hline
		Model & \diagbox{quantity}{$L$} & 8 & 10 & 12 & 16 & 20 & 24 \\
		\hline
		NH AM & $\alpha$ & 4.75 & 4.8 & 4.9 & 4.9 & 5.0 & 5.06 \\
		NH AM & $\beta_c$ & 2.61 & 2.6 & 2.63 & 2.6 & 2.56 & 2.62 \\
	NH	U(1) model & $\alpha$ & 4.1 & 4.1 & 4.3 & 4.4 & 4.46 & 4.61 \\
	NH 	U(1) model & $\beta_c$ & 2.96 & 2.98 & 2.96 & 2.95 & 2.88 & 2.84 \\
		\hline
	\end{tabular}\label{table_ps_c}
\end{table}


\begin{figure}[h]
	\centering
				\subfigure[ $P(s)$ for insulator phase]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig4_1}
		
		\end{minipage}%
	}%
	\subfigure[ $P(s)$ for metal phase and Ginibre ensembles]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig4_2}
		
		\end{minipage}%
	}%
	\caption{Level spacing distribution $P(s)$ in (a) insulator phase ($W=100$) and (b) metal phase ($W=3$) of non-Hermitian (NH) Anderson model (AM) and U(1) model.
		In (b)$P(s)$'s for GinUE, GinOE, GinSE with matrix dimension $D=10^4$ are also shown; 64 samples for GinUE and 640 samples for GinOE and GinSE.
		The distributions for insulator phase are constructed out of the $10\%$ eigenvalues around $E=0$ over $6\times 10^4$ samples for  both models ($L=12$).
	The distributions for metal phase are calculated from $10\%$ eigenvalues around $E=0$ over $6400$ samples for both models ($L=16$). }
	\label{Ps_NH_insulator}
\end{figure}

 \begin{figure}
	\centering
			\subfigure[NH Anderson model at $W=3$ (metal phase); $\beta \approx  2.74$]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig5_1}
		
		\end{minipage}%
	}%
	\subfigure[NH U(1) model at $W=3$ (metal phase); $\beta \approx 2.93 $]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig5_2}
		
		\end{minipage}%
	}%

		\subfigure[NH Anderson model at $W_c=6.3$; $\beta_c \approx  2.62$]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig5_3}
		
		\end{minipage}%
	}%
	\subfigure[NH U(1) model at $W_c=7.16$; $\beta_c \approx 2.84 $]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig5_4}
		
		\end{minipage}%
	}%

	\subfigure[NH Anderson model at $W_c=6.3$; $\alpha\approx 5.06 $]{
	\begin{minipage}[t]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{sfig5_5}
	
	\end{minipage}%
}%
\subfigure[NH U(1) model at $W_c=7.16$; $\alpha\approx 4.61$]{
	\begin{minipage}[t]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{sfig5_6}
	
	\end{minipage}%
}%
	\caption{Level spacing distribution $P(s)$ at small $s$ [(a)-(d)] and large $s$ [(e), (f)]. The blue lines are fitting curves by
    $P(s)\propto s^{\beta_{c}}$ at small $s$ and $P(s)\propto e^{-\alpha s}$ at large $s$. 
    In metal phase [(a), (b); $W=3$], the data comes from $6400$ disorder realizations with $L=16$. 
    At the critical point [(c)-(f); $W=W_c$], we take $W_c=6.3$, $L=24$ and $10^4$ samples for the 
   non-Hermitian (NH) Anderson model, and $W_c=7.16$, $L=24$ and $6400$ samples for the NH U(1) model. 
   Red points are data excluded for the linear fitting. }
	\label{Ps_O1_U1_NH_c}
\end{figure}

\clearpage
\subsection{Level spacing ratio distribution for the Hermitian system}
For the Hermitian case, we consider a distribution of level spacing ratio $r$, that is defined by\cite{Oganesyan07}
\begin{align}
r_i\equiv {\rm min}\Big(\frac{E_{i+1}-E_i}{E_i-E_{i-1}}, \frac{E_{i}-E_{i-1}}{E_{i+1}-E_{i}}\Big). 
\end{align}
Here $\{E_i\}$ are ordered in the ascending order ($E_1<E_2<E_3\cdots$). 
For comparison, we calculate the level spacing ratio distribution, $P(r)$, for the random matrix in 
GOE, GUE and GSE and for the Hermitian AM, U(1) models (FIG. \ref{Pr_Hermitian}). 
A random matrix theory \cite{Atas13} tells that $P(r)$ in the metal phase is given by
\begin{align}
P(r)=\frac{1}{C_{\beta}}\frac{(r+r^2)^{\beta}}{(1+r+r^2)^{1+\frac{3}{2}\beta}}\Theta(1-r). \label{pr_metal}
\end{align}
Here $C_{\beta}$ is a constant, $\beta=1,2,4$ for GOE, GUE, and GSE, respectively, and $\Theta(x)$ is the 
Heaviside step function.

In insulator, $P(r)$ is given by \cite{Atas13}
\begin{align}
P(r)=\frac{2}{(1+r)^2}\Theta(1-r). \label{pr_insulator}
\end{align}
for all the three WD classes. 
FIG. \ref{Pr_Hermitian} shows that $P(r)$ in metal phase of the Hermitian AM and U(1) models are consistent 
with $P(r)$ of GOE and GUE, respectively. It also shows that $P(r)$ in the insulator phase has the same distribution as 
in Eq.~(\ref{pr_insulator}) for  both models.

\begin{figure}[b]
	\centering
	\includegraphics[width=0.7\linewidth]{sfig6}
	\caption{Level spacing ratio distribution $P(r)$ for GOE, GUE, GSE, the Hermitian Anderson model (AM) and U(1) model. 
We set the matrix dimension $D=10$ with $6\times 10^6$ realizations for GOE, GUE, and $D=1000$ with $64000$ realizations for GSE. 
For the Hermitian Anderson model and U(1) model, $P(r)$ is obtained from $10\%$ eigenvalues around $E=0$ calculated 
with the system size $L=10$ and $10^5$ samples.}
	\label{Pr_Hermitian}
\end{figure}

\subsection{Level spacing ratio distribution for non-Hermitian system}
For the non-Hermitian case, we can consider not only the level spacing ratio $r$ but also the angle of $z_i$,
\begin{align}
\theta_i\equiv {\arg} (z_i),
\end{align}
where $z_i$ is defined in Eq. (\ref{z_i}). To see distributions of $r$ and $\theta$ in metal phase, 
we calculate $10\%$ eigenvalues around $E=0$ in the complex 
Euler plane for NH AM and U(1) model at $W=3$. We take the statistics over 6400 disorder realizations, 
to obtain the level spacing ratio distributions, $P(r)$ and $P(\theta)$. 
$P(r)$ and $P(\theta)$ in metal phase of the NH U(1) model are consistent with those of the GinUE. 
However, $P(r)$ and $P(\theta)$ in metal phase of the NH AM behave quite differently from the GinUE. 
This supports the conclusion of Ref. \cite{Hamazaki20} that the metal phase of the 
NH AM and the metal phase in the NH U(1) model belong to two different universality classes. 
Noted that FIG. \ref{Pr_NH_O1_U1_Metal} also shows a small deviation between $P(\theta)$ of the NH U(1) model 
and $P(\theta)$ of the GinUE. We speculate that $\theta$ is more sensitive to the finite-system-size effect 
than $r$, because of boundary effects \cite{Sa20}.

FIG. \ref{Pr_Ptheta_NH_O1_U1_Metal_insulator} (a)-(d) show behaviors of $P(r)$ and $P(\theta)$ from metal phase 
to insulator phase in the NH AM and NH U(1) model.
$P(r)$ becomes linear in $r$ in the insulator phase for both models. This observation is 
consistent with that in Ref. \cite{Sa20}. On the other hand, $P(\theta)$ shows a small peak at $\theta=0$ in the 
insulator phase for both models (FIG. \ref{Pr_Ptheta_NH_O1_U1_Metal_insulator} (e)). This observation is different from Ref. \cite{Sa20}. 
We speculate that the small peak in $P(\theta)$ come from the boundary effect, as pointed out in Ref \cite{Sa20}. 
Namely, those $E_i$ around the boundary of the $10\%$ circular energy window have higher 
chance to give smaller $\theta_i$, because such $E_i$ is apt to find its nearest ($E_{\rm NN}$) and next nearest 
neighbor ($E_{\rm NNN}$) in the same direction (an inner direction of the circular window; toward $E=0$). For 
calculations with the smaller system size, these eigenvalues near the circular boundary have considerable effect, 
causing a small peak at $\theta=0$ in $P(\theta)$. To uphold this speculation, we also calculate all the  
eigenvalues of the NH U(1) model with $L=20$ at $W=100$, and take the statistics over 640 samples.  
$P(\theta)$ thus obtained is flat in $\theta$ as expected (FIG. \ref{Pr_Ptheta_NH_O1_U1_Metal_insulator}(f)).

FIG. \ref{Pr_Ptheta_NH_O1_U1_c} shows critical $P(r)$ and $P(\theta)$ for the two NH models. 
$P(\theta)$ shows some amount of size dependences as $|\theta|$ approaches $\pi$, where 
$P(\theta)$ for the large $\theta$ has a tendency to be larger for larger system size. In other 
words, $P(\theta)$ for the small $\theta$ tends to be smaller for the larger system. We speculate 
that this size dependence also partially comes from the boundary effect mentioned above. 
Note also that $P(r)$ and $P(\theta)$ are almost identical in the two models, except for 
small deviations observed in $P(r)$ at smaller $r$ (FIG. \ref{Pr_Ptheta_NH_O1_U1_c} (f)) and 
$P(\theta)$ at larger $\theta$ (FIG. \ref{Pr_Ptheta_NH_O1_U1_c} (e)). We conclude that 
it is hard to distinguish the two different universality classes in NH AM and NH U(1) models 
in terms of critical distributions of $P(r)$ and $P(\theta)$. 

\begin{figure}[h]
	\centering
		\subfigure[ $P(r)$]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig7_1}
		
		\end{minipage}%
	}%
	\subfigure[ $P(\theta)$]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig7_2}
		
		\end{minipage}%
	}%
	\caption{$P(r)$ and $P(\theta)$ for the Ginibre unitary ensemble, non-Hermitian (NH) Anderson model ($W=3$), 
        NH U(1) model ($W=3$). We take 64 samples for the GinUE, $10\%$ eigenvalues around $E=0$ over 
        6400 samples for the NH Anderson model and U(1) model. It is guaranteed that the $10\%$ eigenvalues 
        around $E=0$ in these two NH models are in the metal phase at $W=3$.}
	\label{Pr_NH_O1_U1_Metal}
\end{figure}

\begin{figure}
	\centering
	\subfigure[$P(r)$ for NH AM in metal phase ($W=1$), critical point ($W=6.3$) and insulator phase ($W=7.4, 100$)]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig8_1}
		
		\end{minipage}%
	}%
	\subfigure[$P(r)$ for NH U(1) model in metal phase ($W=1, 6.24$) and insulator phase ($W=8.72, 100$)]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig8_2}
		
		\end{minipage}%
	}%

\subfigure[$P(\theta)$  for NH AM]{
	\begin{minipage}[t]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{sfig8_3}
	
	\end{minipage}%
}%
\subfigure[$P(\theta)$  for NH U(1) model]{
	\begin{minipage}[t]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{sfig8_4}
	
	\end{minipage}%
}%

	\subfigure[$P(\theta)$  for NH AM and U(1) model]{
	\begin{minipage}[t]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{sfig8_5}
	
	\end{minipage}%
}%
	\subfigure[$P(\theta)$ for NH U(1) model with $L=20$, $W=100$, $100\%$ eigenvalues.]{
	\begin{minipage}[t]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{sfig8_6}
	
	\end{minipage}%
}%
	\caption{$P(r)$ and $P(\theta)$  from metal to insulator phase (a)-(d), and (e) comparison of $P(\theta)$ between non-Hermitian (NH) 
    Anderson model (AM) and U(1) model. We take $10\%$ eigenvalues around $E=0$ over $6\times 10^4$ samples (NH U(1) model), 
    $3\times 10^5$ samples (NH AM at $W=6.3, 7.4$), and $6\times 10^4$ samples (NH AM at $W=1,100$). (f) 
   $P(\theta)$ in insulator phase of the NH U(1) model with $100\%$ eigenvalues.}
	\label{Pr_Ptheta_NH_O1_U1_Metal_insulator}
\end{figure}

\begin{figure}
	\centering
	\subfigure[$P(r)$ of the NH AM at the critical point]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig9_1}
		
		\end{minipage}%
	}%
	\subfigure[$P(r)$ of the NH U(1) model at the critical point]{
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=1\linewidth]{sfig9_2}
		
		\end{minipage}%
	}%

	\subfigure[$P(\theta)$ of the NH AM at the critical point]{
	\begin{minipage}[t]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{sfig9_3}
	
	\end{minipage}%
}%
\subfigure[$P(\theta)$ of the NH U(1) model at the critical point]{
	\begin{minipage}[t]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{sfig9_4}
	
	\end{minipage}%
}%

	\subfigure[critical $P(\theta)$ of the NH AM and U(1) model]{
	\begin{minipage}[t]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{sfig9_5}
	
	\end{minipage}%
}%
\subfigure[critical $P(r)$ of the NH AM and U(1) model]{
	\begin{minipage}[t]{0.5\linewidth}
		\centering
		\includegraphics[width=1\linewidth]{sfig9_6}
	
	\end{minipage}%
}%
	\caption{ $P(r)$ and $P(\theta)$  at the critical point. The statistics for the non-Hermitian (NH) Anderson model (AM) are taken over 
    $M=1.2\times 10^5$, $6\times10^4$, $10^4$ samples for the system size $L=16$, $20$, $24$, respectively. The statistics for 
     NH U(1) model are taken over $M=2.5\times 10^4$, $1.2\times10^4$, $6400$ samples for the system size 
   $L=16$, $20$, $24$, respectively.}
	\label{Pr_Ptheta_NH_O1_U1_c}
\end{figure}
\end{widetext}

\end{document}