Statistical correlations between locally randomized measurements: a toolbox for probing entanglement in many-body quantum states

TL;DR

This paper introduces a comprehensive framework for using statistical correlations of randomized measurements to probe entanglement and quantum states in many-body systems, with applications demonstrated in quantum simulators.

Contribution

It develops a general theoretical approach for measurement protocols using randomized measurements, extending previous work to access state overlaps and tomography.

Findings

Successfully measured second Rényi entropy in a 10-qubit trapped-ion system.

Generalized protocols to access quantum state overlaps.

Discussed scalable quantum state tomography methods.

Abstract

We develop a general theoretical framework for measurement protocols employing statistical correlations of randomized measurements. We focus on locally randomized measurements implemented with local random unitaries in quantum lattice models. In particular, we discuss the theoretical details underlying the recent measurement of the second R\'{e}nyi entropy of highly mixed quantum states consisting of up to $10$ qubits in a trapped-ion quantum simulator [Brydges et al., Science 364, 260 (2019)]. We generalize the protocol to access the overlap of quantum states, prepared sequentially in an experiment. Furthermore, we discuss proposals for quantum state tomography based on randomized measurements within our framework and the respective scaling of statistical errors with system size.

Figures

Equations

$$S_2({\rho }_A)=-{\log }_2\text{Tr}\left[{\rho }_A^2\right].$$

$$\text{ or}\begin{array}{rl}{S}_2\left({\rho }_A\right)& >S_2\left({\rho }_{A\cup B}\right)\\ S_2\left({\rho }_B\right)& >S_2\left({\rho }_{A\cup B}\right),\end{array}$$

$$\text{Tr}\left[{\rho }_A^2\right]=\text{Tr}\left[\mathbb{S}{\rho }_A\otimes {\rho }_A\right]$$

$$\text{Tr}\left[{\rho }_A^2\right]$$

$$=({\mathcal{D}}_A+1)\sum _{\mathbf{s}}\overline{P_U({\mathbf{s}}_A{)}^2}-1,$$

$$\text{Tr}\left[{\rho }_A^2\right]=d^{N_A}\sum _{{\mathbf{s}}_A,{\mathbf{s}}_A^{\prime }}(-d{)}^{-D[{\mathbf{s}}_A,{\mathbf{s}}_A^{\prime }]}\overline{P_U({\mathbf{s}}_A)P_U({\mathbf{s}}_A^{\prime })},$$

$$\text{Tr}\left[{\rho }^2\right]=2\left(\overline{P_U(\uparrow {)}^2+P_U(\downarrow {)}^2-P_U(\uparrow )P_U(\downarrow )}\right)$$

$$\rho =\frac{1}{4}\sum _{\mu ,\nu =0}^3{r}_{\mu \nu }{\sigma }_{\mu }\otimes {\sigma }_{\nu },$$

$$r=\left(\begin{array}[]{c|ccc}1&r_{01}&r_{02}&r_{03}\\ \hline\cr r_{10}&r_{11}&r_{12}&r_{13}\\ r_{20}&r_{21}&r_{22}&r_{23}\\ r_{30}&r_{31}&r_{32}&r_{33}\end{array}\right)\equiv\left(\begin{array}[]{cc}1&\mathbf{w}^{T}\\ \mathbf{v}&R\end{array}\right),$$

$$\text{Tr}\left[{\rho }^2\right]=\frac{1}{4}\left(1+|\mathbf{v}{|}^2+|\mathbf{w}{|}^2+||R|{|}^2\right),$$

$$\begin{array}{rl}\mathbf{v}& \to {\mathbf{v}}^{\prime }=Q_1\mathbf{v}\\ \mathbf{w}& \to {\mathbf{w}}^{\prime }=Q_2\mathbf{w}\\ R& \to R^{\prime }=Q_{1}R{Q}_2^{†},\end{array}$$

$$Z_U^{(1)}$$

$$=P_U(\uparrow \uparrow )+P_U(\uparrow \downarrow )-P_U(\downarrow \uparrow )-P_U(\downarrow \downarrow )$$

$$Z_U^{(1,2)}$$

$$=P_U(\uparrow \uparrow )-P_U(\uparrow \downarrow )-P_U(\downarrow \uparrow )+P_U(\downarrow \downarrow )$$

$$={\left(Q_{1}R{Q}_2^{†}\right)}_{33}.$$

$$\text{Tr}\left[{\rho }^2\right]$$

$$=4\sum _{\stackrel{\mathbf{s},{\mathbf{s}}^{\prime }=}{\uparrow \uparrow ,\uparrow \downarrow ,\downarrow \uparrow ,\downarrow \downarrow }}(-2{)}^{-D[\mathbf{s},{\mathbf{s}}^{\prime }]}\overline{P_U(\mathbf{s})P_U({\mathbf{s}}^{\prime })},$$

$${\int }_{\text{Haar}}\text{d}U=1,$$

$${\int }_{\text{Haar}}\text{d}Uf(VU)={\int }_{\text{Haar}}\text{d}Uf(UV)={\int }_{\text{Haar}}\text{d}Uf(U).$$

$${\Phi }^{(k)}(O)={\int }_{\text{Haar}}\text{d}U{\left(U^{†}\right)}^{\otimes k}O{U}^{\otimes k}.$$

$$W_{\pi }=\sum _{s_1,\dots ,s_k=1}^d|s_{\pi (1)}⟩\cdots |s_{\pi (k)}⟩⟨s_1|\cdots ⟨s_k|.$$

$${\Phi }^{(k)}(O)=\sum _{\pi ,\sigma \in {\mathcal{S}}_k}{C}_{\pi ,\sigma }\text{Tr}\left[W_{\sigma }O\right]W_{\pi },$$

$${\Phi }^{(1)}(O)=\frac{1}{d}\text{Tr}\left[O\right]$$

$$\begin{array}{rl}{\Phi }^{(2)}(O)=\frac{1}{d^2-1}& (1\text{Tr}\left[O\right]+\mathbb{S}\text{Tr}\left[\mathbb{S}O\right]\\ & -\frac{1}{d}\mathbb{S}\text{Tr}\left[O\right]-\frac{1}{d}1\text{Tr}\left[\mathbb{S}O\right]),\end{array}$$

$${\Phi }_{\mathcal{E}}^{(k)}(O)={\int }_{\mathcal{E}}\text{d}U{\left(U^{†}\right)}^{\otimes k}O{U}^{\otimes k},$$

$${\Phi }_{\mathcal{E}}^{(k)}(O)=\frac{1}{|\mathcal{E}|}\sum _{U\in \mathcal{E}}{\left(U^{†}\right)}^{\otimes k}O{U}^{\otimes k}.$$

$${\Phi }_N^{(k)}(O)\equiv \overline{{\left(U^{\otimes k}\right)}^{†}O{U}^{\otimes k}},$$

$${\Phi }_N^{(k)}(O)=\sum _{\pi ,\sigma \in {\mathcal{S}}_k^{\otimes N}}{C}_{\pi ,\sigma }{W}_{\pi }\text{Tr}\left[W_{\sigma }O\right].$$

$${\Phi }_N^{(k)}\left(\underset{i=1}{\overset{N}{⨂}}{O}_i\right)$$

$$=\underset{i=1}{\overset{N}{⨂}}{\Phi }_1^{(k)}(O_i)$$

$$=\underset{i=1}{\overset{N}{⨂}}\sum _{{\pi }_i,{\sigma }_i\in S_k}{C}_{{\pi }_i,{\sigma }_i}{W}_{{\pi }_i}\text{Tr}\left[W_{{\sigma }_i}{O}_i\right]$$

$$\sum _{\mathbf{s},{\mathbf{s}}^{\prime }}$$

$$=\text{Tr}\left[\sum _{\mathbf{s},{\mathbf{s}}^{\prime }}{O}_{\mathbf{s},{\mathbf{s}}^{\prime }}|\mathbf{s}⟩⟨\mathbf{s}|\otimes |{\mathbf{s}}^{\prime }⟩⟨{\mathbf{s}}^{\prime }|\overline{U^{\otimes 2}\rho \otimes \rho {\left(U^{†}\right)}^{\otimes 2}}\right]$$

$$=\text{Tr}\left[O{\Phi }_N^{(2)}(\rho \otimes \rho )\right]$$

Full text

I Introduction

The development of intermediate- and large scale quantum simulators Preskill (2018), consisting of tens of individually controlled quantum particles, requires new tools to probe and verify complex many-body quantum systems Bloch et al. (2012); Blatt and Roos (2012); Browaeys et al. (2016); Gambetta et al. (2017). A key feature of composite quantum systems, in particular quantum lattice models, is bipartite entanglement which can be accessed by the measurement of Rényi entropies Horodecki et al. (2009). In spin models with a few degrees of freedom, Rényi entropies can be determined from tomographic reconstruction of the quantum state of interest Häffner et al. (2005); Gross et al. (2010); Lanyon et al. (2017); Torlai et al. (2018). In systems realizing one-dimensional Bose Hubbard models, the measurement of the second Rényi entropy has been demonstrated in remarkable experiments Islam et al. (2015); Kaufman et al. (2016). Here, two identical copies of the quantum state have been prepared and the second Rényi entropy has been determined from an interference experiment Bovino et al. (2005); Daley et al. (2012); Pichler et al. (2013).

In Ref. Brydges et al. (2019) we have demonstrated in a theory-experiment collaboration a new protocol in which the second order Rényi entropy is inferred from statistical correlations of locally randomized measurements Elben et al. (2018). Here, a spin model in a trapped-ion quantum simulator was realized and the generation of entanglement during quench dynamics was monitored. The measurement protocol to access the second Rényi entropy was based on only local operations on individual spins and a single instance of the quantum state. It required, albeit an exponential scaling with the number of degrees of freedom, a significantly lower number of measurements than standard quantum state tomography Brydges et al. (2019). This measurement protocol is thus immediately applicable in a broad class of quantum simulators, realizing spin models, with single site readout and control. In particular, we have in mind systems based on trapped-ions Zhang et al. (2017); Brydges et al. (2019), Rydberg atoms Zeiher et al. (2017); Barredo et al. (2018); Guardado-Sanchez et al. (2018); Keesling et al. (2019) and superconducting qubits Blumoff et al. (2016); Barends et al. (2016); Otterbach et al. ; Gong et al. (2019) in arbitrary spatial dimensions. Moreover, our protocol can straightforwardly be applied to extended systems such as quantum networks. In this context, it has also been proposed use randomized measurements to detect genuine multipartite correlations Ketterer et al. (2019).

The key ingredient for the protocols described in this paper are statistical correlations between randomized measurements. The goal of this article is to develop a general mathematical formalism to evaluate such correlations and, equipped with this toolbox, elaborate on the theoretical details behind the protocol realized in Ref. Brydges et al. (2019).

In general, a random measurement on a (reduced) quantum state $\rho$ is performed by the application of a random unitary $U$, sampled from an appropriate ensemble (see below), and the subsequent measurement of the expectation value $\langle O\rangle_{U}=\textrm{Tr}\left[U\rho U^{\dagger}O\right]$ of a fixed observable $O$. In this paper, we focus on spin models and locally randomized measurements where the random unitaries are of the form $U=\otimes_{i}U_{i}$ with the $U_{i}$ independent random spin rotations sampled from unitary designs (in particular the circular unitary ensemble (CUE)) Elben et al. (2018); Brydges et al. (2019). Furthermore, we provide an in-depth comparison with protocols based on globally randomized measurements van Enk and Beenakker (2012); Elben et al. (2018) where global random unitaries $U$ are sampled from unitary designs defined on the entire Hilbert space. These global random unitaries can be generated in interacting quantum lattice models with engineered disorder using random quenches Nakata et al. (2017); Elben et al. (2018); Vermersch et al. (2018) and the corresponding protocols are in particular relevant for atomic Hubbard models Elben et al. (2018); Vermersch et al. (2018).

In the second part, we extend the formalism to derive a protocol to measure the overlap $\textrm{Tr}\left[\rho\rho^{\prime}\right]$ of two states $\rho$ and $\rho^{\prime}$ which are prepared sequentially in one experiment. This allows in particular to directly measure the many-body Loschmidt echo Goussev et al. , without implementing time-reversed operations. Finally, we discuss within our formalism a proposal of Ohliger et al. Ohliger et al. (2013) to use randomized measurements, implemented with global random unitaries, to perform full quantum state tomography in atomic Hubbard models. We generalize this protocol to local random unitaries available in spin models and investigate in detail numerically how the required number of measurements to reconstruct the density matrix $\rho$ up to a fixed error scales with system size.

II Measurement of the second Rényi entropy

In this section, we discuss the measurement of the second Rényi entropy in quantum lattice models. After a short review, we focus on protocols utilizing randomized measurements, give explicit recipes and discuss examples. We consider a lattice system $\mathcal{S}$ described by a quantum state $\rho$. The second Rényi entropy $S_{2}(\rho_{A})$ of the reduced density matrix $\rho_{A}=\textrm{Tr}_{\mathcal{S}\backslash A}\left[\rho\right]$ of an arbitrary subsystem $A\subseteq\mathcal{S}$ consisting of $N_{A}$ sites is defined as

[TABLE]

Using $S_{2}(\rho_{A})$, one shows that bipartite entanglement exists between two disjoint subsystems $A$ and $B$ of $\mathcal{S}$ with reduced density matrices $\rho_{A}=\textrm{Tr}_{\mathcal{S}\backslash A}\left[\rho\right]$ and $\rho_{B}=\textrm{Tr}_{\mathcal{S}\backslash B}\left[\rho\right]$ 111We note that this criterium is sufficient but not necessary. if

[TABLE]

where $\rho_{A\cup B}=\textrm{Tr}_{\mathcal{S}\backslash A\cup B}\left[\rho\right]$ is the reduced density matrix of $A\cup B$ Horodecki et al. (2009).

To measure $S_{2}(\rho_{A})$, i.e.  the purity $\textrm{Tr}\left[\rho^{2}_{A}\right]$ of a (reduced) density matrix $\rho_{A}$, various protocols have been proposed and realized which we shortly review in the following: A first option realized in spin models is to perform full quantum state tomography of $\rho_{A}$ Häffner et al. (2005); Gross et al. (2010); Lanyon et al. (2017); Torlai et al. (2018). However, due to the exponential scaling of the number measurement settings, at least $d^{2N_{A}}$ for $N_{A}$ spins and standard tomography Gross et al. (2010), this approach is limited to systems with a few degrees of freedom Häffner et al. (2005). On the contrary, recent efficient tomographic methods require a specific structure of the state of interest Lanyon et al. (2017); Torlai et al. (2018).

A second class of protocols Daley et al. (2012); Pichler et al. (2013) is based on noting that the purity

[TABLE]

can be obtained from measuring the expectation value of the swap operator $\mathbb{S}$ acting on two copies of a quantum state $\rho_{A}$ Bovino et al. (2005); Horodecki (2003). Here, $\mathbb{S}$ is defined by $\mathbb{S}\ket{s_{A}}\otimes\ket{s_{A}^{\prime}}=\ket{s_{A}^{\prime}}\otimes\ket{s_{A}}$ for any two states $\ket{s_{A}},\ket{s_{A}^{\prime}}$. Prerequisite for these protocols is thus the experimental ability to create two identical copies $\rho_{A}\otimes\rho_{A}$ and to perform joint operations on them to measure $\textrm{Tr}\left[\mathbb{S}\rho_{A}\otimes\rho_{A}\right]=\textrm{Tr}\left[\rho_{A}^{2}\right]$. Despite the experimental complexity of this task, the purity of a quantum state of up to six particles, realizing a one-dimensional Bose-Hubbard model, has been measured in remarkable experiments Islam et al. (2015); Kaufman et al. (2016). Recently, the method has been also transferred to trapped-ion quantum simulator, and applied to a one-dimensional system of five qubits Linke et al. (2018). Creating identical copies of a quantum state in larger systems and higher spatial dimensions remains however a significant technological challenge.

Only a single instance of a quantum state is required in a third class of protocols van Enk and Beenakker (2012); Elben et al. (2018); Vermersch et al. (2018); Brydges et al. (2019), utilizing statistical correlations of randomized measurements, which we discuss in the remainder of this paper.

II.1 Second Rényi entropy from statistical correlations of randomized measurements

In this subsection, we describe the protocol to estimate the purity $\textrm{Tr}\left[\rho^{2}_{A}\right]$ of a (reduced) density matrix $\rho_{A}$ from statistical correlation of randomized measurements. While protocols based on global random unitaries are also applicable to atomic Hubbard models Elben et al. (2018); Vermersch et al. (2018), we focus in this paper on spin models consisting of $N$ spins with local Hilbert space $\mathscr{h}$ of dimension $d$ (i.e. $N$ qudits). Here, in addition, experimentally simpler local random unitaries are available (see below). A schematic view of the experimental sequence we have in mind is displayed in Fig. 1 (see also Ref. Brydges et al. (2019); Elben et al. (2018)). An (entangled) quantum state of interest $\rho$ is for instance prepared via quench dynamics originating from a simple initial state. The experimental protocol to measure the purity of the reduced density matrix $\rho_{A}$ of a subsystem $A$ of $N_{A}$ qudits consists then of several steps. First, one applies to $\rho_{A}$ a random unitary $U_{A}$. This can (i) either be a global random unitary sampled from a unitary $2$-design Gross et al. (2007) defined on the entire Hilbert space $\mathcal{H}_{A}=\mathscr{h}^{\otimes N_{A}}$ with dimension $\mathcal{D}_{A}=d^{N_{A}}$ of the subsystem, or (ii) a local unitary of the form $U_{a}=\bigotimes_{i\in A}U_{i}$ where the $U_{i}$ are sampled independently from a unitary $2$-design defined on the local Hilbert space $\mathscr{h}$. We compare both approaches in detail below. Subsequently, a measurement in the computational basis is performed, which is repeated with the same random unitary $U_{A}$ to estimate the occupation probabilities $P_{U}(\mathbf{s}_{A})=\textrm{Tr}\left[U_{A}\rho U_{A}^{\dagger}\ket{\mathbf{s}_{A}}\bra{\mathbf{s}_{A}}\right]$ of computational basis states $\ket{\mathbf{s}_{A}}=\ket{s_{1},\dots,s_{N_{A}}}$ ($s_{i}=1\dots d$ for $i\in A$). In a second step, this is repeated for many different random unitaries, to estimate the average over the ensemble of random unitaries.

Given the set of outcome probabilities $P_{U}(\mathbf{s}_{A})$ for the computational basis states $\mathbf{s}_{A}$, one estimates the purity of $\rho_{A}$ from second-order cross correlations across the random unitary ensemble. For global random unitaries [case (i)], one finds

[TABLE]

where the “global” Hamming distance is defined as $D_{G}[\mathbf{s}_{A},\mathbf{s}^{\prime}_{A}]=0$ if $\mathbf{s}_{A}=\mathbf{s}_{A}^{\prime}$ and $D_{G}[\mathbf{s}_{A},\mathbf{s}_{A}^{\prime}]=1$ if $\mathbf{s}_{A}\neq\mathbf{s}_{A}^{\prime}$. The expression in the second line has first been given in Ref. van Enk and Beenakker (2012). If independent local random unitaries on individual qudits are used [case (ii)], the purity is estimated from

[TABLE]

where the Hamming distance $D[\mathbf{s}_{A},\mathbf{s}_{A}^{\prime}]$ between two states $\ket{\mathbf{s}_{A}}=\ket{s_{1},\dots,s_{N_{A}}}$ and $\ket{\mathbf{s}_{A}^{\prime}}=\ket{s^{\prime}_{1},\dots,s^{\prime}_{N_{A}}}$ is defined as the number of local constituents $i\in A$ where $s_{i}\neq{s}^{\prime}_{i}$, i.e.  $D[\mathbf{s}_{A},{\mathbf{s}_{A}}^{\prime}]\equiv\#\left\{i\in A\,|\,s_{i}\neq{s}^{\prime}_{i}\right\}$. Equation (5) has first been obtained in Ref. Brydges et al. (2019), and represents an explicit version of the recursive formula given in Ref. Elben et al. (2018). In Ref. Brydges et al. (2019), the identity has been established directly using moments of matrix elements of random unitaries. In this paper, we develop a general formalism to evaluate statistical correlations of randomized measurements based on the Weingarten calculus Collins and Nechita (2010). This approach enables a constructive and simple derivation of Eq. (5) and, in addition, it simplifies the development of new protocols (see below).

An intriguing connection to the previous works realizing the swap operator on two physical copies Daley et al. (2012); Pichler et al. (2013); Islam et al. (2015); Kaufman et al. (2016) can be seen as follows: We can rewrite any product of outcome probabilities $P_{U}(\mathbf{s}_{A})P_{U}(\mathbf{s}_{A}^{\prime})=\textrm{Tr}_{\mathcal{H}^{\otimes 2}_{A}}\left[U_{A}^{\otimes 2}\rho_{A}^{\otimes 2}U_{A}^{\dagger\otimes 2}\ket{\mathbf{s}_{A}}\bra{\mathbf{s}_{A}}\otimes\ket{\mathbf{s}_{A}^{\prime}}\bra{\mathbf{s}_{A}^{\prime}}\right]$ as the expectation value of an operator $\ket{\mathbf{s}_{A}}\bra{\mathbf{s}_{A}}\otimes\ket{\mathbf{s}_{A}^{\prime}}\bra{\mathbf{s}_{A}^{\prime}}$ on the doubled Hilbert space $\mathcal{H}^{\otimes 2}_{A}$. Using this, we can intuitively understand the ensemble average over second order cross correlations taken in Eqs. (4) and (5) as an effective construction of the swap operator on two virtual copies of $\rho_{A}$ (see for details Sec. III).

The approaches with global and local random unitaries differ in various aspects: First, the implementation of global random unitaries from a unitary $2$-design acting on the entire many-body quantum state $\rho_{A}$ requires interactions between the particles. It has been proposed to prepare them efficiently in quantum circuits using (random) entangling gates Dankert et al. (2009) or in generic interacting many-body systems using time evolution subject to random quenches Ohliger et al. (2013); Nakata et al. (2017); Elben et al. (2018); Vermersch et al. (2018). On the contrary, local random unitaries, available in spin models, are single “qudit” operations (random spin rotations) which have been demonstrated with high fidelity and repetition rates Brydges et al. (2019). Second, the protocol utilizing local random unitaries allows, from a single experimental dataset obtained from randomized measurements on the subsystem $A$, to estimate the purity $\textrm{Tr}\left[\rho_{A^{\prime}}^{2}\right]$ of the reduced density matrix $\rho_{A^{\prime}}$ of any subsystem $A^{\prime}\subseteq A$. To this end, Eq. (5) is evaluated with occupation probabilities $P_{U}(\mathbf{s}_{A^{\prime}})$ of states $\ket{\mathbf{s}_{A^{\prime}}}$ of the logical basis of $A^{\prime}$. Third, the two protocols differ in their sensitivity to statistical errors. In an experiment, statistical errors of the estimated purity arise from a finite number $N_{U}$ of applied unitaries and a finite number $N_{M}$ of measurements per random unitary (projection noise). The total number of measurements $N_{M}N_{U}$ scales exponentially with the number $N_{A}$ of degrees of freedom in the subsystem $A$, with exponents significantly smaller than in full quantum state tomography Brydges et al. (2019); Elben et al. (2018). As we discuss in detail below, the protocol utilizing global unitaries is, for pure product states, favorable in terms of statistical errors.

To prove Eqs. (4) and (5), we introduce results of the theory of random unitaries in Sec. III. The proof follows then in Sec. IV.1.

II.2 Illustrative examples

In the remainder of this section, we illustrate Eqs. (4) and (5) using simple examples.

II.2.1 Single qubit

The density matrix $\rho=\frac{1}{2}\left({\mathbb{1}}_{2}+\mathbf{v}\cdot\mathbf{\sigma}\right)$ of a single qubit is conveniently represented on the Bloch sphere, with the real Bloch vector $(\mathbf{v})_{i}=\textrm{Tr}\left[\rho\sigma_{i}\right]$ and $\mathbb{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ the Pauli matrices. The purity $\textrm{Tr}\left[\rho^{2}\right]=\frac{1}{2}\left(1+|\mathbf{v}|^{2}\right)$ is thus fully determined by the length of the Bloch vector $|\mathbf{v}|$. Our approach to estimate the length of the Bloch vector consists in applying a random unitary $U$ to the state $\rho$ and measuring the difference of occupation $Z_{U}=P_{U}(\uparrow)-P_{U}(\downarrow)=\textrm{Tr}\left[U\rho U^{\dagger}\sigma_{z}\right]=\left(Q_{U}\mathbf{v}\right)_{3}$ of the computational basis states $\ket{\uparrow},\ket{\downarrow}$. Here, $Q_{U}$ is the unique rotation matrix corresponding to the unitary $U$, i.e.  $U\mathbf{v}\cdot\sigma U^{\dagger}=\left(Q_{U}\mathbf{v}\right)\cdot\sigma$ for all Bloch vectors $\mathbf{v}$ Gamel (2016). This is now repeated with different random unitaries, to sample the distribution of $Z_{U}$ across the circular unitary ensemble. Since the random variable $Z_{U}=\left(Q_{U}\mathbf{v}\right)_{3}$ describes the $z$-component of the rotated Bloch vector $Q_{U}\mathbf{v}$ subject to random rotations $Q_{U}$, it is intuitively apparent that its range contains information about the length of $\mathbf{v}$ (see also Fig. 2).

Formally, we note that moments $\overline{Z_{U}^{n}}$ ($n\in\mathbb{N}$) of the random variable $Z_{U}$ are invariant under unitary transformations due to the invariance properties of the Haar measure [see below, Eq. (14)]. Thus, they must be determined by properties of the Bloch vector $\mathbf{v}$ which are invariant under arbitrary rotations. The squared length $|\mathbf{v}|^{2}$ is its unique second order invariant, and thus we conclude $\overline{Z_{U}^{2}}=\overline{\left(Q_{U}\mathbf{v}\right)^{2}_{3}}\sim|\mathbf{v}|^{2}$. Indeed, using the $n$-design properties, we find that $Z_{U}$ is uniformly distributed, with zero mean and variance $\overline{Z_{U}^{2}}=|\mathbf{v}|^{2}/3$, such that $\textrm{Tr}\left[\rho^{2}\right]=\frac{1}{2}\left(1+3\overline{(P_{U}(\uparrow)-P_{U}(\downarrow))^{2}}\right)$. Inserting that $1=(P_{U}(\uparrow)+P_{U}(\downarrow))^{2}$ we can bring this into a more symmetric form to arrive at

[TABLE]

which corresponds to Eq. (4) [Eq. (5)] for the special case of $\mathcal{D}_{A}=2$ [$d=2$ and $N_{A}=1$].

II.2.2 Two qubits

We consider now the case of two qubits where randomized measurements are implemented using independent local random unitaries $U=U_{1}\otimes U_{2}$. Generalizing the single-qubit case, the two-qubit system is conveniently represented in the basis of Pauli strings

[TABLE]

where $\sigma_{0}={\mathbb{1}}_{2}$ and $\sigma_{i}$ the Pauli matrices ($1\leq i\leq 3$). The real coefficients $r_{\mu\nu}=\textrm{Tr}\left[\rho\sigma_{\mu}\otimes\sigma_{\nu}\right]$ constitute the Bloch matrix $\mathbf{r}$ (generalizing the Bloch vector)

[TABLE]

where $v_{i}=r_{i0}$, $w_{j}=r_{0j}$, and $R_{ij}=r_{ij}$. The vectors $\mathbf{v}$ and $\mathbf{w}$ are the Bloch vectors of the reduced density matrices $\rho_{1}=\textrm{Tr}_{2}\left[\rho\right]$ and $\rho_{2}=\textrm{Tr}_{1}\left[\rho\right]$ of the individual qubits, respectively. The matrix $R$ quantifies the correlations between the two qubits Gamel (2016). Using these definitions, the purity of the density matrix $\rho$ is then given by

[TABLE]

where $||R||^{2}=\textrm{Tr}\left[R^{\dagger}R\right]$.

Our aim is to estimate the purity $\textrm{Tr}\left[\rho^{2}\right]$ using random unitaries of the $U_{1}\otimes U_{2}$ and (collective) measurements in the computational basis $\ket{\uparrow\uparrow},\ket{\uparrow\downarrow},\ket{\downarrow\uparrow},\ket{\downarrow\downarrow}$. First, we note that under unitary transformation $\rho\rightarrow U_{1}\otimes U_{2}\rho U_{1}^{\dagger}\otimes U_{2}^{\dagger}$, the individual elements of the Bloch matrix transform as

[TABLE]

there $Q_{1}$ ($Q_{2}$) is the unique rotation matrix corresponding to $U_{1}$ ($U_{2}$) Gamel (2016). Thus $|\mathbf{v}|^{2}$ and $|\mathbf{w}|^{2}$ can be estimated from single qubit measurements, i.e. from the variances of the distributions

[TABLE]

and $Z^{(2)}_{U}$, respectively. The unique element of the (transformed) correlation matrix accessible from measurements solely in the computational basis is

[TABLE]

Again, due to the Haar average, moments of the random variable $Z^{(1,2)}_{U}$ must be invariant transformations (10) with arbitrary rotations $Q_{1}$ and $Q_{2}$, and a second order invariant of this type is the matrix norm $||R||_{2}$. Indeed, we find $\overline{\left(Z^{(1,2)}_{U}\right)^{2}}=||R||_{2}/9$. Thus we obtain,

[TABLE]

which corresponds to formula (5) for the special case of $d=2$ and $N_{A}=2$. To arrive at the last line, we used that $1=\left(P_{U}(\uparrow\uparrow)+P_{U}(\uparrow\downarrow)+P_{U}(\downarrow\uparrow)+P_{U}(\downarrow\downarrow)\right)^{2}$. We note that these arguments generalize to $N_{A}$-qubit systems whose density matrix can be parametrized by a rank $N_{A}$ “Bloch tensor”. The second moment of a suitably generalized random variable $Z^{(1,\dots,N_{A})}_{U}=\textrm{Tr}\left[\bigotimes_{i=1}^{N_{A}}U_{i}\rho_{A}\bigotimes_{i=1}^{N_{A}}U_{i}^{\dagger}\sigma_{z}^{\otimes N_{A}}\right]$ is connected to the Hilbert-Schmidt norm $||R^{(N_{A})}||^{2}$ of the rank $N_{A}$ “correlation tensor” $R^{(N_{A})}$.

III Random unitaries, diagrammatic calculus and unitary $k$-designs

In this section, we provide an overview over the mathematical tools necessary to describe randomized measurements based on random unitaries. We consider first a single qudit with Hilbert space $\mathscr{h}$ of (arbitrary) dimension $d$, which also corresponds to the case of global random unitaries (with the replacement $\mathscr{h}\rightarrow\mathcal{H}_{A}$, $d\rightarrow\mathcal{D}_{A}=d^{N_{A}}$, see Sec. II.1). In this setting, we introduce the Haar measure, defining the CUE, and derive elementary properties of the unitary twirling channel. Moreover, we introduce a graphical calculus simplifying the calculations. In the last subsection, we extend the framework to composite systems of many qudits where multiple, independent, random unitaries are applied locally.

III.1 Haar randomness

We first provide an overview over Haar random unitaries and introduce the central object of our formalism, the unitary twirling channel. We follow the treatment of Ref. Roberts and Yoshida (2017). An orthonormal basis of $\mathscr{h}$ is denoted with $\left\{\ket{{s}}\right\}$.

The most important ingredient are Haar random unitaries. These are unitary matrices which are distributed according to the probability distribution defined by the Haar measure on the unitary group 222Intuitively, Haar random unitaries are matrices with elements whose real and imaginary parts are independently distributed according to a normal distribution, with additional unitary constraints on the entire matrix.. Here, the Haar measure is the unique probability measure on the group of unitary matrices $\mathcal{U}(\mathscr{h})$ on $\mathscr{h}$ which is both left- and right-invariant, i.e. it satisfies for any function $f$ on $\mathcal{U}(\mathscr{h})$ and any unitary $V\in\mathcal{U}(\mathscr{h})$

[TABLE]

Often, $\mathcal{U}(\mathscr{h})$ equipped with the Haar measure is also called the circular unitary ensemble (CUE) Haake (2010), and we use $\overline{f(U)}\equiv\int\!\textrm{d}Uf(U)$ to denote the ensemble average over the CUE.

In the following, we consider the $k$-fold copy space $\mathscr{h}^{\otimes k}$, $k\in\mathbb{N}$, to calculate higher order moments of random unitaries. We note that this is a purely mathematical construction: $k$th-order products of outcome probabilities of randomized measurements can be viewed as expectation value an operator acting on $\mathscr{h}^{\otimes k}$, realizing thus $k$ “virtual copies”. It is a key property of any measurement protocol presented in this paper that only a single physical instance of a quantum state is required in the experiment. We define on $\mathscr{h}^{\otimes k}$ a quantum channel, the $k$-fold twirl by

[TABLE]

for any operator $O$ on $\mathscr{h}^{\otimes k}$. As a simple consequence of the invariance of the Haar measure [Eq. (14)], $\Phi^{(k)}$ forms a projector $\Phi^{(k)}(\Phi^{(k)}(O))=\Phi^{(k)}(O)$. We show in the following that its image is spanned by permutation operators $W_{\pi}$ , for permutations $\pi=(\pi(1),\dots,\pi(k))\in\mathcal{S}_{k}$ with the symmetric group $\mathcal{S}_{k}$, which are defined as

[TABLE]

These operators permute states between individual copies $W_{\pi}\ket{{s}_{1}}\cdots\ket{{s}_{k}}=\ket{{s}_{\pi(1)}}\cdots\ket{{s}_{\pi(k)}}$. Using that $[W_{\pi},V^{\otimes k}]=0$ for any $\pi\in\mathcal{S}_{k}$, it follows directly that $\Phi^{(k)}(W_{\pi})=W_{\pi}$, i.e. the permutation operators $W_{\pi}$ are invariant under the projection $\Phi^{(k)}$. Indeed, they span the total image of $\Phi^{(k)}$, which is proved with the Schur Weyl duality Roberts and Yoshida (2017). Explicitly, one finds

[TABLE]

where the coefficients $C_{\pi,\sigma}=\operatorname{Wg}(\pi\sigma^{-1})$ constitute the real-valued, symmetric Weingarten matrix $C$ determined by the Weingarten function $\operatorname{Wg}$ Collins and Nechita (2010); Puchała and Miszczak (2017); Roberts and Yoshida (2017). For $k\leq d$, $C$ is invertible, with inverse $Q\equiv C^{-1}$ and $Q_{\pi,\sigma}=d^{\sharp\textrm{cycles}(\pi\sigma)}$ Roberts and Yoshida (2017). For $k=1$ and $k=2$, we find

[TABLE]

and

[TABLE]

with the identity $\mathbb{1}=W_{(1,2)}$ and $\mathbb{S}=W_{(2,1)}=\sum_{{s}_{1}{s}_{2}}\ket{{s}_{2}}\ket{{s}_{1}}\bra{{s}_{1}}\bra{{s}_{2}}$ being the swap operator.

III.2 Diagrammatic calculus

In the previous section, we showed how to evaluate the unitary twirling channel $\Phi^{(k)}$ in terms of permutation operators. Here, we introduce a graphical calculus which enables the evaluation of arbitrary functionals of random unitaries, in particular $\Phi^{(k)}$. We follow and adapt here the treatment of Ref. Collins and Nechita (2010) (see also Ref. Brouwer and Beenakker (1996) for a similar approach). We first note that, for any $k\in\mathbb{N}$, an operator $O$ acting on $\mathscr{h}^{\otimes k}$ can be viewed as a $(k,k)$-tensor $O\in\mathscr{h}^{*\otimes k}\otimes\mathscr{h}^{\otimes k}$ with $k$ covariant and $k$ contravariant indices. As shown in Fig. 4(a), we represent tensors in the following as a box with decorations, where the number of empty (filled) symbols corresponds to the number of contravariant (covariant) indices. Similarly, ket vectors $\ket{a}\in\mathscr{h}$ can be viewed as a $(0,1)$ tensors with a single white decoration, and bra vectors $\bra{a}\in\mathscr{h}^{*}$ are $(1,0)$ tensors with a single black decoration. As discussed in the previous section, permutation operators are of central importance to evaluate twirling channels. Thus, we introduce a special graphical notation displayed in Fig. 4(b), which corresponds directly to their action on the $k$-fold copy space.

In order to evaluate the Haar average of an arbitrary diagram containing random unitaries, it turns out to be useful to consider Bell states $\sum_{i}\ket{i}\ket{i}$ [$\sum_{i}\bra{i}\bra{i}$], which represent special $(0,2)$ [$(2,0)$] tensors [Fig. 5(a). Such states allow the definition of graphs connecting two decorations of identical color [Fig. 5(b)] and, using the identity $(O^{T})_{ij}=O_{ji}=\sum_{k,k^{\prime}}\braket{i}{k}\braket{k^{\prime}}{O}{k}\braket{k^{\prime}}{j}$, the transposed matrix, given as the tensor with interchanged decorations [Fig. 5(c)].

Equipped with these definitions, we are now in the position to evaluate the ensemble average of an arbitrary diagram. We describe a general procedure Collins and Nechita (2010) applicable to any diagram. As an explicit example we present in Fig. 6 the evaluation of $\bra{i}\bra{j}\Phi^{(2)}(U)\ket{i^{\prime}}\ket{j^{\prime}}$ 333Note that for clarity, we omitted here the boxes of basis states, and just kept the decorations.. The first step is to replace all unitaries $U^{T}$ with $U$ and $U^{\dagger}$ with $U^{*}$ using the Bell states and tensors defined in Fig. 5. If in the resulting diagram, the number of unitaries $U$ does not equal the number of complex conjugates $U^{*}$, the diagram evaluates to zero. Otherwise, white decorations of boxes $U$ are connected with white decorations of boxes $U^{*}$ and black decorations of boxes $U$ are (independently) connected with black decorations of boxes $U^{*}$. Subsequently decorations and boxes of the random unitaries are removed. Given $k$ boxes $U$, there exist $(k!)^{2}$ possible ways to draw these connections. The ensemble average is obtained as a weighted sum of all resulting diagrams, with coefficients determined by the Weingarten matrix. To calculate these coefficients, one labels both the boxes $U$ and boxes $U^{*}$ with arbitrary integers $1,\dots k$. Each possibility to connect white (black) decorations is now described by a permutation $\alpha\in S_{k}$ ($\beta\in S_{k}$): If the white decoration of $U$-box $i$ is connected to the white decoration of $U^{*}$-box $j$, then $\alpha(i)\equiv j$. The coefficient of a diagram is obtained as $\operatorname{Wg}(\alpha\beta^{-1})$. In the case of Fig. 6, there exist four possibilities, and the resulting sum of diagrams is the graphical representation of Eq. (19).

III.3 Unitary $k$-designs

In the previous sections, we considered Haar random unitaries drawn from the circular ensemble for which Eq. (17) holds for arbitrary $k\in\mathbb{N}$. In applications, one is however typically interested in moments of random unitaries up to a finite (small) number $k$ Nakata et al. (2017), for instance $k=2$ for the estimation of second order Rényi entropy. Since the preparation of Haar random unitaries using, for instance, random quantum circuits requires an amount of resources scaling exponentially with system size Nakata et al. (2017), simpler ensembles, unitary $k$-designs, have been introduced Gross et al. (2007); Roy and Scott (2009); Dankert et al. (2009). These ensembles approximate Haar random unitaries in the sense that up to $k$-th order moments are identical, i.e. Eq. (17) holds, for an arbitrary operator $O$, up to a finite, fixed $k$ 444 Loosely speaking, up to the $k$th moment, $k$-designs are as random as Haar random unitaries.. To define unitary $k$-designs formally, we introduce the $k$-fold twirl with respect to a continuous ensemble $\mathcal{E}$ of unitary operators by

[TABLE]

and for a discrete ensemble $\mathcal{E}$ with cardinality $|\mathcal{E}|$ by

[TABLE]

One says that $\mathcal{E}$ forms an unitary $k$-design if and only if $\Phi^{(k)}_{\mathcal{E}}=\Phi^{(k)}$ and that $\mathcal{E}$ forms an * $\epsilon$-approximate $k$-design* if and only if $||\Phi^{(k)}_{\mathcal{E}}-\Phi^{(k)}||_{\diamond}<\epsilon$ Dankert et al. (2009). It follows directly that any $k$-design is also an $k^{\prime}$ design, for any $k^{\prime}<k$. A prime example of an exact $3$-design is the Clifford group Roy and Scott (2009). Importantly, $\epsilon$-approximate $k$-design can be prepared efficiently in local random quantum circuits Dankert et al. (2009), and generic interacting quantum simulators Nakata et al. (2017); Elben et al. (2018); Vermersch et al. (2018).

III.4 Composite systems

In this section, we generalize our treatment to composite systems consisting of ${N_{A}}$ qudits with Hilbert space $\mathcal{H}=\mathscr{h}^{\otimes N_{A}}$. For simplicity of notation we drop the subscript $A$. The basis $\{\ket{\mathbf{s}}\}$ denotes a product basis $\ket{\mathbf{s}}=\otimes_{i=1}^{N}\ket{s_{i}}$ for all $\mathbf{s}=(s_{1},\dots,s_{N})$. We consider random unitaries of the form $U=\bigotimes_{i=1}^{N}U_{i}$ where the $U_{i}$ ($i\in\{1,\dots,N\}$), acting on the individual qudits, are sampled independently from the $\text{CUE}(\mathscr{h})$ (a unitary $k$-design) defined on the local Hilbert space $\mathscr{h}$. We define a $k$-fold local twirling channel by

[TABLE]

where $\overline{\vphantom{O}\dots}$ denotes in this context the ensemble average over random unitaries of the form $U=\bigotimes_{i=1}^{N}U_{i}$. Generalizing Eq. (17), we find

[TABLE]

Here, $\pi=\bigotimes_{i=1}^{N}\pi_{i}\in\mathcal{S}^{\otimes N}_{k}$ and $\sigma=\bigotimes_{i=1}^{N}\sigma_{i}\in\mathcal{S}^{\otimes N}_{k}$ are tensor products of permutations and the corresponding operators $W_{\pi}\equiv\bigotimes W_{\pi_{i}}$ act locally on the $k$-fold copy space $\mathscr{h}^{\otimes k}$ of the individual qubits. The coefficients $C_{\pi,\sigma}\equiv\prod_{i=1}^{N}C_{\pi_{i},\sigma_{i}}$ are determined by products of elements of the Weingarten matrix $C$. To proof Eq. (23), we first note that we can expand an arbitrary operator $O$ in a product basis of the $N$ qudits $O=\sum_{\alpha}c_{\alpha}O_{\alpha_{1}}\otimes\dots\otimes O_{\alpha_{N}}$. By linearity, it suffices thus to restrict to $O$ being a tensor product of local operators, i.e. $O=\bigotimes_{i=1}^{N}O_{i}$. We find

[TABLE]

and from the last line, Eq. (23) follows.

A graphical calculus to treat composite systems with multiple independent random unitaries has been introduced in Ref. Vermersch et al. to prove the central identities of a proposal for measuring out-of-time-ordered correlators with local random unitaries. To extend the graphical calculus presented here in Sec. III.2 one introduces decorations of different types (circles, boxes, …) and connects only decorations of the same type. For an example, we refer to Fig. 7 where $\Phi^{(1)}_{2}$ is evaluated.

IV Second order correlations - Purity and overlap of quantum states

In this section, we apply our formalism to derive Eqs. (4) and (5) allowing the estimation of the purity of arbitrary quantum states using the experimental protocol given in Sec. II.1. Furthermore, we discuss an extension of the protocol given in Sec. II.1, to estimate the overlap of two distinct quantum states. Higher order functionals of the density matrix are discussed in the appendix A.

IV.1 Proof of the main result

We now prove Eqs. (4) and (5). We first note that the case (i) of global random unitaries can be viewed as a single qu$\mathcal{D}$it with dimension $\mathcal{D}_{A}$. Then, the first line of Eq. (4) follows directly from Eq. (5) by setting $N_{A}=1$ and $d\rightarrow\mathcal{D}_{A}=d^{N_{A}}$, and the second line by using $\sum_{\mathbf{s}_{A},\mathbf{s}_{A}^{\prime}\neq\mathbf{s}_{A}}P_{U}(\mathbf{s}^{\prime}_{A})=1-P_{U}(\mathbf{s}_{A})$. Thus, we consider in the following the case (ii) of independent local random unitaries applied to a composite system of $N_{A}$ qudits with arbitrary local dimension $d$. For simplicity of notation, we drop the subscript $A$.

We first note that the ensemble average of second order cross-correlation of outcome probabilities of randomized measurements can be rewritten as an expectation value of an operator $O$ acting on two “virtual copies” and the twirled state $\Phi_{N}^{(2)}({\rho^{\otimes 2}})$ [with $\Phi_{N}^{(2)}$ defined in Eq. (23)]. For arbitrary coefficients $O_{\mathbf{s},\mathbf{s}^{\prime}}$, it holds

[TABLE]

where we defined the operator $O=\sum_{\mathbf{s},\mathbf{s}^{\prime}}O_{\mathbf{s},\mathbf{s}^{\prime}}\ket{\mathbf{s}}\bra{\mathbf{s}}\otimes\ket{\mathbf{s}^{\prime}}\bra{\mathbf{s}^{\prime}}$ and used the self-duality of the channel $\Phi^{(2)*}_{N}=\Phi^{(2)}_{N}$. Secondly, we observe that, for an arbitrary quantum state $\rho$, the purity can be rewritten as

[TABLE]

where $\mathbb{S}=\sum_{\mathbf{s},\mathbf{s}^{\prime}}\ket{\mathbf{s}^{\prime}}\bra{\mathbf{s}}\otimes\ket{\mathbf{s}}\bra{\mathbf{s}^{\prime}}=W_{(2,1)^{\otimes N}}=W_{(2,1)}^{\otimes N}$ is the swap operator acting on two “virtual” copies $\mathcal{H}^{\otimes N}\otimes\mathcal{H}^{\otimes N}$ of the Hilbert space of $N$ qudits. Comparing Eqs. (25) and (26), our goal is thus to find coefficients $O_{\mathbf{s},\mathbf{s}^{\prime}}$ of the operator $O$ such that

[TABLE]

Since $\Phi_{N}^{(2)}(o^{\otimes N})=\left(\Phi_{1}^{(2)}(o)\right)^{\otimes N}$ factorizes for an operator $O=\otimes_{i=1}^{N}o$ [Eq. (24)], it is sufficient to find local operators $o=\sum_{s,s^{\prime}=1}^{d}o_{s,s^{\prime}}\ket{s}\bra{s}\otimes\ket{s^{\prime}}\bra{s^{\prime}}$ which fulfill

[TABLE]

Using Eq. (17), this is equivalent to

[TABLE]

Inserting the ansatz $o=\sum_{s,s^{\prime}=1}^{d}o_{s,s^{\prime}}\ket{s}\bra{s}\otimes\ket{s^{\prime}}\bra{s^{\prime}}$ into Eq. (29), we find the following equations to be satisfied by the coefficients $o_{s,s^{\prime}}$

[TABLE]

These are satisfied by the simple choice

[TABLE]

where $D_{G}[s,s^{\prime}]$ is the Hamming distance of the states $s$ and $s^{\prime}$ of the single qudit, i.e. $D_{G}[s,s]=0$ and $D_{G}[s,s^{\prime}]=1$ if $s\neq s^{\prime}$. On the composite system, we then simply choose

[TABLE]

where $D[\mathbf{s},\mathbf{s}^{\prime}]=\sum_{i=1}^{N}D_{G}[s_{i},s^{\prime}_{i}]$ is the Hamming distance of the states $\mathbf{s}$ and $\mathbf{s}^{\prime}$ of $N$ qudits. This leads directly to Eq. (5).

IV.2 Scaling of statistical errors

In the previous sections, we have shown how to access the purity of an arbitrary quantum state $\rho_{A}$ from the ensemble average over cross correlations of outcome probabilities of randomized measurements. Here, we discuss the statistical errors arising in an experiment due to a finite number $N_{U}$ of unitaries to estimate the ensemble average and a finite number $N_{M}$ of measurements per random unitary to estimate the probabilities $P_{U}(\mathbf{s}_{A})$.

For the protocol utilizing global random unitaries, we found in Ref. Vermersch et al. (2018) analytically a scaling law of the typical statistical error of the estimated purity of a density matrix $\rho_{A}$ in a Hilbert space with dimension $\mathcal{D}_{A}$

[TABLE]

where $c_{1}$ and $c_{2}$ are constants of $\mathcal{O}(1)$ which are largest for pure states. Thus the number of measurements per random unitary required to estimate the purity up to an error of $1/\sqrt{N_{U}}$ scales as $N_{M}\sim\sqrt{\mathcal{D}_{A}}$. The scaling behavior with $N_{U}$ is hereby a direct consequence of the central limit theorem, whereas the scaling with $N_{M}$ is directly related to probability of finding doublons when sampling with replacement from a finite probability distribution (the so-called birthday paradox Blinder (2013)). Note, that we use unbiased estimators to infer the squared outcome probabilities $P_{U}(\mathbf{s}_{A})^{2}$ from a finite number of measurements $N_{M}$ Vermersch et al. (2018). The analytical results are supported by numerical simulations presented in Fig. 8, panels a) and c) where the average statistical error, extracted from 100 numerical experiments, is shown, for a pure and a mixed state. The lines are calculated from the scaling law (31). Clearly, the statistical error of the mixed state is smaller, which is explained by the fact that fluctuations across the unitary ensemble are reduced for mixed states (with vanishing fluctuations for the maximally mixed state $\rho_{\text{max}}\sim\mathbb{1}/\mathcal{D}_{A}$).

For the protocol utilizing local random unitaries, we find, from numerical simulations, see Fig. 8, panels b) and d), for a pure product state of $N_{A}$ qubits, a scaling law

[TABLE]

where $c_{3}=\mathcal{O}(N_{A})$. The scaling with $N_{U}$ follows the behavior expected from the central limit theorem. The number of measurements to estimate the purity up to an error $\sim 1/\sqrt{N_{U}}$ scales as $N_{M}\sim 2^{0.75N_{A}}$, and is thus larger than for the global protocol. However, in contrast to the global protocol, the statistical errors for entangled states of $N_{A}$ qubits are reduced which is explained by the fact that here the reduced density matrices for subsystems are mixed, and thus fluctuations across the unitary ensemble are locally reduced. Similar as in the global protocol, this holds true if $\rho_{A}$ itself is mixed.

We note that on the one hand, the local protocol (ii) is more prone to statistical errors compared to the global protocol (i). On the other hand, we can obtain by restriction from the occupation probabilities $P_{U}(\mathbf{s}_{A})$ of basis states $\ket{\mathbf{s}_{A}}$ of the (sub-)system $A$, the occupation probabilities $P_{U}(\mathbf{s}_{A^{\prime}})$ of basis states $\ket{\mathbf{s}_{A^{\prime}}}$ of an arbitrary subsystem ${A^{\prime}}\subseteq A$. This is not possible for the global protocol, where the applied random unitary $U_{A}$ randomizing the entire Hilbert space $\mathcal{H}_{A}$ of $\rho_{A}$ predefines the (sub-) system $A$ of interest.

IV.3 Measurement of the overlap of quantum states

In this subsection, we discuss a natural extension of the protocol presented in Sec. II.1 to estimate the overlap $\textrm{Tr}\left[\rho_{1}\rho_{2}\right]$ of two distinct quantum states $\rho_{1}$ and $\rho_{2}$ of $N$ qudits in Hilbertspace $\mathcal{H}$ of dimension $\mathcal{D}=d^{N}$. To obtain the overlap $\textrm{Tr}\left[\rho_{1}\rho_{2}\right]$, one applies the experimental sequence described in Sec. II.1 twice with the same set of random unitaries, starting first with the state $\rho_{1}$ and secondly with the state $\rho_{2}$. Repeated for many random unitaries $U$, this provides the set of occupation probabilities $P^{(1)}_{U}(\mathbf{s})=\textrm{Tr}\left[U\rho_{1}U^{\dagger}\ket{\mathbf{s}}\bra{\mathbf{s}}\right]$ and $P^{(2)}_{U}(\mathbf{s})=\textrm{Tr}\left[U\rho_{2}U^{\dagger}\ket{\mathbf{s}}\bra{\mathbf{s}}\right]$ of the computational basis states $\ket{\mathbf{s}}$. From cross correlations, the overlap is estimated. Generalizing Eqs. (4) and (5), one finds, if (i) global random unitaries have been used

[TABLE]

and (ii) for local random unitaries

[TABLE]

These equations follow directly from the proof presented in Sec. IV.1, using that $\textrm{Tr}\left[\rho_{1}\rho_{2}\right]=\textrm{Tr}\left[\mathbb{S}\rho_{1}\otimes\rho_{2}\right]$.

We note that this protocol enables a measurement of the Loschmidt echo $\left|\langle\psi_{0}|e^{iH_{2}t/\hbar}e^{-iH_{1}t/\hbar}|\psi_{0}\rangle\right|^{2}$ Goussev et al. without the necessity of implementing time reversed operations or ancilla degrees of freedom; the protocol is outlined in the following: In a first experiment, $\ket{\psi_{0}}$ is evolved forward in time with Hamiltonian $H_{1}$ and after the application of a random unitary $U$, the probabilities $P^{(1)}_{U}(\mathbf{s})=|\braket{\mathbf{s}}{Ue^{-iH_{1}t/\hbar}}{\psi_{0}}|^{2}$ for the basis states $\ket{\mathbf{s}}$ are measured. This is then repeated with $H_{2}$, and same random unitary $U$, to obtain $P^{(2)}_{U}(\mathbf{s})=|\braket{\mathbf{s}}{Ue^{-iH_{2}t/\hbar}}{\psi_{0}}|^{2}$. The overlap $\left|\langle\psi_{0}|e^{iH_{2}t/\hbar}e^{-iH_{1}t/\hbar}|\psi_{0}\rangle\right|^{2}$ is finally inferred from cross-correlations over many random unitaries, according to Eqs. (33) and (34). Thus, there is no the necessity to implement the time reversed evolution operator $e^{iH_{2}t/\hbar}$ in the experiment.

We further remark that this protocol can be used to check the stability of an experiment against drifts, by measuring the overlap of two quantum states $\rho_{1}$ and $\rho_{2}$ which are prepared in the same way, but at different instances of time.

V Randomized quantum state tomography

In this section, we describe a protocol to perform full quantum state tomography, based on statistical correlation of randomized measurements. For global random unitaries, this protocol was first described in Ref. Ohliger et al. (2013) in the context of atomic Hubbard models. Here, we focus on spin models and extend the protocol to composite systems of many qudits where local unitaries are applied, and investigate in detail the scaling of the required number of measurements to estimate the density matrix up to a fixed statistical error with the Hilbert space dimension (the number of constituents).

We consider a quantum state $\rho_{A}$, which can be a reduced state $\rho_{A}=\textrm{Tr}_{\mathcal{S}\backslash A}\left[\rho\right]$ of a subsystem $A\subseteq\mathcal{S}$, defined in the Hilbert space $\mathcal{H}_{A}=\mathscr{h}^{\otimes N_{A}}$ with total dimension $\mathcal{D}_{A}=d^{N_{A}}$. The use of randomized measurements to perform quantum state tomography is based on the observation that

[TABLE]

where $S$ is the swap operator and the partial trace is taken over the second “copy’. Using the results of Sec. IV.1, this gives immediately rise to a measurement protocol to perform quantum state tomography using randomized measurements. Employing the ensemble average over global random unitaries, randomizing the entire Hilbert space $\mathcal{H}_{A}$, we find

[TABLE]

where the expression in the last line was first obtained in Ref. Ohliger et al. (2013). If we consider instead local random unitaries $U=\bigotimes_{i=1}^{N_{A}}U_{i}$ where the $U_{i}$ are sampled independently from $\text{CUE}(\mathscr{h})$, we find

[TABLE]

The experimental protocol reads thus as follows. Given the quantum state $\rho_{A}$ a random unitary $U_{A}$, being either local or global, is applied. Subsequently, a measurement in the computational basis is performed. This is repeated with the same random unitary, to access the occupation probabilities $P_{U}(\mathbf{s}_{A})$ of the computational basis states $\ket{s}_{A}$. Additionally, the random unitary $U_{A}$ is stored as a matrix in the computational basis. Finally, the ensemble average is performed, as in Eqs. (36) and (37).

We emphasize that, in contrast to the previously discussed protocols, the tomographic reconstruction of $\rho_{A}$ from randomized measurements requires the explicit knowledge of the applied unitaries in a specific basis, i.e. not only their property of being sampled from an appropriate random matrix ensemble. The tomographic reconstruction is thus prone to experimental imperfections which lead to a (random) mismatch between the “applied” random unitary and the one “stored” to reconstruct $\rho_{A}$. In general, such decorrelation will appear as depolarizing noise, i.e. induces a bias towards mixed state. This is not the case for protocols which detect properties of $\rho_{A}$ which are invariant under unitary transformations, such as the purity estimation.

A crucial aspect for any measurement scheme providing tomographic reconstruction of a quantum state is the scaling of the required number of measurements with system size. In the protocol presented here, the accuracy of an estimation $(\rho_{A})_{e}$, quantified by the trace distance $\|(\rho_{A})_{e}-\rho\|_{1}$ to the true state $\rho_{A}$, is determined by statistical errors originating from a finite number of random unitaries $N_{U}$ and a finite number of measurements $N_{M}$ per random unitary. In Fig. 9, the scaling behavior of the trace distance of the estimated density matrix $(\rho_{A})_{e}$ for a pure product state $\rho_{A}$ is shown. For both protocols, we find that the numerical data, obtained with random unitaries sampled directly from the CUE Mezzadri (2007), is well described by a scaling law

[TABLE]

where for the global protocol, $a\approx 1.05$, $b\approx 0.5$, and for the local protocol $a\approx 1.0$, $b\approx 0.75$. From Eq. (38) it follows directly that, for pure product states, the required number of random unitaries $N_{U}$ scales as $\mathcal{D}^{2}_{A}=2^{2N_{A}}$ and is thus comparable to the mininmal number of measurement settings in standard tomography Gross et al. (2010). As shown in Fig. 9, the statistical errors of an estimation of $\rho_{A}$ depend in general on the quantum state, where (absolute) errors smaller for mixed states (see also Brydges et al. (2019); Vermersch et al. (2018)). We note that it was shown in Ref. Ohliger et al. (2013) that the above protocols for randomized state tomography can be combined with compressed sensing Gross et al. (2010) to decrease the number of required measurements.

VI Conclusion

We have introduced statistical correlations of randomized measurements as a new tool to probe complex many-body quantum states. While we have focused on measurements protocol accessing bipartite entanglement of quantum states and their tomographic reconstructions, the underlying tools can be applied more generally. In a recent paper Vermersch et al. , we use statistical correlations to design robust protocols to measure out-of-time ordered correlation functions, without the necessity to implement time-reversed operations or ancilla degrees of freedom. A first experimental demonstration of this protocol was presented in Ref. Nie et al. (2019) in a NMR quantum simulator consisting of four qubits. In the future, the paradigm of statistical correlation could be extended to the measurement of order parameters in (symmetry protected) topological phases Chen et al. (2012); Pollmann and Turner (2012); Haegeman et al. (2012); Shapourian et al. (2017) and the estimation of the entanglement spectrum Li and Haldane (2008); Pichler et al. (2016); Dalmonte et al. (2018). It could further be combined with techniques for randomized benchmarking of quantum computation Emerson et al. (2007); Knill et al. (2008); Wallman et al. (2015).

Our protocols can be applied in state-of-the-art quantum simulators with single site read-out and control. As they rely on statistical correlations of many random measurements, they are particularly suitable for systems with high repetition rates such as trapped-ions, superconducting qubits and Rydberg atoms. In atomic Hubbard models, one can take advantage of the possibility to prepare simultaneously many (independent) copies of the quantum systems to reduce the number of experimental runs.

Furthermore, it would be interesting to extend the protocol based on local unitaries to models with (locally) conserved quantum numbers. For instance, local random unitaries with conserved particle number could be generated in atomic Hubbard models, by isolating pairs of sites and applying to each a series of random quenches.

Acknowledgements.

We thank R. van Bijnen, R. Blatt, T. Brydges, C. Kokail, B. Kraus, L. Sieberer, and J. Yu for discussions. Within the European Union’s Horizon 2020 research and innovation programme, this project has received funding under grant agreement No 817482 and from the European Research Council (ERC) under grant agreement No 741541. It is further supported by the ERC Synergy Grant UQUAM, the SFB FoQuS (FWF Project No. F4016-N23) and QTFLAG – QuantERA. Numerical simulations were realized with QuTiP Johansson et al. (2013).

Appendix A Higher order Rényi entropies

In this section, we discuss the estimation of $k$th-order functionals $\textrm{Tr}\left[\rho^{k}\right]$ ($k\in\mathbb{N},k\geq 2$) from statistical correlations of globally randomized measurements, which are directly connected to $k$-th order Rényi entropies $S^{(k)}(\rho)=1/(1-k)\log_{2}\textrm{Tr}\left[\rho^{k}\right]$.

The experimental sequence is the same as described in section II.1. Global random unitaries $U$ sampled from a unitary $k$-design defined on the entire Hilbert space $\mathcal{H}$ are applied to the quantum state $\rho$ and subsequently occupation probabilities $P_{U}(\mathbf{s})$ are measured. In Ref. Vermersch et al. (2018) we showed that the $k$-th moment $\overline{P_{U}(s)^{k}}$ is related to $\textrm{Tr}\left[\rho^{k}\right]$

[TABLE]

where $\mathcal{D}_{k}=\prod_{i=0}^{k-1}(\mathcal{D}+i)$ and $C_{b_{1},\dots,b_{k}}$ denotes the number of permutations $\sigma\in S_{k}$ with $\textrm{typ}(\sigma)=1^{b_{1}}2^{b_{2}}\ldots k^{b_{k}}$ and is given by Bona (2016)

[TABLE]

This result is easily recovered with the help of the results of section III. We obtain

[TABLE]

where we used that $\sum_{\pi\in S_{k}}C_{\pi\sigma}=1/\prod_{i=0}^{k-1}(\mathcal{D}+i)\equiv 1/\mathcal{D}_{k}$ and $\textrm{Tr}\left[W_{\sigma}\ket{s}\bra{s}^{\otimes k}\right]=1$ for any $\sigma\in S_{k}$. Since any unitary $k$-design is also a unitary $l$-design for $l<k$, we can reconstruct $\textrm{Tr}\left[\rho_{A}^{l}\right]$ for $l<k$ from lower order moments $P_{U}(s)^{l}$. Thus, Eq. (42) can be solved recursively to obtain $\textrm{Tr}\left[\rho^{k}\right]$.

In contrast to the case of global unitaries, the generalization of our method to access $\textrm{Tr}\left[\rho^{k}\right]$ with local unitaries is not straightforward, and will be studied in future work.