Randomized measurement protocols for lattice gauge theories

Jacob Bringewatt
Joint Center for Quantum Information and Computer Science, NIST/University of Maryland, College Park, Maryland 20742, USA
Joint Quantum Institute, NIST/University of Maryland, College Park, Maryland 20742, USA
Jonathan Kunjummen
Joint Center for Quantum Information and Computer Science, NIST/University of Maryland, College Park, Maryland 20742, USA
Joint Quantum Institute, NIST/University of Maryland, College Park, Maryland 20742, USA
Niklas Mueller
InQubator for Quantum Simulation (IQuS), Department of Physics, University of Washington, Seattle, WA 98195, USA.

Randomized measurement protocols, including classical shadows, entanglement tomography, and randomized benchmarking are powerful techniques to estimate observables, perform state tomography, or extract the entanglement properties of quantum states. While unraveling the intricate structure of quantum states is generally difficult and resource-intensive, quantum systems in nature are often tightly constrained by symmetries. This can be leveraged by the
symmetry-conscious randomized measurement schemes we propose, yielding clear advantages over symmetry-blind randomization such as reducing measurement costs,
enabling symmetry-based error mitigation in experiments,  allowing differentiated measurement of (lattice) gauge theory entanglement structure, and, potentially, the verification of topologically ordered states in existing and near-term experiments.

§ INTRODUCTION

Measurement in quantum mechanics reveals very limited information regarding the structure of the underlying quantum state. This has major  practical implications, e.g., for variational near-term quantum-classical algorithms <cit.>, the verification of quantum devices <cit.>, or when detecting entanglement <cit.> in quantum simulation experiments.
Randomized measurement protocols, such as randomized benchmarking <cit.>,
classical shadows <cit.>, and
entanglement tomography <cit.> are valuable techniques for addressing this  problem. They allow one to estimate many observables from a few measurements <cit.> or extract non-linear quantities, such as purities $\sim \text{Tr}(\rho^k)$, $k\ge 2$ and entanglement entropies <cit.>, potentially
without the massive overhead of traditional state tomography <cit.>, see e.g. <cit.> for a recent overview. Many  techniques are feasible on noisy, near-term quantum devices <cit.>.

A key application for quantum computing  and randomized measurement protocols is simulating quantum many-body systems, with digital or analog devices based on atomic, molecular and optical (AMO), and solid-state systems <cit.>. Quantum simulation promises to address long standing questions in condensed matter, high energy physics and nuclear physics. Examples include simulating non-equilibrium evolution and thermalization <cit.>, thermal systems <cit.> and quantum phases <cit.>.  One important frontier is the study of lattice gauge theories (LGTs) <cit.>  with intricate entanglement structures <cit.> and, potentially, emergent topological phases that have applications in topological quantum computation <cit.>.
A fundamental problem is that standard randomized measurement protocols do not account for the symmetry structure of LGTs, e.g., by randomizing unnecessarily over unphysical states.

(a) Illustration of our proposed symmetry-conscious randomized measurement scheme, which preserves the symmetry structure of states, compared to symmetry-blind approaches. (b) Applications and benefits include cost reduction when measuring entanglement and finding classical representations of quantum states, allowing a rudimentary symmetry-based error mitigation scheme, and enabling symmetry-resolved  measurement of entanglement structure.

We develop symmetry-conscious  protocols to overcome this issue.
Randomized measurement protocols are based on changing the basis via unitaries $U$ drawn from an appropriate ensemble $\mathcal{E}$, i.e. $\rho \rightarrow U \rho U^\dagger$, followed by measurement in this basis, and classical post-processing or
quantum variational techniques to determine quantities of interest.[Schemes using collective measurement on many copies of  $\rho$ (e.g. shadow tomography <cit.>) do not fit into this class of randomized measurement protocols.]
The effectiveness and cost of a scheme depends on the choice of $\mathcal{E}$ and the particular quantities one wants to compute.
For instance, estimating expectation values, $\text{Tr}[{O_m \rho}]$, of $m$-local operators $O_m$ can be done with qubit-local random rotations, $\mathcal{O}(\log(m))$ measurements, and efficient classical postprocessing <cit.>.
For non-linear quantities, i.e. those that depend on $\rho^k$ ($k\ge 2$),
one approach is to choose an ensemble $\mathcal{E}$ that forms a so-called unitary $k$-design <cit.>.

A (classically known) symmetry $\hat{S}$ of $\rho=\bigoplus_s \rho_s$  ($[\hat{S},\rho]=0$), where $s$ labels $\hat{S}$-eigensectors, is lost after randomization unless also $[U,\rho]=0$ for all $U\in\mathcal{E}$.
This is compactly summarized in Fig. <ref>,
comparing symmetry-ignorant versus the symmetry-conscious deeply-scrambling random circuits $U$  proposed in this manuscript. Advantages include a significant reduction in measurement complexity and allowing for rudimentary symmetry-based error mitigation <cit.>. Further, symmetry-conscious randomization enables a direct measurement of LGT entanglement structure—for instance, it allows one to extract distillable versus symmetry entanglement, or the symmetry-resolved entanglement spectrum of any state. These tasks are otherwise impossible without re-learning these symmetries through full state tomography.

Our ultimate goal is to provide a useful tool for near-term exploration of LGT entanglement structure, including obtaining efficient classical representations of quantum states. A prime application of our approach is verifying topologically ordered states in quantum experiments;  realizable  with sufficient or universal programmability, e.g. <cit.>.
We focus on $\Ztwo$ LGT in $1+1$ dimensions (coupled to staggered matter) and pure $\Ztwo$ LGT in $2+1$ dimensions as prototypical models.

The manuscript is organized as follows: In Section <ref>, we present the basics of our symmetry-respecting randomized measurement scheme with the simple example of a one-dimensional particle number-conserving spin system.  In Section <ref>, we extend our analysis to gauge theories, with a particular emphasis on measuring the entanglement structure in $(1+1)$d and $(2+1)$d LGT systems. Finally, in Section <ref>, we summarize our results and discuss applications. The manuscript is supplemented by several Appendices where we discuss various details of the employed numerical techniques.

§ SYMMETRY-CONSCIOUS RANDOMIZATION

Randomized measurement schemes allow one to extract properties of a quantum state from  limited experimental repetitions. We focus on schemes that are “symmetry-conscious” and are based on generating (approximate) $k$-designs that preserve the symmetry-block  structure of a state $\rho = \bigoplus_s \rho_s$ by constructing circuits with the same structure, $U= \bigoplus_s U_s$, to randomize the measurement basis. Our approach is partly inspired by <cit.>, using the fact that randomization can be performed through generators with the same symmetry structure as the Hamiltonian of a given system.

We illustrate our approach first for spin systems with particle number symmetry, $\hat{S}_{N}\equiv\sum_j (\sigma_j^z+1)/2$, where $j$ labels a lattice site and $\sigma^a_j$ ($a=x,y,z$) are Pauli matrices acting on the site. A physical state, pure or mixed, is structured with respect to one or more non-zero symmetry sectors.
Symmetry-conscious randomization is achieved with two-qubit gates of the form

\begin{align}\label{eq:2qubitbasicgate}
u\equiv \begin{pmatrix}
e^{i\theta} & &\\
& \begin{bmatrix}
\end{bmatrix} & \\
& & e^{i \phi}
\end{pmatrix}\,,
\end{align}

where the rows and columns label (from top down) the $\da\da$ (0-particle), $\ua\da$, $\da\ua$ (1-particle) and $\ua\ua$ (2-particle) sectors; $[u_1]$ indicates the $2\times 2$ ($ZYZ$-decomposition) matrix structure,

\begin{align}\label{eq:1qubbitrotation}
\begin{pmatrix}
e^{i(\alpha+\gamma)}\cos(\beta) & e^{-i(\alpha-\gamma)}\sin(\beta) \\
- e^{i(\alpha-\gamma)}\sin(\beta) & e^{-i(\alpha+\gamma)}\cos(\beta)
\end{pmatrix}\,.
\end{align}

If $\alpha$, $\beta$, and $\gamma$ are selected such that $[u_1]$ is drawn from a circular unitary ensemble (CUE) <cit.>,  $\theta$ and $\phi$ are drawn evenly from $[0,2\pi)$, then eq:2qubitbasicgate is a  block-structured unitary acting Haar-randomly on the blocks of equal particle number.[An overall phase of the $2\times2$ CUE matrix was re-expressed in the 0-, and 2-particle sectors by partially absorbing it into a global phase.]  It has the following circuit realization:

\begin{equation}\label{eq:circuit}
\begin{adjustbox}{width=0.48\textwidth}
\begin{quantikz}[transparent, row sep=0.12cm]
\lstick{$j_1$}\qw & \gate{e^{\;i\frac{\alpha}{2} \sigma^z} \;} \qw & \gate[2]{\mathcal{U}(\beta)}\qw & \gate{e^{\;i\frac{\gamma}{2} \sigma^z}\;}  &\ctrl{1} \qw&\gate{\sigma^x} &\ctrl{1} \qw&\gate{\sigma^x} \\
\lstick{$j_2$} \qw & \gate{e^{-i\frac{\alpha}{2} \sigma^z}} \qw & \qw & \gate{e^{-i\frac{\gamma}{2} \sigma^z}} &\gate{P(\phi)}&\gate{\sigma^x}& \gate{P(\theta)}&\gate{\sigma^x}
\end{quantikz}
\end{adjustbox}
\end{equation}

where $\mathcal{U}_{j_1,j_2}(\beta)\equiv \exp\{i\frac{\beta}{2}( \sigma^y_{j_1} \sigma^x_{j_2} - \sigma^x_{j_1} \sigma^y_{j_2})\}$ and $P(x)= \text{diag}(1,\exp\{ix\})$ is the phase gate.

Random circuit scheme for systems with particle number symmetry. Blue squares connected by black lines represent the unitary $u$ in eq:circuit.

Using this basic circuit as a building block, we can construct an ensemble $\mathcal{E}$ of $N_{\mathcal{E}}$ $n$-qubit symmetry-conscious random circuits, as shown in fig:particle_number_circuit, consisting of $\ell$ layers where $n(n-1)/2$ pairs of qubits ($(n-1)(n-2)/2$ if $n$ is odd) are randomly connected by eq:2qubitbasicgate.

For sufficiently large $\ell$ and $N_{\mathcal{E}}$,
this realizes an approximate unitary $k$-design in every symmetry block $s$. We verify this numerically for $k=2$ by considering the following moments, separately for every $s$,

\begin{align}\label{eq:2designquant}
&(\mathcal{B}^s)^{ i'j'k'l'}_{ijkl}
\equiv \langle U^s_{ij} U^{s*}_{i'j'}U^s_{kl} U^{s*}_{k'l'}\rangle
\nonumber\\
& - \frac{d_s^2}{d_s^2-1} \Big[ (\mathcal{A}^s)^{i'j'}_{ij}(\mathcal{A}^s)^{k'l'}_{kl}+ (\mathcal{A}^s)^{k'l'}_{ij}(\mathcal{A}^s)^{i'j'}_{kl}\Big]\,,
\end{align}

where $\langle \dots
\rangle =(1/N_{\mathcal{E}}) \sum_{U\in \mathcal{E}}\dots$ is the $\mathcal{E}$-ensemble average, $d_s$ is the dimension of sector $s$, and $(\mathcal{A}^s)_{ij}^{kl}\equiv \langle U^s_{ij}U^{s*}_{kl} \rangle$. For a 2-design it holds that

\begin{align}\label{eq:def2design}
& (\mathcal{B}^s)^{i'j'k'l'}_{ijkl} \stackrel{\text{2\text{-des.}}}{=}
-\frac{ \delta_{ii'}\delta_{kk'}\delta_{jl'}\delta_{lj'} + \delta_{ik'}\delta_{ki'}\delta_{jj'}\delta_{ll'} }{d_s(d_s^2-1)}\,,
\end{align}

while for a 1-design $(\mathcal{A}^s)_{ij}^{kl}=\delta_{ij}\delta_{kl}/d_s$ <cit.>.
We numerically simulate an ensemble of $N_{\mathcal{E}}$ circuits and compute their deviation from a 2-design $\epsilon$ as the absolute difference between eq:def2design and eq:2designquant; we average over indices and multiply by $d_s(d_s^2-1)$ to make $\epsilon$ dimension-independent. For details see Appendix <ref>.

(a)  Particle number-respecting random circuits $N_{\mathcal{E}}$ required to approximate a 2-design with precision better than $\epsilon =10^{-2}$ in every symmetry sector ($\ell =128$). Inset: Relative sampling cost reduction compared with a symmetry-ignorant scheme. (b) Measurements $N_M$  required to estimate sector-wise $2$-purities with precision better than $\epsilon=0.05$, for fixed $N_{\mathcal{E}}=1,428$ and $\ell=128$. Here $n$ qubits refer to a subsystem, i.e. $n=N_A$, where $N_A=N/2$ is a bipartition of the lattice.

fig:scalingpn(a) displays the size $N_{\mathcal{E}}$ of a circuit ensemble (with $\ell=128$ layers) to approximate a 2-design better than $\epsilon=10^{-2}$ for all particle number blocks of $n=4,6,8,10$ qubit systems. The sampling complexity scales  with the block-dimension $d_s$ as

\begin{align}\label{eq:scaling2design}
N_\mathcal{E}\sim d_s^{\xi} \,,\quad \xi=2.1\pm 0.2\,,
\end{align}

consistent with $\xi=2$. Error bars are found by comparing with the $\ell\rightarrow \infty$ limit (obtained by directly sampling Haar random unitaries in each symmetry block); the fit error for $\xi$ is determined by varying the fit regime, leaving out the largest few blocks, and by varying the required $\epsilon$ by one order of magnitude.
The inset of fig:scalingpn(a) shows the relative reduction $r_s$ in sampling cost compared to using a symmetry-ignorant scheme, where $r_s:=(d_s/d_{\mathcal{H}})^{\xi}$. Largest gains are found away from half filling, which can be understood by comparing $d_{\mathcal{H}}=2^n$ to $d_s = \binom{n}{s}=n!/(s!(n-s)!)$ and using Stirling's approximation to find that $r_s \sim (d_s/d_{\mathcal{H}})^\xi \approx [ 2 (\frac{s}{n})^{\frac{s}{n}}(1-\frac{s}{n})^{1-\frac{s}{n}}  ]^{-n\xi}/ (2\pi n \frac{s}{n} (1-\frac{s}{n}))^{{\xi}/{2}}$, an exponential (in $n$) cost reduction for $s \ll n/2 $ ($\gg n/2$).

We could not numerically determine the scaling with qubits $n$ of the circuit layer depth $\ell$ necessary to reach this precision; however, in Appendix <ref>  we find, for coarser $\epsilon=5\cdot 10^{-2}$, that $l \propto n$ suggesting a $O(n^2)$ circuit complexity to realize a $2$-design. It is worth noting that our numerical approach can be supplemented by rigorously deriving upper bounds for $N_{{\mathcal{E}}}$ in the $\ell\rightarrow\infty$ limit (and the scaling with $\ell$) using techniques similar to those in Ref. <cit.>. We will return to this in future work.

Next, we explore measuring $k$-purities using a prototypical model with particle number symmetry; we consider the following spin Hamiltonian in $(1+1)$d,

\begin{align}\label{eq:HamiltonianModel1}
H = \frac{1}{2a}\sum_{j=0}^{N-1}( \sigma^+_j\sigma^-_{\scriptsize j+1} + {\rm h.c.}) + m\sum_{j=0}^{N-1}(-1)^j \sigma^+_j \sigma^-_j\,,
\end{align}

on a lattice with $N$ sites labelled by $j\in[0,N-1]$ and with periodic boundary conditions (PBCs);
$\sigma_j^\pm \equiv\frac{1}{2} (\sigma_j^x\pm i\sigma_j^y)$, $a$ is the lattice spacing and $m$ a mass term. eq:HamiltonianModel1 conserves particle number, $[H,\sum_j \frac{1}{2} (\sigma^z_j+1)]=0$. As inputs for the particle number symmetry-conscious scheme in fig:particle_number_circuit, we numerically determine ground states of $H$ via exact diagonalization <cit.>.  Because these are simple, unentangled computational basis states in the limit $1/2a \to 0$, we work in the opposite limit with $m\cdot a  = 0.05$ (the ground states are half-filled).  We focus on a subsystem A with $ N_A = N/2$ sites. Symmetries of $\rho_A$ are particle number $s \equiv n_A\in[0,N_A]$ in the subsystem, noting that if particle number is fixed globally, $\rho_A$ is block-diagonal, $[\rho_A, \sum_{j=0}^{N_A-1} \frac{1}{2}(\sigma^z_j+1)]=0$, i.e.
$\rho_A =\bigoplus_s \rho_{A,s} $.

We extract $k$-purities, $\text{Tr}(\rho_{A,s}^k)$ by measuring  the probabilities $P_U(b,s)$ of bitstring $b$ (and symmetry sector $s$) with $N_M$ shots in the basis defined by the random unitary $U$; $k$-purities for $k\geq 2$ are directly related to the $k$-Rényi entropies $S_{A,s}^{(k)}\equiv [1/(1-k)]\log\text{Tr}( \rho_{A,s}^k)$. We follow the approach taken in Refs. <cit.>,
stochastic moments $\langle P_U(b,s)^k\rangle \equiv (1/N_\mathcal{E}) \sum_{U\in \mathcal{E}} P_U(b,s)^k $ are related to $k$-purities via <cit.>,

\begin{align}\label{eq:kpurities}
\langle P_U(b,s)^k\rangle = \frac{1}{D_k} \sum_{\{ a_i \}_k \in \mathbb{N}_0}C_{\{ a_i \}_k} \prod_{j=0}^k \text{Tr}\left[ \rho_{A,s}^j \right]^{a_j}
\end{align}

where $\{ a_i \}_k\equiv a_1,\dots,a_k\in \mathbb{N}_0$ with $\sum_{j=1}^kja_j=k$, $D_k\equiv\prod_{j=0}^{k-1}(d_s+j)$ and $C_{\{ a_i \}_k}\equiv k!/\prod_{j=1}^k( j^{a_j} a_j!  )$; $\langle P_U(b,s)^k\rangle$ and $\text{Tr}\left[ \rho_{A,s}^k \right]$ refer to the $k$-moments of the probabilities and $k$-purity per symmetry sector $s$, respectively, with $ \sum_{b\in s} P_U(s,b) = \text{Tr}[\rho_{A,s}]\equiv p_s \le 1$ and $\sum_s p_s=1$. We assume an ideal quantum machine, and the total measurement cost is $N_{M}\cdot N_{\mathcal{E}}$.

In fig:scalingpn(b), we show the required shot number $N_M$ to measure $S_{A,s}^{(2)}$ with precision better than $\epsilon=0.05$, for fixed $N_{\mathcal{E}}=$ 1,428 and $\ell = 128$ layers,  comparing different system sizes and symmetry sectors. The fit reveals an approximately linear dependence on $d_s$, error bars indicate standard error of the mean for seven independent trials. (In fig:True_purities_particle_number of Appendix <ref> we investigate the $N_{\mathcal{E}}$ dependence in the infinite shot limit.) Together, our results show that the cost of the symmetry-conscious approach, $N_{\mathcal{E}}\cdot N_{M}$, is proportional to $d_s$ instead of  $d_{A}$  (the size of the subsystem). In many cases $d_s\ll d_{A}$, yielding a significant advantage. We shall see that this advantage is exponential for lattice gauge theories.

§ LATTICE GAUGE THEORIES

Next, we turn towards our main objective: lattice gauge theories (LGTs). Such systems have an extensive number of local constraints in the form of Gauss laws, so that the dimension $d_s$ of a physical sector is exponentially smaller than the Hilbert space size. In this setting, the benefits of symmetry-conscious randomization are obvious but go beyond just improving the sample complexity. We will demonstrate that such schemes also allow one to compute entanglement structure and to verify topologically ordered states without having to reconstruct the symmetries via full state tomography.

(a) Illustration of $\Ztwo$ LGT in $1+1$ dimensions,  depicted is an “even” site, $(-1)^j=1$, where spin up (down) is the presence (absence) of a $\Ztwo$ charge, $\mathcal{Q}_j=+1$ $(0)$. (b) Single layer of a random-measurement circuit with non-local gates $2+|j_1-j_2|$ qubit gates, an extension of the strategy in section <ref>. (c) Near-term strategy based on  3-qubit unitaries. (d) Illustration of the symmetry structure of $\rho_A$.

§.§ Z2 LGT in (1+1)d

We first consider $\Ztwo$ LGT coupled to staggered matter in $1+1$ dimensions, with Hamiltonian,

\begin{align}\label{eq:Z2LGT1d}
H = &\frac{1}{2a} \sum_{j=0}^{N-1}( \sigma^+_j \tilde{\sigma}^x_{j,j+1} \sigma^-_{j+1} + {\rm h.c.} )
\nonumber\\
&+ m\sum_{j=0}^{N-1}\frac{(-1)^j}{2}(1+\sigma^z_j)   + e\sum_{j=0}^{N-1} \tilde{\sigma}^z_{j,j+1},
\end{align}

Gauss laws

\begin{equation}\label{eq:Gausslaw1d}
\hat{G}_j \equiv e^{i\pi\mathcal{Q}_j} \tilde{\sigma}^z_{j-1,j} \tilde{\sigma}^z_{j,j+1},
\end{equation}

and periodic boundary conditions (PBC),
where $j$ labels a site and $[H,\hat{G}_j]=0$; $\sigma^{b}_j$ ($\tilde{\sigma}^b_{j,j+1}$), $b=x,z$, are Pauli operators residing on the sites (links) of the lattice and representing the matter (gauge) degrees of freedom of the theory; $m$ is a mass parameter, $e$ the $\Ztwo$ coupling, $a$ the lattice spacing, $\sigma^\pm_j \equiv (\sigma^x_j \pm i\sigma^y_j)/2$, and $\mathcal{Q}_j\equiv (\sigma^z_j+(-1)^j)/2$ is the $\Ztwo$ charge. This is compactly summarized in fig:overviewLGT1(a) where Gauss law eigensector with $\hat{G}_j | \psi \rangle = +1| \psi\rangle$ are physical.

(a) Sector-wise $k$-purities of the $\mathbb{Z}^{1+1}_2$ ground state measured based on eq:kpurities as a function of $N_\mathcal{E}$, for $N_A=N/2=5$, $m\cdot a=0.1$, $e/m=8$ and $\ell=32$. (b) Von Neumann Entropy from $k$-purities/Rényi entropies using a $4^{\rm th}$ order finite-difference approximation of eq:vonNuenanDef. (c) Symmetry- (symm.), distillable  (dist.) and total entanglement entropies as a function of $e/m$ for  $N_\mathcal{E}$=2000, $m\cdot a=0.1$ and $\ell=32$. Symbols represent random measurement results, dotted lines represent approximating $S_E$ by (exact) Rényi-entropies up to $k=4$, solid lines are exact results.

fig:overviewLGT1(b) and (c) depict random measurement strategies for this model, we investigate a $N_A$ site subsystem with boundary conditions ending in a matter site. The first strategy, (b), is an extension of the circuit in Section <ref>, eq:circuit, made gauge invariant by introducing Wilson lines connecting sites $j_1$ and $j_2$,

\begin{align}
W_{j_1,j_2} \equiv \prod_{j=j_1}^{j_2-1}\tilde{\sigma}^x_{j,j+1}
\end{align}

and with circuit representation,

\begin{equation}\label{eq:circuit_gauge}
\begin{adjustbox}{width=0.49\textwidth}
\begin{quantikz}[transparent, row sep=0.22cm]
\lstick{$j_1$}\qw & \gate{e^{\;i\frac{\alpha}{2} \sigma^z_{j_1}} \;} \qw & \gate[3]{\mathcal{W}_{j_1,j_2}(\beta)}\qw & \gate{e^{\;i\frac{\gamma}{2} \sigma^z_{j_1}}\;}  &\ctrl{2} \qw & \gate{\sigma^x} &\ctrl{2} & \gate{\sigma^x} \\
\qwbundle[alternate]{} & \qwbundle[alternate]{} & \qwbundle[alternate]{} &  \qwbundle[alternate]{}&\qwbundle[alternate]{}& \qwbundle[alternate]{} & \qwbundle[alternate]{}& \qwbundle[alternate]{}
\\
\lstick{$j_2$} \qw & \gate{e^{-i\frac{\alpha}{2} \sigma^z_{j_2}}} \qw & \qw & \gate{e^{-i\frac{\gamma}{2} \sigma^z_{j_2}}} &\gate{P(\phi)} & \gate{\sigma^x} & \gate{P(\theta)} &  \gate{\sigma^x}
\end{quantikz}\,,
\end{adjustbox}
\end{equation}

where $\mathcal{W}_{j_1,j_2}(\beta)\equiv \exp\{i\frac{\beta}{2} (\sigma_{j_1}^y  W_{j_1,j_2}\sigma^x_{j_2}-\sigma_{j_1}^x W_{j_1,j_2} \sigma^y_{j_2} )\}$; the angles $\alpha$, $\beta$, $\gamma$, $\phi$ and $\theta$ are randomly drawn as before.
The abbreviated middle qubit bundle in eq:circuit_gauge refers to qubits representing gauge links (orange in fig:overviewLGT1),  $j_1$ and $j_2$ start and end on matter sites (blue in fig:overviewLGT1). This results in a  $2+|j_1-j_2|$ qubit unitary which is not a feasible strategy near-term. Because of this we will focus on a 3-qubit unitary strategy, fig:overviewLGT1(c), at the cost of somewhat deeper circuits to obtain $k$-designs. We demonstrate numerically that these circuit form approximate 2-designs in Appendix <ref>, and focus here on measuring entanglement.
We note that (b) and (c) generate $k$-designs in every symmetry sector, despite the fact that no explicit 2-qubit entangling operation is performed between gauge sites (orange). This is a consequence of gauge symmetry: For this simple model the gauge link degrees of freedom are not truly independent because
they could have been eliminated using Gauss' law [A caveat is that with PBCs a gauge zero mode cannot be integrated out. A $k$-design acting on the full system, not just a subsystem,  would include a modification which depends on how the remaining bosonic mode is digitized; see e.g. <cit.> for bosonic random measurement schemes.].

Focusing on the entanglement properties of ground states, we consider a bipartition $\rho_A\equiv \text{Tr}_{\bar{A}}(\rho)$ of the lattice with $N_A\equiv N/2$ sites. We first discuss the symmetry structure of $\rho_A$.
As before, $\rho_A$ is block-diagonal in particle number $n_A$ (i.e., $[\rho_A,\sum_{j=0}^{N_A-1} (1+\sigma^z_j)/2]=0$) but additionally has  symmetries beyond those of the non-gauge spin model. They are illustrated in fig:overviewLGT1(d): Because of Gauss' law, on the right (left) boundary $j=j_R\equiv N_A$ ($j=j_L\equiv 0$), out-going (in-going) electric fields can be written as

\begin{align}
\tilde{\sigma}^z_{j,j+1}=\exp\{ i\pi \mathcal{Q}_{j}\}   \tilde{\sigma}^z_{j-1,j }\,.
\end{align}

These operators are symmetries of $\rho_A$ (i.e. $[ \tilde{\sigma}^z_{j_R,j_R+1} , \rho_A]=[ \tilde{\sigma}^z_{j_L-1,j_L} , \rho_A]=0$)  if $\rho$ is physical (i.e. Gauss' law respecting) and result in a four-block symmetry structure. Together, we label the sectors as $s\equiv n_A^{s_L,s_R}$ where $s_{L/R}=\ua/\da$ and $n_A\in [0,N_A]$ is particle number.

In fig:LGT1dInnsbruckScheme(a), we show $k$-purities, $\text{Tr}[\rho_{A,s}^k]$  reconstructed by inverting eq:kpurities, measured for $N=10$ ($N_A=5$), $m\cdot a=0.1$, $e/m$=8 and $\ell=32$. We plot them as a function of  $N_{\mathcal{E}}$; measurements are obtained in the infinite shot  limit ($N_M\rightarrow\infty$).  The $k=1$ result $p_s\equiv \text{Tr}[\rho_{A,s}]$  is recovered exactly by design, higher ($k\ge 2$) purities are recovered for sufficiently large $N_{\mathcal{E}}$ (data is shown up to $N_{\mathcal{E}}=2^{14}$). Not shown is the cost in $N_M$ to obtain constant error which, as in Sec. <ref>,  scales with $d_s$.  In fig:LGT1dInnsbruckScheme(b), we  estimate the von-Neumann entropy per symmetry block,

\begin{align}\label{eq:vonNuenanDef}
S_s\equiv \text{Tr}_s[\bar{\rho}_{A,s}\log(\bar{\rho}_{A,s})] = - \lim_{k\rightarrow 1^+} \frac{{\rm d}}{{\rm d} k} \text{Tr}_s[\bar{\rho}_{A,s}^k],
\end{align}

where $\bar{\rho}_{A,s}\equiv {\rho}_{A,s}/p_s$. We make use of a $4^{\rm th}$ order finite-difference approximation of the derivative to derive eq:vonNuenanDef from the measured $k$-purities;
error bars are obtained from comparing a $4^{\rm th}$ and $3^{\rm rd}$ order derivative.
Finally, fig:LGT1dInnsbruckScheme(b) shows the decomposition of the von-Neumann entropy, $S=-\text{Tr}_A[\rho_A \log(\rho_A)]$ into a symmetry- (`classical' entanglement) and distillable component, $S^{  {\rm symm.}}+S^{\rm dist.}$,

\begin{align}\label{eq:entcomponents}
%  S=S^{  {\rm symm.}}+S^{\rm dist.}\equiv -\sum_s \{ p_s \log(p_s)+ p_s S_s \}\,,
S^{  {\rm symm.}} \equiv -\sum_s  p_s \log(p_s)\,, \quad
S^{\rm dist.}\equiv \sum_s p_s S_s \,,
\end{align}

as a function of $e/m$ and for $m\cdot a=0.1$.

(a) Bottom: Shadow-reconstructed  symmetry-resolved Schmidt spectrum, $m\cdot a=0.05$, $e/m=1$, $\ell=32$ and $N_{ \mathcal{S}}=2^{17}$. Inset: Full spectrum. Top: Probability per sector $p_s$. (b) Sector-wise relative entropy $S(\bar{\rho}_{A,s}||\bar{\sigma}_{A,s})$ between  (normalized) exact $\bar{\rho}_{A,s}$ and shadow-reconstructed reduced density matrices $\bar{\sigma}_{A,s}$, as a function of shadow number $N_{\mathcal{S}}$ with $m\cdot a=0.05$, $e/m=1$ and $\ell=32$.

Next, we  explore classical representations of $\rho_A$, starting with the classical shadow formalism of <cit.>. The basic idea is to randomize $\rho_A \rightarrow U \rho_A U^{\dagger}$ and perform a computational basis measurement yielding a bitstring $b$ from which the sector $s$ can be read off. A symmetry-conscious shadow is $U^\dagger_s | b,s \rangle \langle b,s | U_s$, where the subscript $s$ indicates that one works in block $s$ with dimension $d_s$. The ensemble of random rotations yields a CUE-random quantum channel
$\mathcal{M}[\rho_{A}]=\bigoplus_s\mathcal{M}_s[\rho_{A,s}]$, and with many measurements and, consequently, many shadows one obtains classical sector-wise state representations by taking the expectation value

\begin{align}
\bar{\sigma}_{A,s} \equiv \mathbb{E}\left[ \mathcal{M}_s^{-1} \big( U^\dagger_s | b,s \rangle \langle b,s | U_s\big)\right]\,,
\end{align}

where $\mathbb{E}[\dots]\equiv (1/N_{\mathcal{S}})\sum_{i=0}^{N_{\mathcal{S}}-1}[\dots] $ is the  $N_\mathcal{S}$ shadow average and $\mathcal{M}_s^{-1}(X)\equiv (d_s+1)X - \text{Tr}_s[X] \, \mathbb{I}_s$; $\text{Tr}_s$ and $\mathbb{I}_s$ are the sector-wise trace and identity, respectively. The bar indicates normalization, i.e. $\bar{\sigma}_{A,s}\equiv {\sigma}_{A,s}/p_s$,
$\text{Tr}_s[\bar{\sigma}_{A,s}]=1$, $p_s\equiv \text{Tr}_s[{\sigma}_{A,s}]$ is simply the number of shadows measured in one sector relative to $N_{\mathcal{S}}$.

fig:shadows1d compactly summarizes the results of this analysis, showing the symmetry-resolved Schmidt spectrum of $\rho_A$ in (a), comparing exact results (empty squares) versus shadows (filled triangles) for $m\cdot a=0.05$, $e/m=1$, $\ell=32$ and $N_{ \mathcal{S}}=2^{17}$. The eigenvalue spectrum is well reproduced down to $P_{s,\lambda} \approx 10^{-2}-10^{-3}$; the inset shows the full spectrum and the top of (a) the accurately recovered probability for each block, $p_{s}$.
fig:shadows1d(b) shows the relative entropy, $S(\bar{\rho}_{A,s}|| \bar{\sigma}_{A,s})\equiv  \text{Tr}_s[\bar{\rho}_{A,s}\log(\bar{\rho}_{A,s})-\bar{\rho}_{A,s}\log(\bar{\sigma}_{A,s})]$, between exact $\bar{\rho}_{A,s}$ and  shadow-reconstructed $\bar{\sigma}_{A,s}$, as a function of shadows per sector $p_s\cdot N_{\mathcal{S}}$. The fit (black curve) indicates approximate power law $\sim (p_s N_{\mathcal{E}})^{-\frac{1}{2}}$ convergence. We note that automatically fewer shadows are  sampled in less important sectors.

While $\mathbb{Z}_2^{1+1}$ LGT serves as a useful case study, we next consider $2+1$ dimensions where our ability to extract LGT entanglement structure provides a window to studying topologically ordered systems.

§.§ Detecting topological order from random measurements: Z2 LGT in (2+1)d

(a) Illustration of  $\Ztwo^{2+1}$ LGT, including Hilbert space and Gauss law constraints $G_j$. (b) Random measurement circuits: We work with an even-odd alternating layers consisting of randomly placed `electric' rotations $R_z(\alpha)$, $R_z(\gamma)$ and plaquette rotations $U_\square(\beta)$ (orange squares), approximating $k$- designs for sufficient circuit depth.  (c) The symmetry structure of $\rho_A$ originates from Gauss laws at entanglement boundaries and a non-local `ribbon' operator $V_x^A$ spanning the two entanglement cuts.

Finally, we consider $\Ztwo$ LGT in (2+1)d spacetime dimensions ($\Ztwo^{2+1}$),  graphically illustrated in Fig. <ref>(a), consisting of spin 1/2 degrees of freedom placed on the links $(j,b)$ of a two-dimensional rectangular $N_x\times N_y$  lattice, where $j=(j_x,j_y)$ and $b=x,y$ is the direction of a link. The Hamiltonian is given by

\begin{align}\label{eq:Z22p1Hamilt}
H= -K \sum_{j}\sigma^x_{j,x}\sigma^x_{j+\hat{x},y}\sigma^x_{j+\hat{y},x}\sigma^x_{j,y} -g \sum_{j,b} \sigma^z_{j,b}\,,
\end{align}

and superselection sectors are determined by Gauss laws ($[H, \hat {G}_j]=0$),

\begin{align}\label{eq:glaw_z2_2d}
\hat G_j\equiv \prod_{l\in +(j)} \sigma^z_l,
\end{align}

where the product is over neighboring links to each lattice site $j$, see fig:2d_setup(a) for illustration where $l$ runs over links labelled (1,2,3,4); physical states obey $\hat{G}_j\ket{\psi}=+1\ket{\psi}$. We assume periodic boundary conditions along the $x$-direction, and fixed or periodic boundary conditions (BCs) in the $y$-direction; for fixed BCs, eq:glaw_z2_2d involves three links at the $y$-boundary.

We focus on entanglement properties of $\rho_A\equiv\mathrm{Tr}_{\bar{A}}(\rho)$, where $A$ is a bipartition obtained by separating the lattice along the $x$-direction, with entanglement cuts at $j_x=0$ and $j_x=N_x^A-1$ so that $N_x=N_x^A+N_x^{\bar{A}}$. The boundary is such that system $A$ contains the  $y$-direction links at $j_x=0$ and $j_x=N_x^A-1$, see fig:2d_setup(c) or <cit.>
for more details, with
$N_\square=(N_x^A-1)\times N_y$ plaquettes ($2N_x^AN_y - N_y - N_x^A$ qubits) in $A$.
We seek a gauge invariant, symmetry-respecting family of random circuits that form $k$-designs acting on every symmetry sector of $\rho_A$.
Our ansatz is illustrated in fig:2d_setup(b), made of alternating even-odd half-layers, consisting first of $R_z(\alpha)\equiv \exp\{ i \alpha \sigma^z_i\}$
rotations randomly placed at one side $ i \in \{ a,b,c,d\}$ of a plaquette (orange squares), followed by $U_\square \equiv U_{\square}(\beta) \equiv\{ i \beta \sigma^x_a \sigma^x_b \sigma^x_c\sigma^x_d \} $, and again by a  $R_z(\gamma)\equiv \exp\{ i \gamma \sigma^z_i\}$, placed at the same random $ i \in \{ a,b,c,d\}$. For every plaquette, the angles $\alpha,\beta,\gamma$ are drawn according to a $ZXZ$ decomposition of a CUE matrix,

\begin{equation}
U_{\rm CUE}\equiv e^{i\delta}\begin{pmatrix}
e^{i(\alpha+\gamma)}\cos\beta & ie^{-i(\alpha-\gamma)}\sin\beta\\
ie^{i(\alpha-\gamma)}\sin\beta & e^{-i(\alpha+\gamma)}\cos\beta\\
\end{pmatrix},
\end{equation}

where the phase $\delta$ is irrelevant. In Appendix <ref>, we verify numerically that these random circuits form an approximate unitary 2-design ($k$-design).

Before continuing, we discuss the symmetries of $\rho_A\equiv \text{Tr}_{\bar{A}}[\rho]$, depicted in Fig. <ref>(c).   The Gauss laws at entanglement boundaries allow us to write

\begin{align}\label{eq:symmop2d}
\sigma^z_4 = \sigma^z_1 \sigma^z_2\sigma^z_3\,.
\end{align}

Here, $\sigma^z_4$ is the `electric field' operator just outside A, and $\sigma^z_1 \sigma^z_2\sigma^z_3$ is just inside A. For all $2N_y$ boundary sites, the operator $\sigma^z_1 \sigma^z_2\sigma^z_3$ is a symmetry of $\rho_A$, i.e. $[\sigma^z_1 \sigma^z_2\sigma^z_3,\rho_A]=0$, if $\rho$ is physical (i.e. Gauss law respecting). We label simultaneous eigensectors of all eq:symmop2d as $s\in\{\ua,\da\}^{2^{2N_y}}$; an example  is $s=\ua\ua\,  \da\ua$ (for $N_y=2$) where the first $N_y$ bits are the eigenvalues of eq:symmop2d at $j_x=0$ and the other $N_y$ at $j_x=N_x^A-1$. Additionally, a `ribbon' operator

\begin{align}
V^A_x\equiv \prod_{l \in \mathcal{C}} \sigma^z_l
\end{align}

commutes with $\rho_A$, $[\rho_A,V^A_x]=0$, where
$l\in \mathcal{C}$ indicates the links intersected by a contour $\mathcal{C}$ through the centers of plaquettes, from one boundary to the other in an arbitrary path, see the green dashed line in fig:2d_setup(c). For fixed $y$-BCs $V^A_x$ is not  independent  (it is determined by fixing all sectors eq:symmop2d), but for $y$-PBC its eigensectors are independent and labelled by an additional $\ua/\da$, so that $s\in \{ \ua,\da\}^{2^{2N_y+1}}$.

We focus again on entanglement, display in Fig. <ref>(a)  $k$-purities $\mathrm{Tr}[\rho_{A,s}^k]$ ($k=1,2,3,4$) of the $\mathbb{Z}_2^{2+1}$ ground state at $\epsilon\equiv g/K=0.1$, following the approach outlined in section <ref>, as a function of $N_{\mathcal{E}}$ and in the
infinite shot limit  ($N_M\rightarrow\infty$), with fixed $y$-BCs,  $N_x\times N_y = (3+5)\times 2$ and $\ell=64 $. As before, the $k=1$ results, $p_s =\text{Tr}[\rho_{A,s}]$, are exact by design while $k=2,3,4$-entropies are reproduced with increasing  $N_{\mathcal{E}}$. In fig:LGT2dEHT(b), we also show the Schmidt spectrum $P_{s,\lambda}$ for $\epsilon=0.2$, comparing shadow-reconstructed ($N_{\mathcal{S}}=2^{16}=65,536$, triangles) versus exact results (squares). Different colors represent symmetry sectors $s$ with weights $p_{s} \equiv \text{Tr}_s [\rho_{A,s}] $  in the top panel. Large Schmidt values are well reproduced  down to approximately $10^{-2}$, beyond which we observe significant deviations.

Aiming at higher precision, we explore an alternative approach, Entanglement Hamiltonian Tomography (EHT) <cit.>. The basic idea is to parameterize the reduced density matrix by an Entanglement Hamiltonian (EH),

\begin{align}
H_A=-\log{ [\rho_A]}\,.
\end{align}

A (heuristic)  parameterization of ground state EHs, inspired by the Bisognano-Wichmann (BW) theorem <cit.>, is

\begin{align}
H_A \equiv H_A[\{\beta_{\mathcal{O}} \}] = \sum_{\mathcal{O}} \beta_{\mathcal{O}} H_{\mathcal{O}}\,,
\end{align}

where $ H_{\mathcal{O}}$ are the local operators comprising the physical Hamiltonian, eq:Z22p1Hamilt, and $\beta_{\mathcal{O}}$ a local `temperature', varying with the distance of $\mathcal{O}$ from the entanglement cut(s) (also depending on $j_y$ if translation invariance in $y$ is broken, i.e. for fixed $y$-BCs). The applicability and accuracy of this ansatz was investigated for $\mathbb{Z}_2^{2+1}$ in <cit.>.

(a) Sector-wise k-purities of $\rho_A$ of the $\mathbb{Z}_{2}^{2+1}$ ground state, from eq:kpurities, for $N_x\times N_y = (3+5)\times 2$,  $\epsilon $=0.1 and $\ell$=64 layers and fixed BC in $y$. (b) Bottom: Shadow- ($N_{\mathcal{S}}=2^{16}$) versus BW-EHT-reconstructed ($N_{\mathcal{E}}=50$, $N_{M}=1024$)  Schmidt spectrum $P_{s,\lambda}$; with $N_x\times N_y = (3+5)\times 2$, $\epsilon $=0.2 and $\ell$=64. Top: probability per sector $p_s$; inset: entanglement spectrum $\xi_{s,\lambda}=-\log(P_{s,\lambda})$. (c) Entanglement gap $\Delta_{\xi}$ between  `low- and high-energy' parts of the ES as a function of $\epsilon$, reconstructed using the BW-EHT scheme, for $N_x\times N_y = (3+3)\times 2$ (PBC in $y$), $N_{\rm BW}=50$, $N_{\rm shots}=1024$, and $\ell=64$. A horizontal red line is the infinite volume limit, $\epsilon_c=0.33\pm0.01$. Right: Symmetry-resolved Entanglement spectra $\xi_{s,\lambda}$ for $\epsilon=0.075,03,0.5$.

To extract the $\{\beta_{\mathcal{O}} \}$, we follow <cit.> and first measure probabilities $P_U(b,s)$ in $N_{\mathcal{E}}$ random bases and with $N_M$ shots each. We then minimize the following functional via classical post-processing,

\begin{align}
\chi^2 = \sum_{b,s} \langle \Big(  P_U(b,s) - \text{Tr}_s \big[  \rho_{A,s}U_s | b,s \rangle \langle b,s | U^\dagger_s  \big] \Big)^2 \rangle_{ \mathcal{E}}
\end{align}

where $ \rho_{A,s}\equiv \rho_{A,s}[\{\beta_{\mathcal{O}} \} ] \sim\exp\{ -H_{A,s}[\{\beta_{\mathcal{O}} \} ]\}$, normalized so that $\text{Tr}_s[\rho_{A,s}[\{\beta_{\mathcal{O}} \} ]=p_s$,  $H_{A,s}[\{\beta_{\mathcal{O}} \} ]\equiv \sum_{\mathcal{O}}\beta_{\mathcal{O}} H_{s,\mathcal{O}}$ is the EH with $H_{s,\mathcal{O}}$ restricted to a symmetry sector $s$. Because the BW optimization is performed sector-wise, matrices of size $d_s$ are involved, versus  the dimension of $A$, $d_{A}$. In the infinite measurement limit,  the optimization will yield one universal set $\{\beta_{\mathcal{O}}\}$ for all $s$, but in practice we work with normalized $\bar{\rho}_{A,s}\equiv {\rho}_{A,s}/p_s$ and $\bar{P}_U(b,s) = {P}_U(b,s) / p_s$, so that the extracted $\{\beta_{\mathcal{O},s}\}$ depend on that normalization and differ from $\{\beta_{\mathcal{O}}\}$,  see Appendix <ref> for details.

Results of the BW-EHT-optimization for the Schmidt spectrum $P_{s,\lambda}$ are displayed as diamond symbols in fig:LGT2dEHT(b) along with shadow results, for $N_{\mathcal{E}}=50$  and $N_M=1024$. Despite comparable cost ($N_{\mathcal{E}}\cdot N_M=51,200$) relative to the classical shadow approach ($N_{\mathcal{S}}=2^{16}=65,536$),  BW-EHT  reproduces the eigenvalue spectrum much more accurately; values as small as $10^{-6}$ are approximately recovered and even eigenvalues as small as $10^{-11}$ are not far off. The inset of fig:LGT2dEHT(b) shows the entanglement spectrum (ES), i.e. the spectrum of the EH, which is also well reproduced.
The apparent advantage of the BW-EHT approach comes at the expense of generality, it is tailored for ground states (it can be extended to non-equilibrium states <cit.>) while classical shadows work regardless of the state.

Enabled by the performance of the BW-EHT optimization, we focus on a practical application: detecting topological order (TO) of quantum states. Ground states of $\mathbb{Z}_2^{2+1}$ are separated (in the infinite volume limit) into topologically ordered, $\epsilon < \epsilon_c$, and trivial states, $\epsilon > \epsilon_c$, with a phase transition at a critical coupling $\epsilon_c$. Li and Haldane's entanglement-boundary conjecture  <cit.>,
asserts that TO states are `entanglement-gapped', i.e. their ES has separated low energy (large Schmidt values)
and a high energy (small Schmidt values) parts. Further, the low lying part is (up to rescaling) identical to the  spectrum of a conformal field theory (CFT) describing gapless excitations at the edge of the system. We focus here on measuring the existence of an entanglement gap $\Delta\xi$ to detect TO, which has been shown as very a robust order parameter for the TO transition in this model even for very small systems <cit.>.

Without loss of generality, to reduce finite size effects, we focus on periodic boundary conditions in $y$ (a torus) for the BW-EHT analysis. In fig:LGT2dEHT(c) we show the entanglement gap $\Delta \xi$ for $N_x \times N_y = (3+3)\times 2$,  $N_{\mathcal{E}}=50$, $N_M=1024$, and $\ell=64$ (black diamonds), compared to exact results (black solid line). Error bars represent the combined statistical error due to finite $N_{\mathcal{E}}$ and $N_M$, see Appendix <ref> for details. A vertical red line indicates the infinite volume extrapolated value $\epsilon_c = 0.33\pm 0.01$ <cit.>. We also show result for $N_x \times N_y = (3+3)\times 3$ (green dashed line) and $N_x \times N_y = (3+3)\times 4$ (blue dashed line), taken from <cit.> and approaching the infinite volume limit to within less than 10%. Side panels show the BW-EHT-reconstructed sector-wise ES for $\epsilon=0.075,0.3,0.5$, demonstrating the closing of the entanglement gap (gray shaded area) at $\epsilon_c$ where our results approximately reproduce the phase transition.

We could not numerically simulate systems larger than $N_x \times N_y=(3+3)\times 2$. For example, a $N_x \times N_y=(3+3)\times 4$ lattice of 48 qubits (20 qubits in the subsystem) exhausts our classical computational resources [The results for $N_x \times N_y=(3+3)\times 4$ and $(3+3)\times3$ were obtained using exact diagonalization and working with dual formulations of $\mathbb{Z}_2^{2+1}$ <cit.>.]. However,  the classical (shadow- or BW-EHT-) analysis is simple for such a system if it were prepared in experiment because, while $d_{A}=2^{2N_x^A N_y-N_y}=2^{20}=1,048,576$, the analysis is restricted to symmetry blocks of only $d_s=2^{N_x^A N_y-N_y}=2^{8}=256$ states. This is a significant (in fact exponential) reduction in the space over which the state is randomized, but $d_s$ still grows exponentially with the subsystem size, albeit much slower than $d_A$.

§ CONCLUSION AND OUTLOOK

In this manuscript, we proposed randomized measurement protocols for lattice models that leverage symmetries, focusing primarily on LGT entanglement structure exploration. We devised deep-scrambling circuits that realize symmetry-conscious $k$-designs and illustrated their use in simple gauge and non-gauge model examples. Our approach is intuitive and, therefore, easily generalizable: by examining the physical Hamiltonian one can readily identify basic symmetry-preserving interactions which can be used as the generators of a randomized measurement scheme. Consequently, if a particular physical Hamiltonian can be realized, so can our measurement scheme.

Our goals were twofold. Firstly, we aimed at finding a randomization scheme lower in cost than existing ones by avoiding randomizing over non-relevant Hilbert space parts. Our ansatz has a   sampling cost (to realize a $2$- ($k$-)design) that scales with block size $d_s$, instead of  Hilbert space dimension $d_{\mathcal{H}}$. This reduction can be exponential, e.g., for particle number conserving systems away from half-filling, or for LGTs where it comes from randomizing only over the physical versus unphysical Hilbert space parts.

A second goal was to provide a practical scheme for measuring LGT entanglement structure, a potential useful route e.g.,
for quantum simulating high energy and nuclear physics <cit.>, e.g. to understand Quantum Chromodynamics (QCD) where entanglement is largely unexplored <cit.>, or detecting topologial order.
We illustrated our approach in a simple (1+1)d LGT example, $\mathbb{Z}_2$ coupled to staggered matter, where we extracted symmetry-resolved $k$-purities and von Neumann entanglement entropies, and separated their symmetry- and distillable components.

We then focused on $\mathbb{Z}_2$ in $2+1$ dimensions where the intricate structure of gauge symmetric states can lead to topologically ordered (TO) phases. These are currently receiving great attention, including experimental  realizations in AMO and solid-state platforms <cit.>, motivated by applications such as fractional quantum Hall effect states <cit.> or fault-tolerant quantum computation and storage <cit.>. A difficulty is that TO cannot be probed by measuring local operators, a serious impediment for its experimental verification. Our approach to overcome this is based on measuring the entanglement structure of such systems. While the importance of entanglement as a robust indicator of topological order was realized long ago <cit.>, we developed a concrete random-measurement scheme, following the logic of Li and Haldane <cit.>, that uses a state presentation in terms of Entanglement Hamiltonians (EH) and is based on measuring entanglement gaps of their (symmetry-resolved) spectrum using a tomographic protocol based on the Bisognano-Wichmann theorem <cit.>. Remarkably, performing random measurements on very small  subsystems as small as $N_x^A \times N_y = 3\times 2$, we observe a relatively sharp TO-to-trivial phase transition. While our focus was on $\mathbb{Z}_2^{2+1}$, the protocol
can be easily generalized other LGT(-like) systems.

A benefit of our approach, not explicitly explored in the main text, is that symmetry-conscious randomization allows for a rudimentary, but useful, near-term error-mitigation strategy similar to that discussed in <cit.>. A feature of symmetry-conscious randomization is that symmetries of the input states are not lost and can be measured. Thus, machine errors that violate those symmetries are detectable after randomization, suggesting that e.g., a postselection of measurement results can improve the computation (at the cost of reduced statistics).

There are many future extensions of our work. For example, while our approach significantly reduces algorithmic costs compared to a symmetry-ignorant scheme, extracting entanglement entropies and structure still relies on classical post-processing which ultimately scales exponentially with system size, an issue which can be addressed with quantum variational <cit.> and machine learning techniques <cit.>. In future work, we will also work towards providing analytical performance guarantees for circuit depth and sampling complexity, which we have explored only numerically. We also note that realizing an approximate $k$-design is a sufficient, but not necessary, condition for randomized measurement protocols, and more studies regarding optimal randomization for certain observables are needed. We also emphasize that, while there are encouraging indications <cit.>, the robustness of our scheme against experimental imperfections and noise should be investigated.

Finally, we point out related work <cit.>, following a similar idea for fermionic systems with particle number symmetry, and also demonstrating a significant cost advantage. While not programmable enough to realize $k$-designs, we think this is a very useful approach for fermionic entanglement tomography.
We also point out Ref. <cit.>, proposing randomized measurement schemes that take advantage of a priori knowledge about  observables of interest to improve sampling complexity. and  Ref. <cit.>, which, building on Ref. <cit.>, proposes randomized measurement schemes to extract symmetry-resolved purities. Their scheme is based on local random unitary transformations, in contrast to the symmetry-preserving unitaries we consider here.

§ ACKNOWLEDGMENTS

We thank Alexander F. Shaw for early discussions leading to the formulation of this project.
We also thank
Zohreh Davoudi,
Jonas Helsen, Hsin-Yuan (Robert) Huang,
Martin Savage, Guilia Semeghini, Torsten Zache and Peter Zoller for discussions.
J.B. was funded in part by the Heising-Simons Foundation, the Simons Foundation, and National Science Foundation Grant No. NSF PHY-1748958 and by the U.S. Department of Energy (DOE) ASCR Accelerated Research in Quantum Computing program (award No. DE-SC0020312), DoE QSA, NSF QLCI (award No. OMA-2120757), DoE ASCR Quantum Testbed Pathfinder program (award No. DE-SC0019040), NSF PFCQC program, AFOSR, ARO MURI, AFOSR MURI, and DARPA SAVaNT ADVENT. N.M. acknowledges funding by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, InQubator for Quantum Simulation (IQuS) (<https://iqus.uw.edu>) under Award Number DOE (NP) Award DE-SC0020970; and, during early stages, by the U.S. Department of Energy’s Office of Science, Office of Nuclear Physics under Award no. DE-SC0021143 for quantum simulation of gauge-theory dynamics on near-term quantum hardware. This work was enabled, in part, by the use of advanced computational, storage and networking infrastructure provided by the Hyak supercomputer system at the University of Washington <cit.>.

§ DETAILS OF PARTICLE NUMBER ANALYSIS

In this appendix, we provide additional details for the particle number symmetry analysis, presented in Sect. <ref> of the main text.  In Fig. <ref>, we plot the error between matrix elements of $\mathcal{B}^s$, Eq. (<ref>), generated by random circuits of depth $\ell=128$, to those for a 2-design given in Eq. (<ref>) as a function of $N_{\mathcal{E}}$, demonstrating numerically that our random circuit generation scheme presented in fig:particle_number_circuit forms a 2-design for sufficiently deep circuits and number of samples $N_\mathcal{E}$. Importantly, the matrix elements in Eq. (<ref>) have a systematic dependence on the Hilbert space dimension and, in order to have a sensible comparison between different systems,
we compute the relative error between any matrix element of $\mathcal{B}^s$, Eq. (<ref>), normalized by $d_s(d_s-1)$ where $d_s$ is the dimension of the relevant symmetry block. Because of the otherwise prohibitive cost, we compute the error by only averaging over a sample of all matrix indices of $\mathcal{B}^s$, and we vary the index sample size until we see convergence; for 4, 6, 8, and 10 qubits,  we uniformly randomly sampled 600, 900, 1200, and 2000 non-zero matrix elements in every block.
More precisely, we sample over non-zero matrix elements of $\mathcal{B}^s$, but we have explicitly verified that sampling over all index combinations, including where Eq. (<ref>) yields zero, gives nearly identical results.

The circuits used to generate Fig. <ref>, with $\ell = 128$, are in the large layer limit, safely past the point required to sampling faithfully from a 2-design on each of the symmetry blocks.
We have confirmed this by comparing with unitaries directly drawn from CUE matrices in every symmetry-block. Fig. <ref> illustrates this for $N=10$ lattice sites, comparing $\ell=4,16,32,128$. While perfect scaling is observed for $\ell=16,32,128$, the $\ell=4$ circuit is too shallow, showing saturation even when $N_{\mathcal{E}}$ is further increased. Because we only have data available for $n=2,4,6,8,10$, the precise scaling with qubits $n$ of required circuit depth is difficult to estimate numerically; one needs to work essentially in the $N_{\mathcal{E}}\rightarrow \infty$ limit instead of relying on extrapolation in the scaling regime, see the linear extrapolation in fig:scalingvslayers which fails for small $\ell$. The numerical cost to do so was prohibitive, thus in fig:layerdep2 we show instead results of a coarser analysis where
$\epsilon=5\cdot 10^{-1}$ (horizontal dotted lines in fig:scalingvslayers), not relying on extrapolation and allowing us to determine the saturation value of $\ell$.  We plot in fig:layerdep2 the circuit depth required, averaged over all symmetry blocks, for a given number of qubits $n$. While results for small $n$  appear approximately block-independent,  we see significant block-dependence for larger $n=8,10$, where smaller sectors require fewer layers, represented by the error bars in fig:layerdep2. Whether this is due to an actual (sub-leading) $d_s$ dependence cannot be answered from our limited system size data, but
strongly motivates analytically deriving performance guarantees in future work, similar to symmetry-ignorant designs where they are known.

Average error $\epsilon$ in 2-design matrix elements (Eq. (<ref>)) per particle number symmetry sector, normalized by $d_s(d_s^2-1)$, versus number of unitaries $N_\mathcal{E}$ sampled from the particle number symmetry-respecting random circuits with different numbers of layers $\ell\in\{4,16,32,128\}$ (triangles). We compare to directly sampling from CUE within each symmetry sector (squares). Results shown are for of $N=10$ sites and the error is averaged over 2000 non-zero matrix elements. Dimensions of the sectors are given in parenthesis. Extrapolation of the fits  is used to produce fig:scalingpn of the main text for $\ell=128$.

(a) Circuit layers $\ell$ required to obtain a 2-design with error $5\cdot10^{-1}$, as a function of $n$ and approximately in the infinite sample limit. Error bars indicate the deviations between symmetry blocks for fixed $n$. The dotted line is not a fit, but indicates a linear behavior $\sim n$,  typically expected also from symmetry-ignorant designs.

We note that for $\ell \rightarrow \infty$, the scaling with $N_\mathcal{E}$ is independent of the circuit construction, as every block simply represents a random CUE matrix in this limit. Because of this, the scaling observed with $d_s$ equally applies if $d_s$ is replaced by the full Hilbert space $d_{\mathcal{H}}$, i.e. the size relevant for a symmetry-ignorant scheme.   We use this fact and the fit in Eq. <ref> to generate the inset of fig:scalingpn(a) showing the relative gain $r_s$ in number of circuit samples for the symmetry-conscious over a symmetry-ignorant scheme.

The sample cost scaling shown in Fig. <ref>(b) is entirely due to the number of shots per circuit $N_M$. In Fig. <ref>(b) we fix the number of ensembles ($N_\mathcal{E} =$ 1,428) for each data point and work in the large layer limit, $\ell=128$. Fig. <ref> demonstrates why this is a reasonable approach. In Fig. <ref>(a), we plot the number of ensembles $N_\mathcal{E}$ required to estimate the 2- and 3-purities to a relative error of $\epsilon = 0.05$ in the infinite shot limit $N_M \to \infty$. The necessary number of ensembles appears to peak and then decrease as a function of block dimension (though this trend cannot continue indefinitely, as one must always implement at least one random unitary). We expect this behavior to saturate for large enough block dimension, as the variance of the infinite shot purity estimator approaches a constant in the large Hilbert space dimension limit <cit.>. More practically, Fig. <ref> shows that our choice of ($N_\mathcal{E} =$ 1,428) in Fig. <ref>(b) is well beyond the number of ensembles needed to predict the 2-purity to within 5 percent for all cases considered, so we can be confident in attributing the trend in sample cost to the number of shots required per ensemble.

As a sanity check on the $N_\mathcal{E}$ scaling, in Fig. <ref>(b) we plot the true $k$-purities for the states used to create Fig. <ref>(b) and Fig. <ref>(a). Note that the states, even for large Hilbert space dimension, have purities of order one. As the cost in $N_\mathcal{E}$ for purity estimation is expected to be largest in the pure state case <cit.>, this shows that the trend in Fig. <ref>(a) is not simply because states at large $N$ (large $d_s$) are less pure.

(a) Actual $k$-purities per symmetry block, labeled by particle number $n$, for the (normalized) states used in Fig. <ref>b. (b) The number of ensembles $N_{\mathcal{E}}$ required, in the infinite shot limit,  to estimate the 2-purity and 3-purity of a single block to fixed percent error of 5 percent.

§ DETAILS OF THE (1+1)D Z2 LGT EXAMPLE

Average error $\epsilon$ in 2-design matrix elements (Eq. (<ref>)), normalized by $d_s(d_s^2-1)$, versus number of random unitaries $N_\mathcal{E}$  sampled from either the symmetry-respecting random circuits of Section <ref> of depth $\ell=32$ (triangles) or from direct sampling from the CUE within each symmetry sector (squares) for $\mathbb{Z}_{2}^{1+1}$ LGT for a subsystem of size $N_A=5$ ($N=10$). Number of indices sampled is 900.

In the main text, we presented circuits forming symmetry-conscious $k$-designs for $\Ztwo$ LGT in (1+1)d with matter. We demonstrated that these circuits allow to  measure $k$-purities and von Neumann entropies within each symmetry sector, as well as separately extracting the symmetry and distillable entanglement, and the symmetry-resolved Schmidt spectrum using classical shadows. A sufficient, but not necessary condition for such randomized measurement schemes to be successful is that the randomizing circuits form approximate $k$-designs. In this appendix, we explicitly demonstrate that the circuits in question do, in fact, form a sector-wise $2$-design. In particular, we show that they reproduce the correct $2$-design matrix elements (see eq:def2design) for sufficiently deep circuits.

The results are shown in Fig. <ref> for a subsystem of size $N_A=5$ (9 qubits) of a $N=10$ (matter) site system (20 qubits). As described in the main text,  Gauss laws at the entanglement boundaries lead to symmetries of  $\rho_A$ in the subsystem. We demonstrate that, within each sector, the random circuits described in Sect. <ref> form a 2-design, by computing the matrix elements of $\mathcal{B}^s$ in Eq. (<ref>) for all non-trivial symmetry sectors calculated with respect to $N_\mathcal{E}$ random circuits with $\ell=32$ layers.[Filling sectors $n_A=0$ and $n_A=N^A$ are trivial as they have unit block size and are not shown.]
We average over 900 random non-zero matrix elements and compare these results to the matrix elements we obtain from directly sampling $N_\mathcal{E}$ random unitaries CUE-randomly within each symmetry sector. The error is normalized by a dimensional factor $d_s(d^2_s-1)$ where $d_s$ is the dimension of the corresponding symmetry sector.
The results show good agreement, indicating that our circuits do indeed form 2-designs.

§ DETAILS OF THE (2+1)D Z2 LGT EXAMPLE

In this appendix we provide details of the random measurement analysis of $\mathbb{Z}_{2}^{2+1}$, discussed in the main text. The determination of $k$-purities follows exactly that in (1+1) spacetime dimensions.
We explicitly show that the (2+1)d circuits in fig:2d_setup(b) explicitly realizes a $2$-design by repeating the analysis of section <ref>. The results are compactly summarized in fig:2designtestz22d showing the deviation  from a $2$-design, eq:2designquanteq:def2design, for every symmetry sector $s$, demonstrating again convergence with increasing sample size $N_{\mathcal{E}}$, for a subsystem of size $3\times 2$ and $\ell$ layers.

Average error $\epsilon$ in 2-design matrix elements (Eq. (<ref>)), normalized by $d_s(d_s^2-1)$, versus number of random unitaries $N_\mathcal{E}$  sampled from either the symmetry-respecting random circuits of Section <ref> of depth $\ell=64$ (triangles) or from direct sampling from the CUE within each symmetry sector (squares) for $\mathbb{Z}_{2}^{2+1}$ LGT for a subsystem of size $3\times 2$ (with fixed boundary conditions in $y$). Number of indices sampled is 900.

Relative entropy between exact and shadow-reconstructed symmetry-resolved $\rho_{A,s}$ for the $\mathbb{Z}_{2}^{2+1}$ LGT ground state at $\epsilon=0.2$, with $N_x\times N_y=(3+5)\times 2 $ and fixed BC in $y$, $\ell=64$ circuit layers. Data for $N_{\mathcal{S}}=2^{16}$ are shown in fig:LGT2dEHT(b) of the main text.

The classical shadow analysis for $\mathbb{Z}_{2}^{2+1}$ ground states follows the previously discussed (1+1)d case. fig:app:KL shows
the sectorwise relative entropy,

\begin{align}\label{eq:relentr:app}
S(\bar{\rho}_{A,s}|| \bar{\sigma}_{A,s})\equiv -\text{tr}_s[ \bar{\rho}_{A,s} (\log(\bar{\rho}_{A,s})-\log(\bar{\sigma}_{A,s}) ]
\end{align}

where $\bar{\rho}_{A,s}$ and $\bar{\sigma}_{A,s}$ are the exact and shadow-reconstructed reduced density matrices (projected onto symmetry block $s$) of the $\mathbb{Z}_{2}^{2+1}$ ground state at $\epsilon=0.2$, with $N_x\times N_y=(3+5)\times 2 $, fixed BC in $y$, and $\ell=64$ layers; bars indicate normalization i.e. $\bar{\rho}_{A,s}={\rho}_{A,s}/p_s$ where $p_s = \text{tr}_s[{\rho}_{A,s}]$; $\text{Tr}_s$ denotes the trace over sector $s\in\{ \ua/\da\}^{2^{2N_y}}$. An accuracy of up to $10^{-2}-10^{-3}$ is achieved for the largest samples (where $p_s \cdot N_{\mathcal{S}}\gtrapprox 10^4$); the BW-EHT ansatz at similar cost typically reaches a precision better than $10^{-5}-10^{-6}$ for the same configuration. Convergence with increasing shadow number $N_{\mathcal{S}}$ of the shadow-reconstructed density matrix towards the exact one is evident and shows a power-law behavior consistent with the scaling of the (1+1)d case within error bars.

Left column: Symmetry resolved Schmidt spectrum $P_{s,\lambda}$, reconstructed using BW-EHT, for $\epsilon=0.075,0.3,0.5$ and with $N_x\times N_y=(3+3)\times 2$ and periodic boundary conditions in $y$ (and $x$),  $\ell=64$, $N_{\mathcal{E}}=50$, $N_{M}=1024$. Right column: Symmetry resolved Entanglement spectrum.

Finally, our Bisognano-Wichmann theorem based entanglement Hamiltonian tomography (BW-EHT) protocol follows <cit.>, except that we perform the optimization in every symmetry-sector $s$ separately. The approach is based on a representation of the reduced density matrix in terms of an Entanglement Hamiltonian,

\begin{align}
H_A = -\log[\rho_A]
\end{align}

with Schmidt representation $\rho_A = \sum_{\lambda} P_{\lambda} | \lambda \rangle \langle \lambda| $, where $P_\lambda \equiv \exp\{-\xi_\lambda\}$ and $\xi_\lambda$ are the eigenvalues of $H_A$. Because $\rho_A = \bigoplus_s \rho_{A,s}$ it follows that also $H_A = \bigoplus_s H_{A,s}$, $s$ labeling `quantum numbers' of the spectrum $\xi_{\lambda,s}$ of $H_{A,s}$. An example of the results of this analysis is shown in fig:app:BW2d, where we show the symmetry-resolved Schmidt  $P_{s,\lambda}$ and entanglement spectrum $\xi_{\lambda,s}$ of the $\mathbb{Z}_{2}^{2+1}$ ground state at $\epsilon=0.075,0.3,0.5$, for $N_x\times N_y=(3+3)\times 2$, periodic boundary conditions in $y$ and  $\ell=64$. We compare exact results (empty squares) with an ideal BW analysis  (stars) in the $N_{\mathcal{E}},N_M\rightarrow \infty$ limit, obtained by directly minimizing the sector-wise relative entropy eq:relentr:app between the exact state and the BW ansatz, where

\begin{align}
\bar{\sigma}_{A,s}[\{ \beta_{\mathcal{O},s}\}] \equiv \exp\{ -\sum_{\mathcal{O}}\beta_{\mathcal{O},s} H_{\mathcal{O},s}\}/Z_s
\end{align}

where $H_{\mathcal{O},s}$ are local operators of the physical Hamiltonian projected onto sector $s$, and $Z_s\equiv \text{tr}_s[\exp\{ -\sum_{\mathcal{O}}\beta_{\mathcal{O},s} H_{\mathcal{O},s}\}]$ so that $\text{tr}_s[\bar{\sigma}_{A,s}]=1$. Note that, for the symmetry-global reduced density matrix we could also have written,

\begin{align}
\sigma_A[\{ \beta_{\mathcal{O}}\}]\equiv \exp\{ -\sum_{\mathcal{O}}\beta_{\mathcal{O}} H_{\mathcal{O}}\}/Z\,,
\end{align}

where $H_{\mathcal{O}}$ are un-projected operators  and $Z\equiv \text{tr}_A[\exp\{ -\sum_{\mathcal{O}}\beta_{\mathcal{O}} H_{\mathcal{O}}\}]$, see e.g. <cit.> where a symmetry-ignorant BW analysis of this model was performed. Noting that $\sigma_A = \bigotimes_s \sigma_{A,s}$,
a BW parameterization of $\sigma_{A,s}$ should only be in terms of one global set $\{\beta_{\mathcal{O}}\}$. For convenience however, we work with a BW optimzation based on normalized $\bar{\sigma}_{A,s}={\sigma}_{A,s}/p_s$, so that in practice we obtain (slighly) different $\{ \beta_{\mathcal{O},s}\}$ in every sector. Further, fig:app:BW2d also shows results of the BW analysis for $N_{\mathcal{E}}$ and in the infinite shot limit, $N_M\rightarrow\infty$, (filled diamonds) versus
$N_M=1024$ (empty diamonds), the latter are shown in the main text. The combined statistical error due to finite $N_{\mathcal{E}}$ and $N_M$ (displayed as error bars in fig:LGT2dEHT(c)) is small for the parameters discussed in this manuscript. While the realization of a $k$-design is not required for both shadow and BW analysis, we have  verified through direct comparison with sampling from a CUE matrix in every sector $s$ that this is the case, we work in a regime where the circuit depth $\ell$ is deep enough to ensure sufficient $k$-design randomization over every symmetry sector $s$.