\subsection{Experimental settings}

The $^{87}$Rb BEC is prepared in an optical dipole trap by evaporative cooling for $\sim 18~\text{s}$. The multiple discrete momentum states are coupled with multi-frequency Bragg laser pairs. The different frequency components are imprinted by two acoustic optical modulators. One shifts the frequency of the incoming beam by $-100~{\rm MHz}$, and another shifts it by $100~{\rm MHz}-\sum_n \nu_n/2\pi$ ($n\ge 1$) with $\nu_n = 4(2n-1)\hbar k^2/2m$ (see the main text). As a consequence, the transition between the two momentum states, $|n-1\rangle\leftrightarrow |n\rangle$, can be resonantly triggered by the $\{\omega_+,\omega_n\}$ laser pair.

After the system evolves for a finite time $\tau_e$, we directly resolve the populations in each momentum state by letting the atoms fall freely in space for 20 ms with all lasers switched off, before the atoms are imaged by a CCD camera. Atoms with different momenta get separated in the $x$-direction along which the Bragg beams are applied (see Fig.~1 in the main text). To obtain the relative populations in each state, we integrate the image in the $y$-direction, and then fit the data with a $10$-peak Gauss function, $\mathcal{A}(x)=\sum_{n=-1}^{8} A_n \text{exp}\left[-\left(\frac{x-nd}{a}\right)^2\right]$. Normalizing the resulting amplitude $\mathcal{A}_n$ by $\sum_{n=-1}^{8} \mathcal{A}_n$, we finally get the atom probability distribution on each site, $P_n$.

\subsection{Effective Hamiltonian with off-resonant terms}

Following the theory of light-atom interaction in Ref.~\cite{Gadway2015}, we obtain the effective time-dependent full Hamiltonian
\begin{equation}\label{eqs2}
 H_{\rm eff} = \sum_n\sum_i\hbar\frac{\Omega_+\Omega_i}{4|\Delta|}e^{i[(\omega_+-\omega_i-4(2n-1)\hbar k^2/2m)t+(\phi_+-\phi_i)]}|n-1\rangle\langle n| + {\rm H.c.}
\end{equation}
with $\phi_+$ and $\phi_i$ the phases of the $\omega_+$ beam and $\omega_i$ component respectively (see Fig.~1 in the main text). We simply let $\phi_+=0$, and $\phi_i$ be the modulated phase relative to $\phi_+$ from the AOM. As we choose $\omega_i = \omega_+ - 4(2i-1)\hbar k^2/2m$, the simplified ideal model can be obtained by considering only the resonant terms, i.e., letting $i=n$, as
\begin{equation}\label{eqs1}
 H^{(0)} = \sum_n \hbar t_n |n-1\rangle \langle n| + {\rm H.c.}
\end{equation}
with $t_n = e^{-i\phi_n}\Omega_+\Omega_n/4|\Delta|$. This gives the general tight-binding form for a momentum-state chain. If we simply treat the $n \ge 2$ part as an effective reservoir, and apply the second-order perturbation with $t_2=t_z$ and $t_{n\neq 2}=t$, the loss rate of site $|1\rangle$ should approximately be $\gamma_0 \sim t_z^2/t$ \cite{Lapp2019}. Then we obtain the dissipative two-level Hamiltonian in the main text.

\begin{figure*}[]
\includegraphics[width=0.5\textwidth]{figS1}
\caption{Population distribution after a $6$-state chain evolves for $0.5$ms. A comparison between the experimental results and those from numerical simulations shows the effect of the non-resonant Bragg diffractions. Simulations using the ideal Hamiltonian (\ref{eqs1}) and the full Hamiltonian adopt the same coupling parameters $t_n=2\pi\times 1.50(2)~{\rm kHz}$.}\label{figS1}
\end{figure*}

Then, clearly, the $\ell$th-order non-resonant terms, responsible for the transition $|n-1\rangle\leftrightarrow|n\rangle$, can be induced by $\{\omega_+, \omega_{n-\ell}\}$ and $\{\omega_+, \omega_{n+\ell}\}$ laser pairs with detunings of $\mp 8\ell\hbar k^2/2m$ respectively. These terms are given by
\begin{equation}
 H^{(\ell)} = \sum_n \hbar t_{n\pm \ell} e^{\pm i8\ell(\hbar k^2/2m)t} |n-1\rangle\langle n| + \text{H.c.},
\end{equation}
leading the full Hamiltonian $H_{\rm eff}=\sum_\ell H^{(\ell)}$. In our experiment, $8\hbar k^2/2m$ corresponds to $\sim 2\pi\times 16.2~\text{kHz}$.

To demonstrate the effect of non-resonant terms, we simulate a $0.5$ms time evolution with the time-dependent effective full Hamiltonian (\ref{eqs2}), and make a comparison to the experimentally measured populations, and to simulation results from the effective ideal Hamiltonian (\ref{eqs1}). As shown in Fig.~\ref{figS1}, the full Hamiltonian describes our system better, and the states $|-1\rangle$ and $|6\rangle$ indeed get populated due to the non-resonant couplings, even when they are not resonantly coupled.

\subsection{Quantum Zeno subspace}

\begin{figure}[]
\includegraphics[width=0.5\textwidth]{figS2}
\caption{\label{figS2}
Tailoring the Hilbert space with quantum Zeno effects. We work with the $8$-site momentum lattice illustrated in Fig.~\ref{fig2}(c) with the tunneling strength $t=2\pi\times 1.03(2)~{\rm kHz}$ except at locations of strong Zeno couplings: (a) $t_3=3t$, (b) $t_4=3t$ and (c) $t_5=3t$. The measured time evolutions of $P_n$ are confined within the left-hand subspaces, with small leakage to the outside due to finite Zeno coupling strengths. Here we fix the kick frequency $\Omega/t=10$, and $\tau_0/T = 1$ for all cases.
}
\end{figure}

In this section, we demonstrate the extendability of our experimental configuration by observing quantum Zeno effects withs an enlarged dissipative subspace. At the upper-left corner of the phase diagram in Fig.~\ref{fig4}, $\bar{\gamma}$ approaches zero, as loss to the reservoir is suppressed due to the quantum Zeno effect.
In this deep PTB regime, the Hilbert space is tailored to the quantum Zeno subspace and the system can be described by an approximately Hermitian Hamiltonian~\cite{Facchi2008}. Building upon this understanding, in Fig.~\ref{figS2}, we show a series of examples of tailoring the quantum Zeno subspace by enlarging the dissipative subspace.
Specifically, when we shift the Zeno coupling $t_z$ to $|2\rangle\leftrightarrow |3\rangle$, the effective dissipation rate for site $|2\rangle$ in Fig.~\ref{figS2}(a) should be $\gamma_0\sim 9t$. Now the oscillations take place only in a subspace $\mathcal{H}_z$ spanned by $\{|0\rangle, |1\rangle\}$, and the site $|2\rangle$ becomes a boundary as the quantum Zeno effect guarantees it maintaining its initial condition, $P_2=0$. The imperfect Rabi oscillation originate from several sources, for instance, finite strength of $\gamma_0$, high-order non-resonant coupling terms, and many-body interactions. The size of the Zeno subspace can be easily tuned by shifting the position of the strong Zeno coupling. As shown in Fig.~\ref{figS2}(b), we let $t_{4}/t=3$ and find that, the time evolutions are confined within the subspace $\mathcal{H}_z=\{|0\rangle, |1\rangle, |2\rangle\}$. Furthermore, the subspace gets enlarged with $|3\rangle$ included once we shift Zeno coupling to $|4\rangle\leftrightarrow |5\rangle$; see Fig.~\ref{figS2}(c).
These experiments suggest that, combined with the flexible control offered by momentum lattice,
quantum Zeno effects can be used to map the whole ensemble into multiple Zeno subspaces with arbitrary sizes.
The approximately coherent dynamics within the Zeno subspace provide opportunities to explore a variety of non-classical phenomena, for example, the generation of Sch\"odinger cat states~\cite{Barontini2015}. We also note that, while the multilevel systems studied in this section do not possess PT symmetry, the scalability of our system as well as the demonstration of quantum Zeno subspaces pave the way for future studies of the correspondence between PT symmetry and quantum Zeno effects in higher dimensions.








\clearpage
\end{widetext}
\clearpage
\end{document}

































