

\UseRawInputEncoding
\documentclass[letter,prl, amsmath,amssymb,
 reprint,%reprint,%
 author-year,%
superscriptaddress,
longbibliography
]{revtex4-2}
% \documentclass[aps,prl,reprint,groupedaddress]{revtex4-2}

% Packages
\usepackage{amsmath}
\usepackage{braket}
\usepackage{derivative}
\usepackage{siunitx}
\usepackage{tabularx}
\usepackage{graphicx}
\usepackage{multirow}

\setcounter{secnumdepth}{3}

\renewcommand{\thesection}{S\arabic{section}}
\renewcommand{\thesubsection}{\arabic{subsection}}

\makeatletter
\let\old@subsection\subsection
\renewcommand{\subsection}[1]{%
  \old@subsection{#1}% headings untouched
  \protected@edef\@currentlabel{\thesection.\arabic{subsection}}% references
}
\makeatother


\newcommand\Vatsan[1]{{[\color{blue} V: #1]}}
\newcommand\AhanaC[1]{{[\color{blue} AC: #1]}}
\newcommand{\AhanaW}[1]{\textcolor{red}{#1}}
\newcommand{\Shiv}[1]{{[\color{purple} SP: #1]}}
\newcommand\jp[1]{{[\color{magenta} #1]}}
\newcommand\rv[1]{{[\color{red} RV: #1]}}
\newcommand\jordan[1]{
{[\color{green} JH: #1]}}



\begin{document}



\title{Supplementary Material: Universal monitored dynamics in
multimode bosonic systems}

 \author{Shivam Patel}
  \affiliation{Department of Physics and Astronomy, 
Rutgers University, Piscataway, NJ 08854, USA}
\affiliation{Center for Materials Theory,  
Rutgers University, Piscataway, NJ 08854, USA}
\author{Catherine McCarthy}
\thanks{co-first author}
\affiliation{Department of Theoretical Physics, University of Geneva, 24 quai Ernest-Ansermet, 1211 Gen\`eve, Switzerland}
 \author{Ahana Chakraborty}
 \affiliation{Department of Physics and Astronomy, Louisiana State University, Baton Rouge, Louisiana 70803, USA}
  \affiliation{Department of Physics and Astronomy, 
Rutgers University, Piscataway, NJ 08854, USA}
\affiliation{Center for Materials Theory,  
Rutgers University, Piscataway, NJ 08854, USA}
\author{Jordan Huang}
  \affiliation{Department of Physics and Astronomy, 
Rutgers University, Piscataway, NJ 08854, USA}
\author{Thomas Dinapoli}
  \affiliation{Department of Physics and Astronomy, 
Rutgers University, Piscataway, NJ 08854, USA}
\author{Romain Vasseur}
\affiliation{Department of Theoretical Physics, University of Geneva, 24 quai Ernest-Ansermet, 1211 Gen\`eve, Switzerland}
\author{J. H. Pixley}\email{jed.pixley@physics.rutgers.edu}
  \affiliation{Department of Physics and Astronomy, 
Rutgers University, Piscataway, NJ 08854, USA}
\affiliation{Center for Materials Theory,  
Rutgers University, Piscataway, NJ 08854, USA}
\affiliation{Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA}
\author{Srivatsan Chakram}\email{schakram@physics.rutgers.edu}
  \affiliation{Department of Physics and Astronomy, 
Rutgers University, Piscataway, NJ 08854, USA}
\date{\today}



\definecolor{byzantine}{rgb}{0.74, 0.2, 0.64}
\newcommand{\catherine}[1]{{\color{byzantine} CM: #1}}

%\maketitle must follow title, authors, abstract, and keywords
\maketitle

% \tableofcontents
\section{Additional numerical results for monitored bosonic circuits} \label{sec:additional}

\subsection{Initial state preparation protocol} \label{sec:stateprep}

\begin{figure*}
    \centering
    \includegraphics[width=0.4\linewidth]{si-figs/checkerboard2.pdf}
    \hspace{0.025\linewidth}
    \includegraphics[width=0.4\linewidth]{si-figs/checkerboard0.pdf}
    \caption{Ancilla entanglement entropy $S_R$ for monitored dynamics with BSFP gates both with (left) and without (right) on-site Hubbard gates.  The initial state is prepared by beginning from the ``checkerboard'' state $\ket{\Psi} = \frac{1}{\sqrt 2} \big[ \ket{0101...} \ket{0} + \ket{1010...} \ket{1} \big]$ and scrambling with $2L$ layers of BSFP and on-site Hubbard gates, before running $2L$ monitored layers.  Both cases reproduce the results beginning from a Haar-random state shown in the main text for experimentally accessible system sizes.}
    \label{fig:checkerboard}
\end{figure*}


fig:numericsHere we demonstrate the validity of the state preparation protocol discussed in the main text.  The numerics in the main text mostly start from a Haar-random initial state, which is not efficiently experimentally preparable.  In the experiment we instead begin with a ``checkerboard'' initial state $\ket{\Psi} = \frac{1}{\sqrt 2} \big[ \ket{0101...} \ket{0} + \ket{1010...} \ket{1} \big]$ and then scramble with $2L$ layers of interacting gates without measurements before beginning the monitored dynamics for $2L$ layers of gates.  For monitored dynamics both with and without on-site Hubbard gates, the resulting $S_R$ curves show the same trends as the data in the main text (Fig. 2). 



\subsection{Learnability protocol and data for BSFP gates} \label{sec:learning}


\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{si-figs/learn0.pdf}
    
    \caption{\textbf{Learnability numerics for BSFP monitored dynamics.} Accuracy $A(p)$ of an exact classical decoder for a monitored circuit composed only of BSFP gates.  In line with the main text numerics, the data is consistent with a phase at high measurement rates in which the ``decodability time'' is $O(L)$.}
    \label{fig:learn0}
\end{figure}



In the main text, we proposed an experimental protocol to detect a MIPT by computing the entanglement entropy $S_R$ of an ancilla qubit globally coupled to an $L$-mode system. In the near term, post-selection is experimentally feasible because current circuit-QED implementations are limited to system sizes of at most $L \sim 10$ modes~\cite{huang2025fastsidebandcontrolweakly}; however, post-selection quickly becomes untenable with even modest increases in $L$. Recent works have sought to bypass this post-selection bottleneck by viewing the MIPT as a ``learnability'' transition accessible to an exact classical decoder~\cite{Barratt_2022, Agrawal_2024, singh2025, kim2025}. In this approach, the order parameter is the decoder’s accuracy in predicting which initial state generated a given measurement record~\cite{Barratt_2022}. Here we outline such a learnability protocol and present corresponding numerical results. This protocol also has the advantage of not requiring an ancilla qubit, thereby reducing potential sources of experimental error.

In this protocol, we choose two Haar-random initial states $\ket{\psi_0}$ and $\ket{\psi_1}$ such that $\braket{\psi_0|\psi_1}=0$. We then run a standard monitored circuit starting from $\ket{\psi_\alpha}$, where $\alpha\in\{0,1\}$ is chosen at random. After the circuit is executed, the measurement record $\vec m$---including the locations and outcomes of all measurements---along with the complete gate list and a full description of the states $\ket{\psi_0}$ and $\ket{\psi_1}$, is passed to a classical decoder. The decoder's task is to predict $\alpha$, \emph{i.e.}, which initial state generated the record $\vec m$.

To make this prediction, the decoder simulates the full monitored evolution starting from $\ket{\psi_0}$ and $\ket{\psi_1}$ to compute the Born probabilities $P(\vec m|\alpha)$ of obtaining the record $\vec m$ given the initial state $\ket{\psi_\alpha}$. It then predicts the label corresponding to the larger probability. Averaged over $N$ measurement records, the decoder accuracy is
\begin{equation}
A=\frac{1}{N}\sum_{\vec m}\max\{P(\vec m|0),\,P(\vec m|1)\}.
\end{equation}
The decoder accuracy undergoes a transition from an ``encoded'' phase to a ``decodable'' phase~\cite{Barratt_2022}.

While this protocol avoids the need for experimental post-selection, it is not truly scalable because the decoder must simulate the full evolution of the state, which remains exponentially costly in $L$. The main advantage of the learnability approach is therefore that the exponential difficulty is shifted from an experimental task to a numerical one.

For monitored dynamics with on-site Hubbard gates ($U=2$) and parity measurements, we numerically find a transition in the decoder accuracy consistent with $z=1$ near $p_c \approx 0.25$, as discussed in the main text. In addition, for purely BSFP dynamics ($U=0$), the decoder accuracy is consistent with the main-text trends: at high measurement rates, the system appears to enter a phase in which the ``decodability time'' (analogous to the purification time) scales linearly with system size (Fig.~\ref{fig:learn0}).




\subsection{Number measurements} \label{sec:den}

\begin{figure*}
    \centering
    \includegraphics[width=0.4\linewidth]{si-figs/du2.pdf}
    \hspace{0.01\linewidth}
    \includegraphics[width=0.4\linewidth]{si-figs/du0.pdf}
    \caption{\textbf{Ancilla entanglement entropy for monitored dynamics with number measurements.}  Data is shown for monitored dynamics with ($U=2$, left) and without ($U=0$, right) on-site Hubbard gates with density measurements for a circuit of depth $4L$.  The circuit with on-site Hubbard gates produces a phase transition at critical measurement rate $p_c \approx 0.1$; in contrast, there is no transition in the $U=0$ circuit and the ancilla qubit purifies.}
    \label{fig:den}
\end{figure*}

The main text presented results for parity measurements; here, we present additional results for measurements of photon number.  As was the case for parity measurements, the behavior of the ancilla entanglement entropy curves $S_R(p,t)$ depend on whether the monitored circuit is composed of only BSFP gates or not (Fig.~\ref{fig:den}).  For a circuit of only BSFP gates ($U=0$), the data in Fig.~\ref{fig:den} indicates that there is no MIPT with $z=1$ and that the ancilla purifies for any measurement rate $p>0$.  In contrast, for circuits with on-site Hubbard gates ($U=2$), there appears to be a transition with $z=1$ at the critical measurement rate $p_c \approx 0.1$.

It is interesting to compare the impact of turning on $U>0$ for parity and number measurements.  Despite the fact that neither measurement type features a transition when $U=0$, the ancilla is in a different phase depending on the measurement type-- the ancilla always purifies for density measurements but appears to be in either a purified or critical phase for parity measurements.  Additionally, while there is a transition that appears when $U \neq 0$ for both measurement types, the mechanism through which this occurs appears to be qualitatively different.  Since performing a number measurement locally collapses the state onto a single Fock basis state, the standard picture of MIPTs as a competition between unitary dynamics that hide information and measurements that reveal information is accurate for number measurements.  This intuition is consistent with the fact that for number measurements, turning on on-site Hubbard gates allows the unitary dynamics to hide information more efficiently, which results in the emergence of a mixed phase.  This qualitative picture is entirely different than that of the parity measurement, where increasing the scrambling power of the gates helps the purification process (see Sec. \ref{sec:scrambling}) and produces a phase with a $\mathcal O(1)$ purification time.

Performing measurements onto the mod-$n$ subspace for a generic choice of $n$ is natural to implement on circuit-QED platforms, and understanding the impact of tuning $n$ on the MIPT remains an interesting question.  In practice, it is difficult to experimentally or numerically realize system sizes large enough to meaningfully study cases apart from the parity ($n=2$) and density ($n=Q+1$) cases, and so future investigations of this subject will likely be largely theoretical.


\subsection{Bipartite entanglement entropy} \label{subsec:bipartite}


\begin{figure*}
    \centering
    \includegraphics[width=0.4\linewidth]{si-figs/half-cutS-0.pdf}
    \hspace{0.01\linewidth}
    \includegraphics[width=0.4\linewidth]{si-figs/half-cutS-2.pdf}
    \caption{\textbf{Bipartite entanglement entropy} Bipartite entanglement entropy $S$ for systems with (left) and without (right) on-site Hubbard gates shown for system sizes $L=4-14$.  Both cases appear to feature $S$ that scales logarithmically with system size at large $L$ for numerically-accessible system sizes (insets).}
    \label{fig:bipartite}
\end{figure*}


Here we briefly present results for the bipartite entanglement entropy $S = -\mathrm{Tr}[ \rho_A \ln \rho_A]$ with parity measurements for monitored dynamics with and without on-site Hubbard gates.  In contrast to the ancilla entanglement entropy data presented in the main text, the bipartite entanglement entropy $S$ looks similar for both $U=0$ and $U>0$ (\ref{fig:bipartite}).  Interestingly, the data for both values of $U$ appear consistent with a phase at high measurement rates in which $S \sim \log L$, although it is hard to accurately verify this scaling with only numerical access to small system sizes.

Considering the striking differences in the ancilla entanglement entropy data, it is interesting that the bipartite entanglement entropy looks similar for both cases.  The data is possibly consistent with a decoupling between the ancilla entanglement entropy $S_R$ for the $U=2$ , which purifies at high measurement rates, and the half-cut entanglement entropy, which appears to have a critical-like $\sim \log L$ scaling in the same regime.  Another possibility that we raise here is that this is a matter of protocol-- following the procedure in \cite{zabalo2020critical}, we calculate $S$ immediately after applying gates but before measuring.  Given that the parity measurements project onto an entire subspace, rather than individual Fock basis states, it is possible that choosing this protocol obscures the effect of the measurements.


\subsection{Alternative values of the dynamical exponent} \label{sec:alternative}

\begin{figure*}
    \centering
    \includegraphics[width=0.4\linewidth]{si-figs/z=1.66.pdf}
    \includegraphics[width=0.4\linewidth]{si-figs/z=1.4.pdf}
    \caption{\textbf{Alternative values of $z$.}  Alternative values of $(z,p_c)$ that yield an apparent collapse in the ancilla entanglement entropy $S_R(p_c,t/L^z)$ for BSFP gates with parity measurements.  Both $p_c=0.15, \ z=5/3$ (left) and $p_c=0.2, \ z=1.4$ (right) produce reasonable collapses-- without an independent prediction for either the dynamical exponent or the critical measurement rate, determining if there is a pair of $p_c,\ z$ from numerics alone is challenging.}
    \label{fig:zneq1}
\end{figure*}

In the main text, we presented numerical results based on the assumption that any MIPT would have the scaling exponent $z=1$.  The basis for this assumption is the theoretical argument for $z=1$ in the $d \rightarrow \infty$ limit with no symmetries \cite{Skinner_2019, Jian_2020} as well as numerical evidence showing $z=1$ in the qubit limit \cite{zabalo2020critical}.  The BSFP numerical data was consistent with the emergence of a novel critical phase at high measurement rates with a purification time that scales $\sim L$.  In principle, however, there is no reason to rule out the existence of a MIPT from a phase with an exponentially-long purification time to one with an $O(1)$ purification time with $z \neq 1$.  Furthermore, since there is no way to independently pinpoint the location of the transition with access only to the ancilla entanglement entropy, and since numerics are limited to small system sizes, there are multiple possible combinations of $(z,p_c)$ that are consistent with the numerics, two of which are illustrated in Fig. \ref{fig:zneq1} for $U=0$ (though the same holds true for the $U=2$ case as well).

Without an independent prediction of either $p_c$ or $z$, it is difficult to determine which, if any, combination of $(p_c \ z)$ corresponds to a true transition.  The approach taken in the main text is to assume $z=1$ and look for a transition consistent with this scaling exponent.  Regardless of  whether the system exhibits an anomalous transition with $z>1$ or if the system is in a phase with an $O(L)$ purification time at high measurement rates, the numerics demonstrate that the ancilla purifies faster when on-site Hubbard gates are used in addition to BSFP gates.



\subsection{Qualitative similarity of $z=1$ MIPT away from BSFP gates}\label{sec:qualitative}

\begin{figure}[h!]
    \centering
    \includegraphics[width=0.49\linewidth]{si-figs/u0.2.pdf}
    \includegraphics[width=0.49\linewidth]{si-figs/u2.pdf}
    \caption{\textbf{Stability of crossing for different $U$.} Data is shown for circuits of depth $t_f=4L$ with parity measurements for circuits of $U=0.2$ (left) and $U=2$ (right).  Tuning $U \neq 0$ does not appear to qualitatively affect the nature of the transition.}
    \label{fig:varyU}
\end{figure}


\begin{figure}[h!]
    \centering
    \includegraphics[width=0.49\linewidth]{si-figs/u0.pdf}
    \includegraphics[width=0.49\linewidth]{si-figs/bsrp.pdf}
    \caption{\textbf{Beam-splitter gates with fixed vs random phases.} Comparison of ancilla entanglement entropy $S_R(p, t \sim L)$ for BSFP (left) and BSRP (right) gates.  Choosing random phases produces a crossing in the $S_R$ curves.}
    \label{fig:bsrp}
\end{figure}

In the main text, we claim that there is a sharp distinction between monitored dynamics generated with purely BSFP gates and more generic circuits; here, we present numerical results validating this claim for two cases.  In particular, the presence of a crossing in the $S_R$ curves for the monitored dynamics generated by a combination of BSFP and on-site Hubbard gates remains when varying $U > 0$ (Fig. \ref{fig:varyU}); for a very small $0 < U \ll 1$, we would expect the crossing to emerge for sufficiently large circuit depths.  The precise mechanism through which this takes place remains an interesting open question.  Additionally, choosing beam-splitter unitaries with random phases (BSRP gates) restores a crossing, in comparison to the BSFP gates discussed in the main text (Fig. \ref{fig:bsrp}).



\section{Scrambling and purification with parity measurements} \label{sec:scrambling}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{si-figs/scramble.pdf}
    \caption{\textbf{Scrambling with parity measurements.} This figure illustrates a scenario in which increasing the scrambling power of the unitary dynamics can \emph{decrease} a system’s entanglement entropy. We consider a two-mode system with a conserved total photon number $Q=2$, initialized in a superposition within the even-parity subspace, and evolve it for one timestep under a unitary $U$. At $t=1$, we perform parity measurements on both modes. Evolving with the identity deterministically yields the measurement record $\vec{m}_{t=1}=\{0,0\}$, whereas evolving with a generic random unitary can yield $\vec{m}_{t=1}=\{1,1\}$, thereby reducing the entanglement entropy between the two modes.}
    \label{fig:scramble}
\end{figure}



As discussed in the main text, performing local measurements that project the state onto a subspace of the on-site Hilbert space, as opposed to onto a particular Fock basis state, leads to an effect in which scrambling helps the state purification process.  This effect originates from the fact that scrambling local degrees of freedom allows an observer to measure the same site multiple times and obtain different measurement outcomes, which results in a decrease in entanglement entropy.  Below, we will provide a concrete example of this effect taking place.

Consider a two-mode system with $Q=2$ photons.  The global two-mode state is spanned by the states $\{ \ket{20}, \ket{02}, \ket{11} \}$.  Consider time-evolving the initial state $\ket{\psi_{t=0}} = \frac{1}{\sqrt 2} \big[ \ket{20} + \ket{02}\big]$ by some unitary $U$ for one timestep before parity measurements are performed on both modes at $t=1$ (Fig.~\ref{fig:scramble}). For a generic random gate $U$, it is possible to either measure the outcome $\vec m_{t=1} =\{0,0\}$ or $\vec m_{t=1} =\{1,1\}$; if the outcome $\vec m_{t=1} = \{1,1\}$ is measured, the entanglement entropy between the two modes will decrease to zero.  However, choosing a less scrambling gate (e.g. the identity, or a constrained unitary that only allows for photons to hop in pairs) deterministically leads to the outcome $\vec m_{t=1} = \{0,0\}$ after the application of the gate.  In this case, the entanglement between the two modes will remain the same.  In contrast, for a situation in which number measurements are performed on both sites at both $t=0$ and $t=1$, the gate chosen does not impact the entanglement entropy between the two modes since the number measurements definitively collapse the state onto a single Fock basis state.

While somewhat contrived, the above example clearly illustrates that for measurements which project onto a subspace of the full Hilbert space and dynamics that mix between those sets of subspaces, increasing the scrambling power of the gate can lead to a decrease in the entanglement entropy in the system.  This intuition is consistent with the numerical results showing that BSFP and on-site Hubbard gates purify the ancilla faster than BSFP gates alone when using parity measurements.  Note that the argument above does not answer the interesting question of exactly how scrambling a gate has to be in order to produce a purified phase-- it only serves to illustrate the concept that increasing scrambling can favor purification in the presence of projective measurements onto subspaces.

\section{Analysis of the SNAIL coupler} \label{sec:AnalysisoftheSNAILcoupler}
In this section we derive the Hamiltonian of a Superconducting NonlineAr Inductive eLement (SNAIL) tunable coupler which we use to perform beam-splitter operations on our cavity modes.

A SNAIL consists of three identical Josephson junctions with Josephson energy $E_J$ in parallel with a smaller junction of energy $\alpha E_J$, forming a superconducting loop threaded by an external magnetic flux $\Phi_{\mathrm{ext}}$, as shown in Fig.~\ref{fig:SNAIL_schematic}. The superconducting phase across the smaller junction, $\varphi_s$, serves as a quantum degree of freedom, while the external flux tunes the circuit's effective nonlinearity and resonance frequency.

The potential energy of the SNAIL is 
\begin{equation}\label{SNAIL potential}
    U_S(\varphi_S) = -E_J[\alpha\cos{\varphi_S} + 3\cos{\frac{\varphi_s - \varphi_{\rm ext}}{3}}],
\end{equation}
where $\varphi_{\mathrm{ext}}=2\pi\frac{\Phi_\mathrm{ext}}{\Phi_0}$ and $\Phi_0$ is the flux quantum. The first term in the potential arises from the smaller junction while the second term corresponds to the three larger junctions. 

To form a tunable coupler, the SNAIL is connected in series with a capacitor of capacitance $C$ (charging energy $E_C=\frac{e^2}{2C}$), where $e$ is the electron charge, and a linear inductor of inductance $L$ (inductive energy $E_L=\frac{\varphi_0^2}{L}$ where $\varphi_0=\Phi_0/2\pi$). Denoting the phases across the capacitor, inductor, and SNAIL by $\varphi_C, \varphi_L$, and $\varphi_S$, respectively, Kirchhoff's loop rule gives 
\begin{equation}
    \varphi_c = \varphi_L + \varphi_S.
\end{equation}

The Lagrangian of the circuit is 
\begin{equation}
    \mathcal{L} = \frac{1}{2} C \varphi_0^2 \Dot{\varphi}_C^2 -
            \frac{1}{2} E_L \varphi_L^2 -
            U_S(\varphi_S).
\end{equation}
The absence of a kinetic term for $\varphi_S$ yields the current conservation equation,
\begin{equation}
E_L(\varphi_c - \varphi_S) = \frac{\partial U_S}{\partial \varphi_S},
\end{equation}
which implicitly defines $\varphi_S$ as a nonlinear function of $\varphi_C$.  

We treat the nonlinearity perturbatively by first linearizing the SNAIL's potential about its minimum, $U_S \approx \frac{1}{2} E_J c_2 \varphi_S^2$, where $c_2$ is the second order Taylor coefficient. Plugging this expression into the current conservation equation and performing some algebra gives
\begin{equation}
\varphi_s = p\varphi_c, \qquad
p = \frac{x_J}{c_2 + x_J},
\end{equation}
where $x_J=\frac{E_L}{E_J}$, and $p$ is the participation ratio of the SNAIL which quantifies the fraction of the total inductive energy contributed by the SNAIL. The remaining $1-p$ contribution comes from the linear inductor.


\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{si-figs/SNAIL_schematic.pdf}
    \caption{\textbf{Schematic of the SNAIL superconducting circuit used to activate beams-splitters between modes.} The SNAIL consists of three Josephson junction with junction energy $E_J$ in parallel with a junction with smaller energy $\alpha E_J$. The loop is threaded with an external flux $\Phi_\text{ext}$. We connect the SNAIL to a capacitor with charging energy $E_C$ and an inductor with inductive energy $E_L$ to make a SNAILmon circuit.}
    \label{fig:SNAIL_schematic}
\end{figure}

Using this linear approximation, we incorporate the SNAIL's nonlinearity by defining an effective nonlinear potential as a function of $\varphi_C$,
\begin{equation}
U^{\mathrm{NL}}(\varphi_c) =
\frac{1}{2} E_L (\varphi_c - \varphi_s)^2 + U_S^{\mathrm{NL}}(\varphi_s),
\end{equation}
where $U_S^{\mathrm{NL}}=U_S(\varphi_S)-\frac{1}{2}E_Jc_2\varphi_S^2$ is the nonlinear part of the SNAIL's potential. Expanding $U^{\mathrm{NL}}$ in powers of $\varphi_C$ gives effective Taylor coefficients $\tilde{c}_n$ that are related to those of the bare SNAIL potential,
\begin{equation}
    \tilde{c}_n = -x_J \frac{d^{n-1}\varphi_S}{d \varphi_C^{n-1}},
\end{equation}
providing a systematic way to dress the SNAIL nonlinearity by the linear inductor.

The Hamiltonian of the SNAILmon tunable coupler circuit is then
\begin{equation}
    H_{\mathrm{TC}} = 4E_CN^2 + E_J\sum_{n=2}^{\infty}\frac{\tilde{c}_n}{n!}\varphi^n,
\end{equation}
where $N=\frac{Q}{2e}$ is the dimensionless charge on the capacitor and we choose to drop the subscript from $\varphi_C$. We promote $\varphi$ and $N$ to operators and impose canonical quantization using $[\hat{\varphi},\hat{N}]=i$ and 
\begin{align}
   \hat {\varphi} &= \varphi_{\mathrm{zpf}}(\hat{c}^\dag + \hat{c}) \\
   \hat{N} &= i N_{\mathrm{zpf}}(\hat{c}^\dag - \hat{c}),
\end{align}
where $\varphi_{\mathrm{zpf}}$ and $N_{\mathrm{zpf}}$ are the zero-point fluctuations of the phase and charge, respectively, and $\hat{c}^\dag,\hat{c}$ are the creation and annihilation operators of the SNAILmon tunable coupler. Then the Hamiltonian takes the form
\begin{equation} 
    \hat{H}_{\mathrm{TC}} = \omega_c \hat{c}^\dag \hat{c} + \sum_{n=3}^{\infty} g_n (\hat{c}^\dag + \hat{c})^n,
    \label{eq:SNAILmon Hamiltonian}
\end{equation}
where $\omega_C=\sqrt{8 \tilde{c}_2 E_C E_J}$ is the resonant oscillator frequency of the tunable coupler circuit and the coefficients $g_n=\varphi^n_{\mathrm{zpf}}\frac{\tilde{c}_n}{n!}$.




\subsection{SNAIL-mediated beam-splitter interactions} \label{sec:Beam-splitter Hamiltonian}
The SNAILmon tunable coupler circuit is engineered to have a large third-order nonlinearity, $\frac{\tilde{c}_3}{\tilde{c}_n} \gg 1$, which leads strong three-wave mixing enabling fast high fidelity beam-splitter operations between cavity modes coupled to it \cite{chapman_high--off-ratio_2023,lu2023high, frattini2017}.

We consider a pair of cavity modes coupled to an externally driven SNAILmon tunable coupler. The system Hamiltonian is,
\begin{align}
    \hat{H}_{\mathrm{sys}} &= \hat{H}_{\mathrm{TC}} + \hat{H}_{\mathrm{cavity}} + \hat{H}_{\mathrm{coupling}} + \hat{H}_{\mathrm{drive}},\\
    \hat{H}_{\mathrm{cavity}} &= \omega_a \hat{a}^\dag \hat{a} + \omega_b \hat{b}^\dag \hat{b},\\
    \hat{H}_{\mathrm{coupling}} &= i \sum_{k\in\{a,b\}} g_{k} (\hat{k}^\dag+ \hat{k})(\hat{c}^\dag - \hat{c}),\\
    \hat{H}_{\mathrm{drive}} &= (\epsilon e^{-i(\omega_d t - \varphi)} - \epsilon^* e^{i(\omega_d t + \varphi)})(\hat{c}^\dag - \hat{c}).
\end{align}
$\hat{H}_{\mathrm{TC}}$ is given in Eq.~\ref{eq:SNAILmon Hamiltonian} (up to third order) with $\hat{c}^\dag$ / $\hat{c}$ being the tunable coupler's creation/annihilation operator. $\hat{H}_{\mathrm{cavity}}$ is the Hamiltonian of the cavity oscillator mode Alice(Bob), with frequency $\omega_{a(b)}$ and creation and annihilation operators $\hat{a}^\dag(\hat{b}^\dag)$ and $\hat{a}(\hat{b})$. $\hat{H}_{\mathrm{coupling}}$ describes the capacitive coupling of the tunable coupler to the cavity mode Alice(Bob) with coupling strength $g_{a(b)}$. Finally, $\hat{H}_{\mathrm{drive}}$ describes an external microwave tone of amplitude $\epsilon$ and frequency $\omega_d$ driving the charge on the SNAILmon capacitor, and phase $\varphi$. 

To see the beam-splitter operation in the interaction picture, we first diagonalize $\hat{H}_{\mathrm{coupling}}$ via a Schrieffer-Wolff transformation
\begin{equation}
    \hat{U}_S = e^{i\sum_{k \in a,b} \frac{g_k}{\Delta_k}(\hat{k}\hat{c}^\dag + \hat{k}^\dag\hat{c})},
\end{equation}
where $\Delta_k = \omega_c -\omega_k$. This transformation hybridizes the tunable coupler mode with Alice and Bob, $\hat{c} \rightarrow \hat{c} - i \sum_{k \in a,b} \frac{g_k}{\Delta_k} \hat{k}$. Next, we eliminate the oscillator terms and external drive by moving to a displaced, rotating frame using
\begin{align}
    \hat{U}_D &= e^{\tilde{\zeta}(t) \hat{c}^\dag - \tilde{\zeta}^*(t) \hat{c}},\\
    \hat{U}_R &= e^{-i(\omega_a \hat{a}^\dag \hat{a}
                +\omega_b \hat{b}^\dag \hat{b}
                +\omega_c \hat{c}^\dag \hat{c})t},
\end{align}
where $ \tilde{\zeta}(t) = \zeta e^{-i \omega_d t}, \zeta = \frac{\epsilon}{-\Delta_d} e^{i\varphi},$ and $\Delta_d = \omega_d - \omega_c$. After these transformations the dominant contribution arises from the cubic nonlinearity. Setting the drive frequency to the difference of the Alice and Bob modes, $\omega_d = \omega_a - \omega_b$, and taking a rotating-wave approximation yields the beam-splitter Hamiltonian,
\begin{equation}
    \hat{H}_{\mathrm{BS}} = g_{\mathrm{b s}}(e^{i\varphi}\hat{a}\hat{b}^\dag + \mathrm{h.c.}),
\end{equation}
where the beam-splitter rate $g_{\mathrm{b s}}=6 g_3 \frac{g_a}{\Delta_a} \frac{g_b}{\Delta_b} \frac{\epsilon}{\Delta_d}$. 

\section{Transmon-mediated Parity Measurement} \label{sec:Parity Measurement}
We describe how to perform a projective measurement of the photon number parity of a cavity mode using a transmon ancilla. Later, we generalize this protocol to enable full photon number measurements. 
The cavity mode is dispersively coupled to the transmon via
\begin{equation}\label{eq:Dispersive Hamiltonian}
    \hat{H}_\mathrm{disp} = \frac{\chi}{2} \hat{\sigma}_z \hat{n},
\end{equation}
where $\chi$ is the dispersive shift, $\hat{\sigma}_z$ is the Pauli-Z operator acting on the two lowest transmon states $\{\ket{g}, \ket{e}\}$, and $\hat{n}$ is the number operator of the cavity mode. This dispersive interaction causes the transmon transition frequency to shift depending on the number of photons in the cavity. As a result, the cavity's photon number information can be encoded in the phase of the transmon and subsequently readout using a Ramsey sequence. 

Take the initial state of the transmon, cavity system to be $\ket{\psi}=\ket{g,n}$, where $n$ is the photon number in the cavity mode. A resonant $\frac{\pi}{2}$-pulse on the transmon prepares it in a superposition, bringing the state to $\ket{\psi}=\frac{1}{\sqrt{2}}(\ket{g,n} + \ket{e,n}).$ The system is left to evolve freely under its dispersive interaction for a time $T$. During this interval, the transmon's Bloch vector precesses around the equator at a rate $\chi n$. This evolution imprints a photon-number-dependent relative phase between the transmon basis states,
\begin{equation}
    \ket{\psi(T)} =  \frac{1}{\sqrt{2}}(e^{i \frac{\chi}{2}  n T}\ket{g,n} + e^{-i \frac{\chi}{2}  n T}\ket{e,n}).
\end{equation}
A second $\frac{\pi}{2}$-pulse is then applied to the transmon after which the probability of the transmon being in $\ket{g}$ or $\ket{e}$ is given by
\begin{align}
    P_g = \sin^2 \left(\frac{\chi}{2} n T \right), \\
    P_e = \cos^2 \left(\frac{\chi}{2} n T \right).
\end{align}
If we set the idle time between $\frac{\pi}{2}-$pulses to $T=\frac{\pi}{\chi}$, the argument in the above expressions becomes $\frac{n\pi}{2}$. In this case reading out the transmon will deterministically project it to $\ket{e}$ or $\ket{g}$ for even or odd photon numbers, respectively. Thus, this Ramsey sequence performs a projective measurement of the cavity's photon number parity.

\section{Photon counting through repeated measurements of binary digits}\label{sec:Repeated Measurement Protocol For Photon Counting}
\label{sec: parity measurement}
\begin{figure*}\label{fig:Measurement protocol}
    \centering
    \includegraphics[width=\linewidth]{si-figs/Photon_count_protocol.png}
    \caption{\textbf{Protocol for measuring the photon number of an oscillator mode using a transmon.} To measure the photon number of a storage mode, the mode is swapped into the buffer cavity using the SNAILmon (blue line). The photon number is then determined bitwise in the binary representation using consecutive Ramsey sequences (dashed box). In each sequence, the wait time between the $\pi/2$ pulses is adjusted and the phase of the final pulse is set using feedforward from previous bit measurements, and the transmon is reset to its ground state before the next sequence. The Bloch spheres show the evolution of the transmon’s Bloch vector during measurements of the 0th and 1st bits, conditioned on whether the bit is 0 (green) or 1 (orange).}
    \label{fig:Photon_count}
\end{figure*}
In the previous section, we showed a Ramsey sequence which performs a parity measurement of the cavity's photon population, which corresponds to measuring the $0^{\mathrm{th}}$ bit in the binary expansion of $n$. We now generalize this protocol to measure up to an arbitrary bit of $n$, thereby enabling a fully projective photon number measurement \cite{Wang_2020}. 

To access higher order bits, we modify the idle time $T$ between the $\pi/2$-pulses  and by introducing a phase shift $\phi$ to the axis about which we drive the final $\pi/2$-pulse. With this modification, the probabilities to measure the transmon in $\ket{g}$ or $\ket{e}$ become
\begin{align}
    P_g = \sin^2 \left(\frac{\chi}{2} n T + \phi \right), \\
    P_e = \cos^2 \left(\frac{\chi}{2} n T + \phi \right).
\end{align}
As was the case for the parity measurement, the goal is to choose the idle time $T$ and phase $\phi$ such that the argument in the above expressions is an integer multiple of $\pi$, causing the measurement outcome of the transmon to be conditioned on the bit we wish to measure.

To see how this can be achieved, we consider the binary expansion of $n$,
\begin{equation}
    n = \sum_{i=0}^{K} b_i 2^i, \qquad b_i \in \{0,1\}.
\end{equation}
Suppose we wish to measure the $m$-th bit $b_m$. Dividing $n$ by $2^m$ separates the contribution of this bit from others: bits $>m$ make an even contribution to $\frac{n}{2^m}$ while all bits < $m$ make a total fractional contribution strictly less than one. Consequently, the parity of $\lfloor\frac{n}{2^m}\rfloor = \frac{n}{2^m} - \sum_{j=0}^{m-1} b_j \frac{2^j}{2^m}$ is determined by the value of $b_m$. 

In the transmon-cavity system, dividing by $2^m$ is implemented by setting the idle time $T=\frac{\pi}{2^m\chi}.$ The fractional contribution from the preceding bits is canceled by setting $\phi=-\frac{\tilde{n}\pi}{2^m}$, where $\tilde{n}$ denotes the photon number revealed by measuring all bits before the $m$-th. With these parameters, the measurement probabilities of the transmon state become,
\begin{align}
    P_g = \sin^2 \left(\frac{n - \Tilde{n}}{2^m} \pi \right), \\
    P_e = \cos^2 \left(\frac{n - \Tilde{n}}{2^m} \pi \right),
\end{align}
which deterministically map $\ket{e}$ to $b_m=0$ and $\ket{g}$ to $b_m=1$. 

Importantly, measuring the $m$-th bit requires prior measurements of all preceding bits since $\phi$ depends on $\tilde{n}$, necessitating real-time feedback. By performing these Ramsey sequences sequentially to measure all bits, we build up to a full photon number measurement. This repeated measurement protocol, shown in Fig.~\ref{fig:Photon_count}, is enabled by the long lifetimes of the cavity modes, which ensures photon loss remains minimal over the duration of the sequence.

\section{Preparing an Entangled State of the Reference Mode and the Multimode Cavity Register} \label{sec:Preparation Sequence For Entangled State of the Ancilla and the Multimode Cavity Register}
To probe the bosonic MIPT through ancilla purification, we prepare an initial state in which a dedicated ancillary storage mode is maximally entangled with two orthogonal half-filled configurations of the system storage modes: $\frac{1}{\sqrt{2}}(\ket{0101...}\ket{0} + \ket{1010...}\ket{1})$. This section details a gate sequence using universal sideband control on the RAQM architecture to construct this entangled "checkerboard" state ~\cite{huang2025fastsidebandcontrolweakly}.

We label the Hilbert space as $\ket{s_1s_2s_3..., \mathrm{ancilla}, \mathrm{transmon};ij}$, where $s_k$'s denote storage modes, the ancilla is a dedicated storage mode, the transmon is indexed by its three lowest energy states, $\ket{g}, \ket{e},\ket{f}$, and $i,j$ label the buffer cavity modes. Starting from $\ket{000...,0;g;0,0}$, we sequentially apply gates labeled as $\Theta^I_G$, where $\Theta$ is the pulse-angle, $G$ specifies the state manifold supported, and $I$ indicates the subsystems targeted: $q$ for transmon and $i,j$ for the buffer modes. For example, $\pi^i_{f0g1}$ denotes a $f0-g1$ sideband on buffer mode $i$, while $\frac{\pi}{2}^q_{ef}$ denotes a $\frac{\pi}{2}-$pulse between the $\ket{e}$ and $\ket{f}$ states of the transmon.

The first set of gates prepares an entangled superposition of single photon excitation states of the buffer modes to the $\ket{g}, \ket{e}$ states of the transmon,
\begin{equation}\label{eq:Shelving_protocol}
  \pi^j_{f0g1} \cdot \pi^q_{ge} \cdot \pi^q_{ef} \cdot \pi^i_{f0g1} \cdot \frac{\pi}{2}^q_{ef} \cdot \pi^q_{ge}
\end{equation}
Restricting attention to the transmon buffer subspace, $\ket{\mathrm{transmon; ij}}$, the state evolves as
\begin{align}\label{eq:State_evol_shelving}
    \ket{g;00} &\xrightarrow{\pi^q_{ge}} \ket{e;00} \notag \\
    &\xrightarrow{\frac{\pi}{2}^q_{ef}} \frac{1}{\sqrt{2}}(\ket{e;00} + \ket{f;00})  \notag\\ &\xrightarrow{\pi^i_{f0g1}} \frac{1}{\sqrt{2}}(\ket{e;00} + \ket{g;10})  \notag\\ 
    &\xrightarrow{\pi^q_{ef}} \frac{1}{\sqrt{2}}(\ket{f;00} + \ket{g;10})  \notag\\ 
    &\xrightarrow{\pi^q_{ge}} \frac{1}{\sqrt{2}}(\ket{f;00} + \ket{e;10})  \notag \\
    &\xrightarrow{\pi^j_{f0g1}} \frac{1}{\sqrt{2}}(\ket{g;01} + \ket{e;10}). 
\end{align}
Using the tunable coupler, we swap buffer mode $i$ into storage mode $s_1$ and buffer mode $j$ into $s_2$, resulting in the system state $\frac{1}{\sqrt2}(\ket{1000...000,0;e;00} +\ket{010...000,0;g;00})$.

To populate the next pair of storage modes, we perform a $\pi^q_{ef}$-pulse to bring the state to $\frac{1}{\sqrt{2}}(\ket{g;00} + \ket{f;00}$, then repeat the gate sequence in Eq.~\ref{eq:Shelving_protocol}. This brings our transmon-buffer cavity subsystem to the same final state in Eq.~\ref{eq:State_evol_shelving}. The repeated $\pi^q_{ef}$ and $\pi^q_{ge}$ pulses serve to "shelve" the $\ket{g}$ branch of the transmon's superposition in the $\ket{e}$ manifold, making it invisible to the proceeding $f0g1$ sideband, rendering it insensitive to off-resonant driving during the proceeding $f0g1$ sideband and preserving the coherence of the superposition state~\cite{huang2025fastsidebandcontrolweakly}. We again swap buffer mode $i$ into $s_3$ and buffer mode $j$ into $s_4$, resulting in the system state 
$\frac{1}{\sqrt2}\left(\ket{0101...000,0;g;00} + \ket{1010...000,0;e;00}\right)$.

Repeating this second set of gates populates the storage cavity in an alternating pattern, $\frac{1}{\sqrt{2}}\left(\ket{010101...10,0;g;00} +\ket{1010101...01,0;e;00}\right)$. To end the protocol we perform $\pi_{ef}-$pulse followed by a $\ket{f0}\!-\!\ket{g1}$ sideband to transfer the excitation in the $\ket{e}$ branch into buffer mode $i$, then swap the excitation into the ancilla mode resulting in the desired checkerboard state: $\frac{1}{\sqrt2}\left(\ket{010101...10,0;g;00} +\ket{1010101...01,1;g;00}\right)$.

We take the gate time for a $\pi_{ge}$-pulse to be 20 ns, $\pi_{ef}$-pulse to be 30 ns, and a $\pi_{f0g1}$ sideband to be 100 ns. The total time for state preparation is therefore
\begin{equation}
    T_\mathrm{state\,prep}=1.9\,\mu\text{s} + (L/2) \times 1.3\mu s,
\end{equation}
where $L$ is the system size. 



\section{Hubbard-like interactions via Selective number dependent arbitrary phase gates} \label{sec:Hubbard-like interactions via Selective number dependent arbitrary phase gates}
In this section we show how to achieve a Hubbard-like onsite interaction via selective number dependent arbitrary phase (SNAP) gates enabled by a transmon dispersively coupled to a cavity mode.

The dispersive interaction discussed in Sec.~\ref{sec: parity measurement} shifts the transmon's transition frequency depending on the cavity's photon population, allowing for independent control over the transmon within each photon number subspace. By driving the transmon at the photon number resolved frequencies and steering its Bloch vector along a closed path, an arbitrary phase can be imprinted on each Fock state,
\begin{equation}
    \hat{S} = \sum_n e^{i\theta_n} \ket{n}\bra{n}.
\end{equation}
The SNAP gate is performed in two stages. First, a selective multitone drive applies photon-number-dependent phases phase shifts to the state while exciting the transmon. Second, a fast, nonselective $\pi-$pulse returns the transmon to its ground state. 

Optimal selectivity favors longer pulse times, though longer pulse times increase incoherent errors, while faster times increase coherent errors. Recent work has shown that SNAP gates can instead be executed rapidly, on timescales $\mathcal{O}(\frac{1}{\chi})$, where coherent errors can be systematically suppressed through pulse optimization \cite{landgraf2024fastquantumcontrolcavities}. The multitone drive takes the form
\begin{equation}
    \Omega(t)
= \sum_n \lambda_n \, e^{i\left(\omega_n t + \alpha_n - \Delta \omega_n T / 2\right)},
\end{equation}
where $\lambda_n$, $\omega_n$ and $\alpha_n$ are the photon number ($n$) dependent drive amplitudes, frequencies, and phases, respectively, $\Delta \omega_n=\omega_n - \chi n$, and $T$ is the pulse duration. Starting from an initial guess, the drive parameters are iteratively updated according to 
 \begin{equation}
\Delta \lambda_n = -\frac{\epsilon_n^{(L)}}{2T}, 
\qquad
\Delta \alpha_n = -\Delta \theta_n, 
\qquad
\Delta \omega_n = \frac{\pi \epsilon_n^{(T)}}{2T}.
\end{equation}
Here $\Delta\theta_n$ is the phase error relative to the target $\theta_n$, while $\epsilon^{(L)}$ and $\epsilon^{(T)}$ quantify the longitudinal and transverse errors of the transmon's final Bloch vector from the desired state. These errors are extracted by projecting the Bloch vector onto a tangent plane to the $\ket{e}$ state on the Bloch sphere. Iterating the parameters according to the above rules perturbatively suppresses coherent errors while maintaining a fixed short gate time. We refer the reader to ~\cite{landgraf2024fastquantumcontrolcavities} for full details.

Choosing the SNAP phases as
\begin{equation}
    \theta_n=Un^2,
\end{equation}
 we realize the non-Gaussian Hubbard gates of our model. 



 
\section{Error Budget for the experiment} \label{sec:Error Budget for the experiment}

\begin{table*}[t]
\caption{Hamiltonian parameters used in Lindblad simulations.}
\label{tab:SM_hamiltonian_parameters}
\begin{ruledtabular}
\begin{tabular}{llll}
Subsystem & Parameter & Definition & Value \\ 
\hline

\multirow{4}{*}{Cavity--tunable coupler} 

& $g_{1}$ & Coupling strength between cavity mode 1 and coupler & 0.533 GHz \\
& $g_{2}$ & Coupling strength between cavity mode 2 and coupler & 0.592 GHz \\

& $\Delta_{1}$ & Detuning between cavity mode 1 and coupler & 2.261 GHz \\
& $\Delta_{2}$ & Detuning between cavity mode 2 and coupler & 2.590 GHz \\
& $g_{\mathrm{bs}}$ & Effective beam-splitter coupling rate & 0.500 MHz \\

\hline

\multirow{1}{*}{Cavity--transmon (SNAP)} 

& $\chi$ & Dispersive shift between cavity mode and transmon & -5.352 MHz \\

\multirow{1}{*}{Cavity--transmon (Parity)} 
& $\chi$ & Dispersive shift between cavity mode and transmon & -2.0 MHz \\

\end{tabular}
\end{ruledtabular}
\end{table*}

\begin{table*}[t]
\caption{Coherence time scales used in Lindblad simulations.}
\label{tab:SM_noise_parameters}
\begin{ruledtabular}
\begin{tabular}{llll}
Subsystem & Parameter & Definition & Value \\ 
\hline

\multirow{2}{*}{Cavity} 
& $T_1$ & Intrinsic relaxation time & 1.5(10) ms \\
& $n_a$ & Thermal occupation & 0.1\% \\

\hline

\multirow{9}{*}{Cavity--tunable coupler} 
& $T_1^C$ & Coupler energy relaxation time & $50$ $\mu s$ \\
& $T_\phi^C$ & Coupler pure dephasing time & $5$ $\mu s$  \\
& $n_C$ & Coupler thermal occupation & 2\% \\
& $T_{\mathrm{Th}}$ & Drive-induced thermal dephasing time & 2.5 ms \\
& $\tilde{T}_{1}^{\,a_1}$ & Effective relaxation time of cavity mode 1 & $562(824)$ $\mu s$ \\
& $\tilde{T}_{\phi}^{\,a_1}$ & Effective dephasing time of cavity mode 1 &  $7$ $\mu s$\\

& $\tilde{T}_{1}^{\,a_2}$ & Effective relaxation time of cavity mode 2 & $584(874)$ $\mu s$ \\
& $\tilde{T}_{\phi}^{\,a_2}$ & Effective dephasing time of cavity mode 2 & $92$ $\mu s$ \\

\hline

\multirow{3}{*}{Cavity--transmon qubit} 
& $T_1^q$ & Energy relaxation time & $200$ $\mu s$ \\
& $n_q$ & Thermal occupation & 0.5\% \\
& $T_\phi^q$ & Pure dephasing time & $100$ $\mu s$ \\

\end{tabular}
\end{ruledtabular}
\end{table*}

To model the impact of accumulation of errors, we simulate the full circuit by composing a noisy quantum error channel, which includes relevant dissipation and decoherence, for its elementary operations. Each error channel is obtained by numerically simulating Lindbladian evolution on the minimal target subsystem of the operation. We use the QuTiP python package~\cite{qutip5, johansson2013qutip, johansson2012qutip} to numerically time evolve the Lindblad master equation,
\begin{equation}
\label{eq:Lindbladian master equation}
\begin{split}
\dot{\rho}(t) =\;& -\frac{i}{\hbar} [\hat{H}(t), \rho(t)] \\
&+ \sum_n \frac{1}{2}
\left(
2 \hat{C}_n \rho(t) \hat{C}_n^\dagger
- \rho(t) \hat{C}_n^\dagger \hat{C}_n
- \hat{C}_n^\dagger \hat{C}_n \rho(t)
\right).
\end{split}
\end{equation}
where $\rho(t)$ is the system's density matrix at time $t$, $\hat{H}(t)$ is the time dependent Hamiltonian of the system, and $\hat{C}_n = \sqrt{\gamma_n}\hat{A}_n$ are collapse operators that describe coupling to the environment with strength $\gamma_n$.

The quantum error channel is a superoperator $U$ that maps an initial density matrix to a final density matrix under the corresponding noisy operation. A library of such superoperators for each operation is collected to compose a full end-to-end simulation. We consider four primary error channels corresponding to the main operations in our model. All relevant circuit and cavity parameters and coherence time scales used in simulation are included in Tables~\ref{tab:SM_hamiltonian_parameters},~\ref{tab:SM_noise_parameters}.


\subsection{Decay error} \label{subsec:Decay error}
As a baseline, we simulate the photon loss and thermal excitation of a single idle cavity mode, modeled as a truncated harmonic oscillator evolving freely for a time $T$. The decay and heating processes are modeled by $\sqrt{\kappa_a(1+n_a)} \hat{a}$ and $\sqrt{\kappa_an_a} \hat{a}^\dagger$, respectively, where $\kappa_a=1/T_1$ is the cavity's decay rate and $n_a$ is the average thermal occupation. Throughout, we take $T_1=1.5(10)$ ms and $n_a=0.001$.

 
\subsection{Beam-splitter gate error} \label{Beam-splitter gate error}
To model noisy BSFP operations, we simulate two cavity modes evolving under the beam-splitter Hamiltonian, with $\varphi=0$. 

Though we do not explicitly model the SNAILmon tunable coupler in our Lindblad simulations, we incorporate the inherited decay and decoherence of cavity modes from the coupler. The decay rate inherited by the cavity is 
\begin{equation}
    \tilde{\kappa}_a = \kappa_a + \left(\frac{g}{\Delta} \right)^2\kappa_C,
\end{equation}
where $g$ is the coupling strength between the cavity mode and the tunable coupler, $\Delta$ is the detuning between them, and $\kappa_C$ is the decay rate of the coupler. 
In the white noise regime, the inherited dephasing rate by the cavity mode is 
\begin{equation}
    \gamma_a \approx \left(\frac{g}{\Delta} \right)^4 \gamma_C,
\end{equation}
where $\gamma_C$ is the pure dephasing rate of the coupler. However, if the noise spectrum is pink and low frequency noise dominants, the inherited dephasing rate is instead,
\begin{equation}
    \gamma_a \approx \left(\frac{g}{\Delta} \right)^2 \gamma_C.
\end{equation}
As a moderate choice, we assume pink noise, consistent with experimental observations~\cite{chapman_high--off-ratio_2023}. In the presence of an external drive, the dephasing of the tunable coupler manifests as heating, inducing additional dephasing. This thermal dephasing rate is given by
\begin{equation}
    \gamma_{\mathrm{Th}} \approx n_C \kappa_C,
\end{equation}
where $n_C$ is the average thermal population of the coupler. We take the total inherited dephasing rate by cavity mode to be,
\begin{equation}
    \tilde{\gamma}_a = \left(\frac{g}{\Delta} \right)^2 \gamma_C + \gamma_{\mathrm{Th}},
\end{equation}
with the collapse operator taken to be $\hat{a}^\dag \hat{a}$.
We extract a noisy beam-splitter propagator with an average gate fidelity of $0.997(0)$.

Randomizing the phase $\varphi$ amounts to changing the phase of the external charge drive on the SNAIL and thus does not separate a spererate error channel. To obtain a noisy BSRP gate with phase $\alpha$, we simply apply the following unitary operation on the noisy BSFP gate,
\begin{equation}
    \hat{V}_\mathrm{BSRP}(\varphi=\alpha)=e^{i\alpha \hat{a}^\dag\hat{a}} \hat{V}_\mathrm{BSFP}(\varphi=0) e^{-i\alpha \hat{a}^\dag\hat{a}},
\end{equation}
 where $\hat{a}$ is the destruction operator for one of the target modes.

 
\subsection{Parity measurement error} \label{subsec:Parity measurement error}
Parity measurement errors are modeled by simulating a single cavity mode dispersively coupled to a transmon qubit. Intrinsic photon loss and thermal excitation for the cavity are included as before, and transmon decay, thermal excitation, and dephasing are modeled using the collapse operators $\sqrt{\kappa_q(1+n_q)} \hat{\sigma}_-,\sqrt{\kappa_q n_q} \hat{\sigma}_+$, and $\sqrt{\gamma}\hat{\sigma}_z$, respectively. We realize a parity measurement fidelity of $0.997(0)$. In addition to the intrinsic noise of the parity measurement, we model readout error by introducing a misassignment probability $\epsilon=0.004$, corresponding to the fraction of measurements in which the outcome is flipped and the state collapses onto the incorrect parity.

\subsection{SNAP gate error} \label{subsec:SNAP gate error}
The noisy SNAP gates are modeled using the cavity transmon system as in the parity measurements and the same noise parameters. Averaged over random cavity superposition states, we obtain a noisy SNAP gate fidelity of $0.997(0)$.

\subsection{Full circuit simulation} \label{subsec:Full circuit simulation}
Having obtained propagators for cavity dissipation, beam-splitter gates, SNAP (on-site Hubbard) gates, and parity measurements, we now simulate full monitored circuit realizations. This allows us to quantify how different noisy gates affect our ability to identify the transition. These simulations neglect correlated many-body dephasing errors and only include local decay and dephasing errors as described above.

The full circuit consists of four cavity modes, a dedicated ancilla mode, and a transmon qubit. The cavity modes are modeled as harmonic oscillators with a truncated local Hilbert space dimension of 3 so that at half-filling each oscillator can be occupied by a maximum of two photons. The ancilla and transmon are modeled as two-level systems. The total Hilbert space is therefore indexed as $\mathrm{mode}_1\, \otimes \, \mathrm{mode}_2\, \otimes \, \mathrm{mode}_3\, \otimes \, \mathrm{mode}_4\, \otimes \, \mathrm{transmon}\, \otimes \, \mathrm{ancilla}.$
The propagators obtained from each subsystem simulation are embedded into this larger Hilbert space. Since the noisy gates are modeled as superoperators, we require simulations on the doubled Hilbert space of dimension 104,976. 

The system is initialized in the globally entangled checkerboard state described in Sec.5, $\frac{1}{\sqrt{2}} \left( \ket{1010,g,1} + \ket{0101,g,0} \right)$. For any given layer of our circuit, the operations are applied in the following order:

\textbf{Beam-splitter gates}
\begin{itemize}
    \item Apply beam-splitter propagators for a random durations ($\theta$) and phases $\varphi$ (if BSRP) on modes (1,2) and (3,4), alternating with (1,2) to implement the brick-layer geometry with open boundary conditions.
    \item For every application of a beam-splitter gate, apply a decay propagator to the ancilla qubit as well as spectator modes to account for the idle time while performing the operation.
\end{itemize}


\textbf{On-site Hubbard gates}

For every cavity mode:
\begin{itemize}
    \item Apply a 500 ns decay propagator to all modes and the ancilla to simulate idle time during the swap of the target mode to the buffer cavity.
    \item Apply the SNAP gate propagator to the target mode and transmon. 
    \item Apply decay propagators to all spectator modes and the ancilla to simulate idle time during the SNAP operation.
    \item Reset the transmon ancilla to the ground state and apply a decay propagator to all modes and the ancilla qubit to simulate idle time during the reset operation.
    \item Apply another 500 ns decay propagator to all modes and ancilla to simulate idle time for the swap of the target mode back to the storage cavity.
\end{itemize}

\textbf{Parity Measurements}

For each mode selected for measurements:
\begin{itemize}
    \item Apply a 500 ns decay propagator to all modes and the ancilla qubit to simulate idle time during the swap the target mode to the buffer cavity.
    \item Apply the parity measurement propagator to the target mode and transmon.
    \item Apply decay propagators to all spectator modes and the ancilla to account for idle time during the parity measurement operation.
    \item Apply decay propagators to all modes, the transmon, and the ancilla to account for the idle time during readout. 
    \item Collapse the transmon into the ground or excited state according to the appropriate Born probabilities obtained from its reduced density matrix. Apply a 4\% chance to flip the outcome and collapse onto the opposite parity, corresponding to a readout error.
    \item Apply another 500 ns decay propagator too all modes and ancilla to simulate the idle time for the swap of the target mode back to the storage cavity.
    \item Repeat the step above for every mode that has to be measured.
\end{itemize}
Two clarifying remarks: (1) Although we account for the idle time during swaps between storage and buffer cavities, we do not explicitly simulate the buffer mode, as doing so would significantly increase computational cost. (2) Transmon dissipation and decoherence are included in the SNAP gates and parity measurements. During other operations the transmon remains in its ground state and acts as a spectator, so additional decay modeling is unnecessary.

 The circuit begins with $S$ scrambling layers consisting of BSFP or BSRP gates followed by on-site Hubbard gates but no measurements, which prepares a sufficiently random state (see Sec. \ref{sec:stateprep}). This is followed by $M$ monitored layers where random parity measurements are also included at a rate $p$ and on-site Hubbard gates are optional depending on the model.
Including the overhead from swapping between the storage and buffer cavity and transmon reset, the effective experimental duration of a SNAP operation is $T_{\mathrm{SNAP}} \sim 1.32~\mu\mathrm{s}$, while a parity measurement takes $T_{\mathrm{parity}}\sim1.47~\mu\mathrm{s}$. The total experimental wall time is therefore
\begin{equation}
T_{\mathrm{wall}} = L\,\tilde{S}\,T_{\mathrm{SNAP}} + L M p\,T_{\mathrm{parity}} + L\,\frac{\tau_\mathrm{bs}}{2} (S+M),
\end{equation}
where $\tilde{S}=S$ for BSFP and BSRP monitored circuits and $\tilde{S}=S+M$ for monitored circuits with on-site Hubbard gates, and $\tau_\mathbf{bs}=0.25~\mu\mathrm{s}$ is the average duration of a random BSFP gate. For the end-to-end simulations reported here ($L=4$ and $S=M=8$), at $p=1$, this yields an average circuit duration of $\sim104~\mu\mathrm{s}$ (BSFP) and $152~\mu\mathrm{s}$ (BSFP and on-site Hubbard). BSRP circuits share the same time scale as BSFP circuits since randomizing $\varphi$ doesn't cost any additional time.

\begin{figure*}
    \centering
    \includegraphics[width=0.5\linewidth]{si-figs/error_budget_extra_int.pdf}
    \caption{(a) Averaged entropy curves for circuits with BSFP and on-site Hubbard gates with different noise sources active given $T_1=1.5$ ms: beam-splitter (orange) noise, parity measurement noise (red), SNAP gate noise (in scrambling) (purple), $T_1$ decay only (green), all sources of noise active (blue), and the ideal noiseless case (black). (b) Averaged entropy curves for same noise scenarios but with $T_1=10$ ms.}
    \label{fig:Noisy gate analysis int}
\end{figure*}

\begin{figure*}
    \centering
    \includegraphics[width=1.0\linewidth]{si-figs/error_budget_combined.pdf}
    \caption{\textbf{Analysis of individual noisy gates}(Top row) BSFP circuits with no on-site Hubbard gates (a) Averaged entropy curves with different noise sources active given $T_1=1.5$ ms: beam-splitter (orange) noise, parity measurement noise (red), SNAP gate noise (in scrambling) (purple), $T_1$ decay only (green), all sources of noise active (blue), and the ideal noiseless case (black). (b) Averaged entropy curves for same noise scenarios but with $T_1=10$ ms. (c) Bar graph showing the excess residual entropy generated by different noise sources relative to the case where all noise sources are active. (Bottom row) BSRP circuits (d-f) Same plots as top row but for BSRP circuits.}
    \label{fig:Noisy gate analysis BSFP+BSRP}
\end{figure*}

In Fig.\ref{fig:Noisy gate analysis int}, we show the breakdown of the noisy gates on BSFP circuits with on-site Hubbard gates included in the monitored layers, as discussed in the main text. The simulations consider several individual noise channels including: noisy BSFP gates, noisy parity measurements, noisy SNAP gates (during scrambling), $T_1$ decay, and the case where all noise sources are activate. For each noise channel, we report the excess residual entropy relative to the ideal dynamics. Consistent with the analysis discussed in the main text, noisy parity measurements show a small residual entropy of 0.001(9) compared to noisy BSFP gates, noisy SNAP gates, and $T_1$ decay which are the primary sources of residual entropy. For a cavity lifetime of $T_1=1.5$ ms, when individually activated, $T_1$ decay gives a residual entropy of 0.140(8), noisy BSFP gates give 0.099(7), and noisy SNAP gates give 0.090(8). Increasing the lifetime to $T_1=10$ ms suppresses the contribution from the $T_1$ decay, bringing its residual entropy down to 0.031(7), after which the BSFP and SNAP gate imperfections become the dominant error sources. When all noise channels are activated simultaneously, the residual entropy is 0.251(8) for $T_1=1.5$ ms and 0.182(8) for $T_1=10$ ms.

In the top row of Fig.\ref{fig:Noisy gate analysis BSFP+BSRP}, we show the averaged entropy curves (over 1000 random samples) for BSFP circuits with $L=4$ modes and 16 layers (8 scrambling and 8 monitored). BSFP circuits follow similar qualitative lossy behavior to BSFP and on-site Hubbard gate model. The residual from noisy parity measurements is small at 0.001(9) compared to the rest of the noisy gates. The residual entropy from $T_1$ decay decreases from 0.098(8) to 0.022(7) when the lifetime is increased from $T_1=1.5$ ms to 10 ms. These residuals are smaller than those seen in BSFP and on-site Hubbard circuits due to the shorter wall times introduced from fewer SNAP operations. The residual error from noisy BSFP gates is similar to the BSFP with on-site Hubbard gate case at 0.096(7), and noisy SNAP gates (in scrambling) give 0.084(8). When all noise channels are active, the total residual entropy is 0.228(8) for $T_1=1.5$ ms and 0.174(8) for 10 ms. 

In the bottom row of Fig.\ref{fig:Noisy gate analysis BSFP+BSRP}, we show the breakdown of the noisy gates on BSRP circuits where $\varphi$ is randomized in the beam-splitter gates. The residual entropies from $T_1$ decay are consistent with the  BSFP results with no on-site Hubbard gates since the wall times for both circuits are similar. The residuals from the noisy beam-splitter gates, SNAP gates (in scrambling), and parity measurements are slightly higher than their BSFP counterparts at 0.145(9), 0.126(9), and 0.017(9) respectively. When all noise channels are active, the total residual entropy is 0.284(9) at $T_1=1.5$ ms and 0.229(8) at $T_1=10$ ms. 







\bibliography{references}

\end{document}
