Proof. The existence of this form is standard up to permutations, e.g., see [55]. Specialized algorithms exist to compute it [56], although it suffices to recognize that for real skew-symmetric matrices, the real Schur decomposition coincides with this normal form (up to a permutation). ◻

Proof. By triangle inequality, \[\|A^\star - B\| \leq \|A^\star - A\| + \|A - B\| = \|\Sigma(B) - \Sigma(A)\| + \|A - B\|\] where \(\|\Sigma(B) - \Sigma(A)\| = \max_j |\sigma_j(A) - \sigma_j(B)|\). The claim follows from Weyl’s inequality [55]: \[|\sigma_j(A) - \sigma_j(B)| \leq \|A - B\| \quad \forall j. \qedhere\] ◻

Proof. Because \(| \psi \rangle \in \wedge^\eta \mathbb{C}^n\) is already antisymmetrized, we can write \(\Phi_{\mathrm{pas}}(U) | \psi \rangle = U^{\otimes \eta} | \psi \rangle\). Hence \[\mathop{\mathrm{\mathsf{dist}_{tr}}}(\Phi_{\mathrm{pas}}(U) | \psi \rangle, \Phi_{\mathrm{pas}}(V) | \psi \rangle) = \sqrt{1 - |\langle \psi | (U^\dagger V)^{\otimes \eta} | \psi \rangle|^2} \leq \min_{\theta \in \mathbb{R}} \| (U^{\otimes \eta} - e^{i\theta} V^{\otimes \eta}) | \psi \rangle \|.\] By the max–min inequality and telescoping through tensor products, \[\begin{align} \max_{| \psi \rangle \in \wedge^\eta \mathbb{C}^n} \mathop{\mathrm{\mathsf{dist}_{tr}}}(\Phi_{\mathrm{pas}}(U) | \psi \rangle, \Phi_{\mathrm{pas}}(V) | \psi \rangle) &\leq \min_{\theta \in \mathbb{R}} \max_{| \psi \rangle \in \wedge^\eta \mathbb{C}^n} \| (U^{\otimes \eta} - e^{i\theta} V^{\otimes \eta}) | \psi \rangle \| \notag\\ &\leq \min_{\theta \in \mathbb{R}} \|U^{\otimes \eta} - e^{i\theta} V^{\otimes \eta}\|\\ &\leq \min_{\theta \in \mathbb{R}} \eta \|U - e^{i\theta/\eta} V\| = \eta \mathop{\mathrm{\mathsf{dist}_{ph}}}(U, V). \qedhere \end{align}\] ◻

Proof. Constructive algorithms were first given by [50], [51], although the fundamental idea dates back to conventional optics [59]. For completeness, we provide a sketch of the proof. Consider the active case, where without loss of generality we can assume \(W \in \mathrm{SO}(2n)\) (because any \(W' \in \mathrm{O}^{-}(2n)\) can be written as \(W' = XW\) with \(\Phi(X) = \gamma_1\)). We take the QR factorization of \(W = QR\), where \(Q\) is orthogonal and \(R\) is triangular. This can be computed systematically by Givens rotation eliminations in \(\mathcal{O}(n^3)\) time [55]; that is, we determine \(Q = G_L \cdots G_2 G_1\) where \(G_\ell\) are Givens rotations between two adjacent rows and \(L \leq \binom{2n}{2}\). Additionally, since \(W\) is orthogonal, so too must \(R\); this forces \(R\) to be diagonal.

Suppose the Givens rotation \(G\) rotates rows \(q, q+1\) by angle \(\theta\). The desired transformation \[\Phi(G)^\dagger \gamma_p \Phi(G) = \begin{cases} \cos(\theta) \gamma_q - \sin(\theta) \gamma_{q+1} & \text{if } p = q,\\ \sin(\theta) \gamma_q + \cos(\theta) \gamma_{q+1} & \text{if } p = q + 1,\\ \gamma_p & \text{else}. \end{cases}\] can be achieved with \(\Phi(G) = e^{-\frac{\theta}{2} \gamma_p \gamma_{p+1}}\). Meanwhile since \(R\) is diagonal and orthogonal, its diagonal entries can only be \(\pm 1\) so \(\Phi(R)^\dagger \gamma_p \Phi(R) = \pm \gamma_p\). These signs can be implemented by gates of the form \(\gamma_p \gamma_q = i e^{-\frac{\pi}{2} \gamma_p \gamma_q}\), of which there are at most \(n\) (see [60]). Use the homomorphism \(\Phi(W) = \Phi(G_L) \cdots \Phi(G_2) \Phi(G_1) \Phi(R)\) to conclude the gate complexity.

For passive FLOs we can restrict to passive elementary gates because the transformation \[\Phi_{\mathrm{pas}}(G)^\dagger a_p \Phi_{\mathrm{pas}}(G) = \begin{cases} \cos(\theta) a_q - \sin(\theta) a_{q+1} & \text{if } p = q,\\ \sin(\theta) a_q + \cos(\theta) a_{q+1} & \text{if } p = q + 1,\\ a_p & \text{else}. \end{cases}\] can be achieved with \(\Phi_{\mathrm{pas}}(G) = e^{-\theta (a_q^\dagger a_{q+1} - a_{q+1}^\dagger a_q)}\). Also note that in this \(\mathrm{U}(n)\)-representation, the diagonal transformation \(R = \mathop{\mathrm{diag}}(e^{i\alpha_1}, e^{i\alpha_2}, \ldots, e^{i\alpha_n})\) is easily implemented via \(e^{-i \alpha_q a_q^\dagger a_q} a_p e^{i \alpha_q a_q^\dagger a_q} = e^{i\alpha_q \delta_{pq}} a_p\). ◻

Proof. We unpack the notation of [34]. Reducing that result to the \(k = 1\) case yields the expression \[\widehat{\boldsymbol{D}}_{ij} = \langle 1_i | \Phi_{\mathrm{pas}}(\pi_{\boldsymbol{b}} \boldsymbol{V})^\dagger \mathcal{E}_\eta \Phi_{\mathrm{pas}}(\pi_{\boldsymbol{b}} \boldsymbol{V}) | 1_j \rangle,\] where \(\pi_{b} \in \mathbb{Z}^{n \times n}\) is any permutation matrix that maps \(\{i \in [n] : b_i = 1\}\) to \([\eta]\), and the operator \(\mathcal{E}_\eta\) takes the form \[\mathcal{E}_\eta = \sum_{k=1}^n (-1)^{s_k + 1} \binom{n - s_k}{1 - s_k} \binom{n - \eta + s_k}{s_k} | 1_k \rangle\!\langle 1_k |, \quad s_k = \begin{cases} 1, & k \leq \eta,\\ 0, & \text{else}. \end{cases}\] To get the desired expression for \(\widehat{\boldsymbol{D}}\), first recognize that the coefficients in \(\mathcal{E}_\eta\) simplify either to \(n - \eta + 1\) if \(k \leq \eta\), or otherwise to \(-\eta\) for the remaining \(n - \eta\) terms. We can pull \(\pi_{\boldsymbol{b}}\) out of \(\Phi_{\mathrm{pas}}(\cdot)\) and apply it to \(\mathcal{E}_\eta\) to get \[\begin{align} \mathcal{E}_\eta^{(\boldsymbol{b})} &\mathrel{\vcenter{:}}= \Phi_{\mathrm{pas}}(\pi_{\boldsymbol{b}})^\dagger \mathcal{E}_\eta \Phi_{\mathrm{pas}}(\pi_{\boldsymbol{b}})\\ &= (n - \eta + 1) \sum_{k : \boldsymbol{b}_k = 1} | 1_k \rangle\!\langle 1_k | - \eta \sum_{k : \boldsymbol{b}_k = 0} | 1_k \rangle\!\langle 1_k |\\ &= \sum_{k=1}^n E(\boldsymbol{b})_{kk} | 1_k \rangle\!\langle 1_k |. \end{align}\] Finally, use \(\langle 1_k | \Phi_{\mathrm{pas}}(\boldsymbol{V}) | 1_j \rangle = \boldsymbol{V}_{kj}\) to arrive at \[\begin{align} \widehat{\boldsymbol{D}}_{ij} &= \langle 1_i | \Phi_{\mathrm{pas}}(\boldsymbol{V})^\dagger \mathcal{E}_\eta^{(\boldsymbol{b})} \Phi_{\mathrm{pas}}(\boldsymbol{V}) | 1_j \rangle\\ &= \sum_{k=1}^n \langle 1_i | \Phi_{\mathrm{pas}}(\boldsymbol{V})^\dagger | 1_k \rangle E(\boldsymbol{b})_{kk} \langle 1_k | \Phi_{\mathrm{pas}}(\boldsymbol{V}) | 1_j \rangle\\ &= \sum_{k=1}^n [\boldsymbol{V}^\dagger]_{ik} E(\boldsymbol{b})_{kk} \boldsymbol{V}_{kj} = \boldsymbol{V}^\dagger E(\boldsymbol{b}) \boldsymbol{V}. \end{align}\] The claim \(\mathop{\mathrm{\mathbb{E}}}[\widehat{\boldsymbol{D}}] = D\) follows from the fact that classical shadows produces unbiased estimators. ◻

Proof. Since each sample is identically distributed, let us temporarily drop the \(\ell\) subscripts. The bound for \(R\) is given by \[\|\boldsymbol{X}\| \leq \frac{\|\widehat{\boldsymbol{D}}\| + \|D\|}{N} \leq \frac{n + 1}{N},\] which follows from the fact that \(\|D\| \leq 1\) for any state, and from ?? we have \(\|\widehat{\boldsymbol{D}}\| = \|E(\boldsymbol{b})\| = \max\{n+1-\eta, \eta\} \leq n\).

For \(\sigma^2\), we can compute it without needing the third moment of the Haar distribution. Observe that \(\widehat{\boldsymbol{D}}^2 = \boldsymbol{V}^\dagger E(\boldsymbol{b})^2 \boldsymbol{V}\). We can rewrite \(E(\boldsymbol{b})^2\) as \[\label{eq:E40b4195derivation} \begin{align} E(\boldsymbol{b})^2 &= (n+1)^2 \mathop{\mathrm{diag}}(\boldsymbol{b})^2 - 2\eta(n+1) \mathop{\mathrm{diag}}(\boldsymbol{b}) + \eta^2 \mathbb{I}\\ &= (n+1)(n+1-2\eta) \mathop{\mathrm{diag}}(\boldsymbol{b}) + \eta^2 \mathbb{I}\\ &= (n+1-2\eta)\mathopen{}\left[ (n+1) \mathop{\mathrm{diag}}(\boldsymbol{b}) - \eta \mathbb{I}\right]\mathclose{} + \eta(n+1-2\eta) \mathbb{I}+ \eta^2 \mathbb{I}\\ &= (n + 1 - 2\eta) E(\boldsymbol{b}) + \eta(n+1-\eta) \mathbb{I}, \end{align}\tag{4}\] where we used the fact that \(\mathop{\mathrm{diag}}(\boldsymbol{b})^2 = \mathop{\mathrm{diag}}(\boldsymbol{b})\). Then summing over all \(N\) independent samples, we get \[\begin{align} \sum_{\ell=1}^N \mathop{\mathrm{\mathbb{E}}}[\boldsymbol{X}_\ell^2] &= \frac{1}{N^2} \sum_{\ell=1}^N \mathopen{}\left( \mathop{\mathrm{\mathbb{E}}}[\widehat{\boldsymbol{D}}_\ell^2] - \mathop{\mathrm{\mathbb{E}}}[\widehat{\boldsymbol{D}}_\ell]^2 \right)\mathclose{}\\ &= \frac{1}{N} \mathopen{}\left[ (n+1-2\eta) D+ \eta(n+1-\eta) \mathbb{I}- D^2 \right]\mathclose{}. \end{align}\] A crude but sufficient upper bound for the operator norm of this matrix can be found by triangle inequality and dropping all negative terms:⁹ \[\label{eq:variance95op95bound} \mathopen{}\left\| (n+1-2\eta) D+ \eta(n+1-\eta) \mathbb{I}- D^2 \right\|\mathclose{} \leq (n+1) + \eta(n+1) + 1.\tag{5}\] The bound for \(\sigma^2\) follows. ◻

Proof. We adopt the notation from 8 1 and set \(\overline{\boldsymbol{D}} \mathrel{\vcenter{:}}= \frac{1}{N} \sum_{\ell=1}^N \widehat{\boldsymbol{D}}_\ell\). Matrix Bernstein provides the upper bound on the probability in 6 . Invoking 1, which states that \(B \leq \frac{n+1}{N}\) and \(\sigma^2 \leq \frac{(n+1)(\eta+1) + 1}{N}\), we get the bound \(N(\sigma^2 + B\varepsilon/3) \leq 6 n\eta\) for all \(n, \eta \geq 1\) and \(\varepsilon < 1\). Hence \[\begin{align} \Pr\mathopen{}\left( \|\overline{\boldsymbol{D}} - D\| \geq \varepsilon \right)\mathclose{} &\leq 2n \exp\mathopen{}\left( \frac{-N\varepsilon^2}{2N(\sigma^2 + B\varepsilon/3)} \right)\mathclose{}\\ &\leq 2n \exp\mathopen{}\left( -\frac{N\varepsilon^2}{12 n\eta} \right)\mathclose{}. \end{align}\] The claim follows from demanding this be at most \(\delta\). ◻

Proof. Since the Frobenius norm is the entrywise \(2\)-norm, \[\begin{align} \|\Gamma(\psi_1) - \Gamma(\psi_2)\|_F^2 &= 4 \mathopen{}\left\| \begin{pmatrix} \operatorname{Im}[D(\psi_1)] - \operatorname{Im}[D(\psi_2)] & -(\operatorname{Re}[D(\psi_1)] - \operatorname{Re}[D(\psi_2)])\\ \operatorname{Re}[D(\psi_1)] - \operatorname{Re}[D(\psi_2)] & \operatorname{Im}[D(\psi_1)] - \operatorname{Im}[D(\psi_2)] \end{pmatrix} \right\|_F^2\mathclose{}\\ &= 4 \mathopen{}\left( 2\|{\operatorname{Im}[D(\psi_1)] - \operatorname{Im}[D(\psi_2)]}\|_F^2 + 2\|{\operatorname{Re}[D(\psi_1)] - \operatorname{Re}[D(\psi_2)]}\|_F^2 \right)\mathclose{}\\ &= 8 \|D(\psi_1) - D(\psi_2)\|_F^2. \end{align}\] Then use the fact that \(D(\psi_j)\) is a rank-\(\eta\) projector since \(| \psi_j \rangle\) is a Slater determinant. By subadditivity of rank, the rank of the difference \(D(\psi_1) - D(\psi_2)\) is at most \(\min\{2\eta, n\}\). We conclude by chaining the estimate \(\|X\|_F \leq \sqrt{\mathop{\mathrm{rank}}(X)} \|X\|\) with 15. ◻

Proof. The algorithm is described in 1. First we verify its runtime. Recall from 12 that any FLO can be implemented in \(\mathcal{O}(n^2)\) elementary gates, and that this circuit can be determined in \(\mathcal{O}(n^3)\) time. However the time complexity can be accelerated since we are compiling Haar-random FLO circuits, in which case recent work has shown how to draw and construct these circuits in only \(\mathcal{O}(n^2)\) time [61]. Each repetition is therefore dominated by the formation of the estimate \(\overline{\boldsymbol{D}}\), through the cost of matrix multiplication for \(\boldsymbol{V}^\dagger \mathop{\mathrm{diag}}(\boldsymbol{b}) \boldsymbol{V}\) which reduces to multiplying an \(n \times \eta\) matrix with an \(\eta \times n\) matrix. The algorithm then performs one eigendecomposition at the end, an additive \(\mathcal{O}(n^3)\) cost.

Now we show the copy complexity. Let \(\overline{\boldsymbol{D}}\) be the unrounded estimate and \(\boldsymbol{D}^\star\) the rounded RDM, corresponding to a Slater determinant \(| \widehat{\boldsymbol{\psi}} \rangle\). By 6 2, we have \[\mathop{\mathrm{\mathsf{dist}_{tr}}}(| \widehat{\boldsymbol{\psi}} \rangle, | \psi \rangle) \leq \sqrt{\eta} \|\boldsymbol{D}^\star - D\| \leq 2 \sqrt{\eta} \|\overline{\boldsymbol{D}} - D\|.\] Therefore we want the spectral error on \(\overline{\boldsymbol{D}}\) to be at most \(\frac{\varepsilon}{2\sqrt{\eta}}\). This occurs with probability \(\geq 1 - \delta\) as long as we take \(N = \mathopen{}\left\lceil \frac{48 n\eta^2 \log(2n/\delta)}{\varepsilon^2} \right\rceil\mathclose{}\), per 14. ◻

Proof. Use 6 to get \[\| | \widehat{\boldsymbol{u}} \rangle\!\langle \widehat{\boldsymbol{u}} | - | u \rangle\!\langle u | \| \leq 2 \| \overline{\boldsymbol{D}} - | u \rangle\!\langle u | \|,\] hence \[\frac{1}{\sqrt{2}} \| | \widehat{\boldsymbol{u}} \rangle\!\langle \widehat{\boldsymbol{u}} | - | u \rangle\!\langle u | \|_F \leq \| | \widehat{\boldsymbol{u}} \rangle\!\langle \widehat{\boldsymbol{u}} | - | u \rangle\!\langle u | \| \leq 2 \| \overline{\boldsymbol{D}} - | u \rangle\!\langle u | \|.\] Set \(t = \varepsilon/2\) in 18 to arrive at the claim. ◻

Proof. Because the algorithm is the same as 17, the time complexity is also the same. The only difference is that we take \(N = \mathopen{}\left\lceil \frac{384(11 n + 5 \log(2/\delta))}{\varepsilon^2} \right\rceil\mathclose{}\) in 1. To establish correctness, chase the proof of 2 to get \[\mathop{\mathrm{\mathsf{dist}_{tr}}}(| \psi(\widehat{\boldsymbol{u}}) \rangle, | \psi(u) \rangle) \leq \frac{1}{\sqrt{2}} \| | \widehat{\boldsymbol{u}} \rangle\!\langle \widehat{\boldsymbol{u}} | - | u \rangle\!\langle u | \|_F = \sqrt{1 - |\langle \widehat{\boldsymbol{u}} | u \rangle|^2}\] along with 3. Because the uniform POVM is symmetric on \(\mathbb{S}^{n-1}\), and all postprocessing is symmetric with respect to the subspace orthogonal to \(| u \rangle\), the error \(| \boldsymbol{w} \rangle\) is Haar random in that subspace. ◻

Proof. This is essentially proven in [28]. We only address some minor details¹² relevant to our setting. The original statement of applies directly to the passive FLO case over the projective metric. But the distance bounds they use also hold for the non-projective metric, with slightly smaller constants: \[\label{eq:frac95power95nonproj} \|U^{1/p} - V^{1/p}\| \leq \frac{\pi}{p} \|U - V\|\tag{8}\] for any \(p \geq 1\) and \(U = e^X, V = e^Y\) such that \(\|X\|, \|Y\| < \frac{1}{\pi}\) [28]. The rest of their argument can then be used without modification.

The statement also applies to active FLOs. It clearly holds over \(\mathrm{SO}(2n)\) because this is a connected Lie subgroup of \(\mathrm{U}(2n)\). To extend to the full orthogonal group, recall that the coset \(\mathrm{O}^{-}(2n)\) is equal to \(X \cdot \mathrm{SO}(2n)\) for any reflection \(X\). Thus for any two \(U, V \in \mathrm{O}^{-}(2n)\) there exist \(U', V' \in \mathrm{SO}(2n)\) such that \(U = XU'\), \(V = XV'\), and so the distance bounds apply by unitary invariance. Note that in the case that, say, \(U \in \mathrm{SO}(2n)\) but \(V \in \mathrm{O}^{-}(2n)\), we always have \(\min_{s\in\{\pm 1\}} \|U - sV\| = \|U - V\| = 2\); taking \(\varepsilon_0 \leq \frac{1}{10}\) is more than enough to avoid this happening.

Finally we remark on some constants. Because we are not concerned with Item (b) from [28], one can choose a larger constant \(\varepsilon_0 = \frac{1}{50}\) and exponent base for \(\delta_t = \frac{\delta}{2^{T+1-t}}\) to achieve the desired final error and success probability. In fact, for the non-projective metric we can take any \(\varepsilon_0 < \frac{1}{3\pi}\) because in that case, (1) we only need to be within a ball of radius \(r = \frac{1}{\pi}\) to improve the precision per iteration (guaranteed if \(3\varepsilon_0 < r)\) and (2) the tighter bound of 8 implies that the final error is \(\frac{\pi}{p_T} \varepsilon_0 \leq \frac{1}{p_T} \leq \varepsilon\) if \(\varepsilon_0 \leq \frac{1}{\pi}\) (already guaranteed by choosing \(\varepsilon_0 < \frac{1}{3\pi}\)). ◻

Proof. The proof is contained in the second half of [28]. It is applicable to our setting because we have an equivalent form of state tomography using \(\mathrm{U}(n)\)-shadows (19). We track down constants by a alternate version of their argument, which we defer to 9. ◻

Proof. Observe that \(\langle \widetilde{\Psi}(\boldsymbol{\theta}) | X | \widetilde{\Psi}(\boldsymbol{\theta}) \rangle = \langle \Psi(\boldsymbol{\theta}) | \Phi_{\mathrm{pas}}(\boldsymbol{W})^\dagger X \Phi_{\mathrm{pas}}(\boldsymbol{W}) | \Psi(\boldsymbol{\theta}) \rangle\). Expand: \[\Phi_{\mathrm{pas}}(\boldsymbol{W})^\dagger a_1^\dagger a_{n+1}^\dagger \Phi_{\mathrm{pas}}(\boldsymbol{W}) = \sum_{j,k=1}^{n+1} [\boldsymbol{W} \oplus 1]_{1j}^* [\boldsymbol{W} \oplus 1]_{n+1,k}^* a_j^\dagger a_k^\dagger\] and \[\begin{align} \langle \Psi(\boldsymbol{\theta}) | a_j^\dagger a_k^\dagger | \Psi(\boldsymbol{\theta}) \rangle &= \frac{1}{2} \mathopen{}\left( \langle 0^{n+1} | a_j^\dagger a_k^\dagger | 0^{n+1} \rangle + e^{i\boldsymbol{\theta}} \langle 0^{n+1} | a_j^\dagger a_k^\dagger | 1_1 1_{n+1} \rangle \right.\mathclose{}\\ &\hphantom{=~} \mathopen{}\left. + \, e^{-i\boldsymbol{\theta}} \langle 1_1 1_{n+1} | a_j^\dagger a_k^\dagger | 0^{n+1} \rangle + \langle 1_1 1_{n+1} | a_j^\dagger a_k^\dagger | 1_1 1_{n+1} \rangle \right)\mathclose{} \notag\\ &= \frac{1}{2} e^{-i\boldsymbol{\theta}} (\delta_{1j} \delta_{n+1,k} - \delta_{1k} \delta_{n+1,j}), \end{align}\] which are the entries of the \((n+1) \times (n+1)\) matrix \[A(\boldsymbol{\theta}) \mathrel{\vcenter{:}}= \frac{1}{2} e^{-i\boldsymbol{\theta}} \begin{pmatrix} 0_{n \times n} & | 1 \rangle\\ -\langle 1 | & 0 \end{pmatrix}.\] Thus \[\begin{align} \langle \Psi(\boldsymbol{\theta}) | \Phi_{\mathrm{pas}}(\boldsymbol{W})^\dagger a_1^\dagger a_{n+1}^\dagger \Phi_{\mathrm{pas}}(\boldsymbol{W}) | \Psi(\boldsymbol{\theta}) \rangle &= [(\boldsymbol{W}^* \oplus 1) A(\boldsymbol{\theta}) (\boldsymbol{W}^\dagger \oplus 1)]_{1,n+1}\\ &= \frac{1}{2} e^{-i\boldsymbol{\theta}} \mathopen{}\left[ \begin{pmatrix} 0_{n \times n} & \boldsymbol{W}^* | 1 \rangle\\ -\langle 1 | \boldsymbol{W}^\dagger & 0 \end{pmatrix} \right]_{1,n+1}\mathclose{}\\ &= \frac{1}{2} e^{-i\boldsymbol{\theta}} \boldsymbol{W}_{11}^*. \end{align}\] If we denote this by \(\boldsymbol{w}\), then simplifying \(\langle X \rangle = \boldsymbol{w} + \boldsymbol{w}^*\) and \(\langle Y \rangle = i(\boldsymbol{w} - \boldsymbol{w}^*)\) yields the claim. ◻

Proof. First we bound the statistical error. Although we can in principle evaluate the variance of \(X\) and \(Y\), it suffices to simply bound their ranges. The operators \(X\) and \(Y\) are diagonalized by Bogoliubov transformations: \[\label{eq:quadrature95diagonalization} b_1(\varphi) \mathrel{\vcenter{:}}= \frac{a_1 + e^{i\varphi} a_{n+1}^\dagger}{\sqrt{2}}, \quad b_2(\varphi) \mathrel{\vcenter{:}}= \frac{a_{n+1} - e^{i\varphi} a_1^\dagger}{\sqrt{2}}, \quad \text{and } d_j(\varphi) \mathrel{\vcenter{:}}= b_j^\dagger(\varphi) b_j(\varphi),\tag{11}\] so that \[\begin{align} X &= d_1(0) + d_2(0) - \mathbb{I},\\ Y &= \textstyle{d_1(\frac{\pi}{2}) + d_2(\frac{\pi}{2}) - \mathbb{I}}. \end{align}\] Thus the spectra of \(X\) and \(Y\) are \(\{-1, 0, 1\}\). Measuring in the eigenbasis of these quadratures yields random variables with magnitudes bounded by \(1\). Write \(\boldsymbol{W}_{11} = \boldsymbol{r} e^{i\boldsymbol{\xi}}\) per 24 and let \(t \geq 0\). By Hoeffding’s inequality (7), we can compute estimates \(\widehat{\boldsymbol{c}}\) and \(\widehat{\boldsymbol{s}}\) such that \[|\widehat{\boldsymbol{c}} - \boldsymbol{r} \cos(\boldsymbol{\theta} + \boldsymbol{\xi})| \leq t \quad \text{and} \quad |\widehat{\boldsymbol{s}} - \boldsymbol{r} \sin(\boldsymbol{\theta} + \boldsymbol{\xi}))| \leq t\] except with probability \(2p'\), provided that \(N \geq \frac{2 \log(2/p')}{t^2}\).

Next we combine this with the systematic error due to \(\boldsymbol{W}\). By triangle inequality, the necessary cosine/sine bounds for applying 25 are: \[\begin{align} |\widehat{\boldsymbol{c}} - \cos(\boldsymbol{\theta} + \boldsymbol{\xi})| &\leq |\widehat{\boldsymbol{c}} - \boldsymbol{r} \cos(\boldsymbol{\theta} + \boldsymbol{\xi})| + |\boldsymbol{r} \cos(\boldsymbol{\theta} + \boldsymbol{\xi}) - \cos(\boldsymbol{\theta} + \boldsymbol{\xi})| \notag\\ &\leq t + |\boldsymbol{r} - 1|\\ &\leq t + \|\boldsymbol{W} - \mathbb{I}\| \notag\\ &\leq t + \varepsilon , \end{align}\] and similarly for \(|\widehat{\boldsymbol{s}} - \sin(\boldsymbol{\theta} + \boldsymbol{\xi})|\). Use the fact that, for \(0 \leq \tau < \frac{\pi}{2}\), the condition \(\delta \leq \frac{\sin(\tau)}{\sqrt{2}}\) always holds whenever \(\delta \leq \frac{\sqrt{2}}{\pi} \tau\). Taking \(\delta = t + \varepsilon\), one choice of constants that satisfies this linear inequality is \(t = (\sqrt{2} - 1) \varepsilon\) and \(\tau = \pi \varepsilon\). This implies that \(\widehat{\boldsymbol{\theta}}\) satisfies \(|\widehat{\boldsymbol{\theta}} - (\boldsymbol{\theta} + \boldsymbol{\xi})| \leq \pi \varepsilon\), so that \[\begin{align} |e^{i\widehat{\boldsymbol{\theta}}} - e^{i\boldsymbol{\theta}}| &\leq |e^{i(\widehat{\boldsymbol{\theta}} - (\boldsymbol{\theta} + \boldsymbol{\xi}))} - 1| + |e^{i\boldsymbol{\xi}} - 1|\\ &\leq |\widehat{\boldsymbol{\theta}} - (\boldsymbol{\theta} + \boldsymbol{\xi})| + |\boldsymbol{W}_{11} - 1| + |\boldsymbol{r} - 1|\\ &\leq (\pi + 2) \varepsilon, \end{align}\] where we used the inequality \(|e^{ix} - 1| = 2|{\sin(x/2)}| \leq |x|\). ◻

Proof. Run the protocol of 26 with \(p' = \frac{p}{2}\), where copies of \(| \widetilde{\Psi}(\boldsymbol{\theta}) \rangle\) are prepared via \[| \widetilde{\Psi}(\boldsymbol{\theta}) \rangle = \Phi_{\mathrm{pas}}(U) \Phi_{\mathrm{pas}}(\boldsymbol{U}^\star)^\dagger e^{\frac{\pi}{4} (a_1^\dagger a_{n+1}^\dagger - a_{n+1} a_1)} | 0^{n+1} \rangle.\] From the phase estimate \(\widehat{\boldsymbol{\theta}}\) construct \(\boldsymbol{U}^{\sharp} \mathrel{\vcenter{:}}= e^{i \widehat{\boldsymbol{\theta}}} \boldsymbol{U}^\star\), which obeys \[\begin{align} \|\boldsymbol{U}^{\sharp} - U\| &= \|e^{i (\widehat{\boldsymbol{\theta}} - \boldsymbol{\theta})} \boldsymbol{W}^\dagger - \mathbb{I}\| \notag\\ &\leq \|e^{i (\widehat{\boldsymbol{\theta}} - \boldsymbol{\theta})} \boldsymbol{W}^\dagger - \boldsymbol{W}^\dagger\| + \|\boldsymbol{W}^\dagger - \mathbb{I}\|\\ &\leq (\pi + 2) \varepsilon + \varepsilon. \qedhere \end{align}\] ◻

Proof. \(\boldsymbol{V}\) and \(\boldsymbol{G}\) are estimated with error \(\frac{\varepsilon}{175}\) except with probability at most \(\frac{\delta}{4}\) each. By 23, the phase-corrected unitary \(\boldsymbol{U}^\star\) is \(\frac{\varepsilon}{7}\)-close to \(U\) in \(\mathop{\mathrm{\mathsf{dist}_{ph}}}\) error, except with probability \(\frac{\delta}{2}\). Running the phase estimation step, we can conclude the final error bound via 2. ◻

Proof. As 6 indicates, \(\mathsf{PassiveTomo}(\Phi(Q), \frac{1}{10}, \delta)\) makes \(\mathcal{O}(n^2 \log(1/\delta))\) queries. Hence by 21 the bootstrapped process with error \(\varepsilon/n\) makes a total of \(\mathcal{O}(n^3 \log(1/\delta) / \varepsilon)\) queries. According to 10, \(\|\boldsymbol{U}^\sharp - U\| \leq \varepsilon/n\) implies that the diamond distance between the FLOs is at most \(\varepsilon\).

The gate count per experiment is as follows. We use the fact that any FLO can be implemented in \(\mathcal{O}(n^2)\) quantum gates (12). Let \(\boldsymbol{V}_t\) be the current estimate at each iteration \(t = 0, 1, \ldots, T = \lceil \log_2(n/\varepsilon) \rceil\) and \(p_t = 2^t\) (see 2). Synthesizing the unitary \((\Phi_{\mathrm{pas}}(U) \Phi_{\mathrm{pas}}(\boldsymbol{V}_t^\dagger))^{p_t}\) requires \(\mathcal{O}(p_t n^2)\) gates, which is the dominant gate complexity. Thus we use at most \(\mathcal{O}(p_T n^2) = \mathcal{O}(n^3/\varepsilon)\) gates per experiment.

We conclude with the classical cost. The circuit for \(\Phi_{\mathrm{pas}}(\boldsymbol{V}_t^\dagger)\) can be determined in \(\mathcal{O}(n^3)\) time (this only needs to be calculated once per iteration). Within each call to the base tomography, we have:

1. \(\mathsf{PhaselessTomo}\) calls \(\mathsf{SlaterTomo}\) \(n\) times, plus a final SVD, for a total of \(\mathcal{O}(n(n^2 N_t + n^3) + n^3) = \mathcal{O}(n^4 + n^3 \log(n/\delta_t))\) operations (19);

2. Determining the circuits for \(\Phi_{\mathrm{pas}}(F^\dagger)\) and \(\Phi_{\mathrm{pas}}(\boldsymbol{U}^\star)^\dagger\) costs \(\mathcal{O}(n^3)\) operations (12);

3. \(\mathsf{ColumnPhases}\) costs \(\mathcal{O}(n^2)\) operations since we only perform elementwise matrix operations (4);

4. \(\mathsf{PhaseEst}\) costs \(\mathcal{O}(N_{{\mathrm{ph}},t}) = \mathcal{O}(\log(1/\delta_t))\) operations (5).

Summing over \(t\) from \(0\) to \(T\) with \(\delta_t = \frac{\delta}{2^{T+1-t}}\) yields a total classical computational complexity of \[\begin{align} \mathcal{O}\mathopen{}\left( \sum_{t=0}^T \mathopen{}\left( n^4 + n^3 \log(n/\delta_t) \right)\mathclose{} \right)\mathclose{} = \mathcal{O}\mathopen{}\left( n^4 T + n^3 (T^2 + \log(n/\delta) \right)\mathclose{} \end{align}\] where \(T = \mathcal{O}(\log(n/\varepsilon))\). ◻

Proof. For convenience, we diagonalize the covariance matrices: \[\label{eq:covmat95diagonalization} \Gamma_1 = U_1 D U_1^\dagger, \quad \text{where } U_1 = Q_1 \Omega^\mathrm{T}, \;\Omega \mathrel{\vcenter{:}}= \frac{1}{\sqrt{2}} \begin{pmatrix} \mathbb{I}& i\mathbb{I}\\ \mathbb{I}& -i\mathbb{I} \end{pmatrix}, \text{ and } D = \begin{pmatrix} i\mathbb{I}& 0\\ 0 & -i\mathbb{I} \end{pmatrix},\tag{12}\] and similarly \(\Gamma_2 = U_2 D U_2^\dagger\) with \(U_2 = Q_2 \Omega^\mathrm{T}\). By unitary invariance of the operator norm, we have \[\|\Gamma_1 - \Gamma_2\| = \|U_1^\dagger (\Gamma_1 - \Gamma_2) U_2\| = \|DW - WD\|\] where \(W \mathrel{\vcenter{:}}= U_1^\dagger U_2\). Expressing \[W = \begin{pmatrix} W_{11} & W_{12}\\ W_{21} & W_{22} \end{pmatrix},\] we get \[\label{eq:DW95commutator95bound} \begin{align} \|\Gamma_1 - \Gamma_2\| &= \|DW - WD\|\\ &= \mathopen{}\left\| \begin{pmatrix} W_{11} & W_{12}\\ -W_{21} & -W_{22} \end{pmatrix} - \begin{pmatrix} W_{11} & -W_{12}\\ W_{21} & -W_{22} \end{pmatrix} \right\|\mathclose{}\\ &= \mathopen{}\left\| \begin{pmatrix} 0 & 2 W_{12}\\ -2 W_{21} & 0 \end{pmatrix} \right\|\mathclose{}\\ &= 2 \max\{\|W_{12}\|, \|W_{21}\|\}. \end{align}\tag{13}\] In fact, from the structure of \(W = \Omega^* (Q_1^\mathrm{T}Q_2) \Omega^\mathrm{T}\), we have \(W_{21} = W_{12}^*\), hence \(\|\Gamma- \Gamma'\| = 2 \|W_{12}\|\). This can be seen by direct calculation: \[\label{eq:W95symmetry} \begin{align} W &= \Omega^* \underbrace{\begin{pmatrix} A & B\\ C & D \end{pmatrix}}_{Q_1^\mathrm{T}Q_2} \Omega^\mathrm{T}\\ &= \frac{1}{2} \begin{pmatrix} (A + D) + i(B - C) & (A - D) - i(B + C)\\ (A - D) + i(B + C) & (A + D) - i(B - C) \end{pmatrix}\\ &= \begin{pmatrix} W_{11} & W_{12}\\ W_{12}^* & W_{11}^* \end{pmatrix}. \end{align}\tag{14}\] Furthermore, the singular values of the blocks \(W_{11}, W_{12}\) are tightly related. This can be seen from the CSD of \(W \in \mathrm{U}(2n)\), which states that there exist unitaries \(X_1, X_2, Y_1, Y_2 \in \mathrm{U}(n)\) such that [55]: \[\label{eq:csd95W} \begin{pmatrix} X_1 & 0\\ 0 & X_2 \end{pmatrix} \begin{pmatrix} W_{11} & W_{12}\\ W_{21} & W_{22} \end{pmatrix} \begin{pmatrix} Y_1^\dagger & 0\\ 0 & Y_2^\dagger \end{pmatrix} = \begin{pmatrix} C & S\\ -S & C \end{pmatrix}.\tag{15}\] Here, \(C\) is a diagonal matrix containing the singular values of \(W_{11}\) in non-increasing order and \(S = \sqrt{\mathbb{I}- C^2}\) contains the singular values of \(W_{12}\) (in reverse order). Note that the singular values of the blocks all lie within \([0, 1]\) due to unitarity.

Next we compare this to \(\|Q_1 - Q_2 R\|\) for some \(R \in \mathrm{O}(2n) \cap \mathrm{Sp}(2n, \mathbb{R})\). It is useful to express \[\begin{align} \|Q_1 - Q_2 R\| &= \|\mathbb{I}- Q_1^\mathrm{T}Q_2 R\|\\ &= \|\mathbb{I}- W \Omega^* R \Omega^\mathrm{T}\|\\ &= \mathopen{}\left\|\mathbb{I}- W \begin{pmatrix} V^* & 0\\ 0 & V \end{pmatrix} \right\|\mathclose{}, \end{align}\] where the final line comes from parametrizing \[R = \begin{pmatrix} \operatorname{Re}V & -\operatorname{Im}V\\ \operatorname{Im}V & \operatorname{Re}V \end{pmatrix}\] for some \(V \in \mathrm{U}(n)\) and conjugating it with \(\Omega^* (\cdot) \Omega^\mathrm{T}\) (the calculation is analogous to 14 ). Take the left polar decomposition of \(W_{11}\): \[W_{11} = HZ, \text{ where } H \mathrel{\vcenter{:}}= \sqrt{W_{11} W_{11}^\dagger}, \;Z \in \mathrm{U}(n).\] The minimum over all \(V \in \mathrm{U}(n)\) is at most the value at any particular point, so we can get an upper bound by evaluating the norm at \(V = Z^\mathrm{T}\): \[\begin{align} \min_{R \in \mathrm{O}(2n) \cap \mathrm{Sp}(2n, \mathbb{R})} \|Q - Q'R\| &= \min_{V \in \mathrm{U}(n)} \mathopen{}\left\|\mathbb{I}- W \begin{pmatrix} V^* & 0\\ 0 & V \end{pmatrix} \right\|\mathclose{}\\ &\leq \mathopen{}\left\|\mathbb{I}- W \begin{pmatrix} Z^\dagger & 0\\ 0 & Z^\mathrm{T} \end{pmatrix} \right\|\mathclose{}\\ &= \mathopen{}\left\| \begin{pmatrix} \mathbb{I}- W_{11} Z^\dagger & W_{12} Z^\mathrm{T}\\ W_{12}^* Z^\dagger & \mathbb{I}- W_{11}^* Z^\mathrm{T} \end{pmatrix} \right\|\mathclose{}\\ &= \mathopen{}\left\| \begin{pmatrix} \mathbb{I}- H & W_{12} Z^\mathrm{T}\\ W_{12}^* Z^\dagger & \mathbb{I}- H^* \end{pmatrix} \right\|\mathclose{}. \end{align}\] We proceed by triangle inequality, splitting the matrix into its diagonal and off-diagonal blocks. For the diagonal blocks, since \(H\) is Hermitian, \(H^* = H^\mathrm{T}\) has the same spectrum. Thus it suffices to consider the norm of \(\mathbb{I}- H\). By construction, the eigenvalues of \(H\) are the singular values of \(W_{11}\), so from the CSD we have that \(0 \preceq H \preceq \mathbb{I}\) and hence \[\begin{align} \|\mathbb{I}- H\| &= 1 - \sigma_{\min}(H)\\ &= 1 - \langle n | C | n \rangle\\ &= 1 - \sqrt{1 - \langle 1 | S | 1 \rangle^2}\\ &\leq \langle 1 | S | 1 \rangle^2 = \|W_{12}\|^2. \end{align}\] Meanwhile, the off-diagonal blocks clearly have operator norm \(\|W_{12}\|\). Altogether, combining these facts with 13 we get \[\begin{align} \min_{R \in \mathrm{O}(2n) \cap \mathrm{Sp}(2n, \mathbb{R})} \|Q_1 - Q_2 R\| &\leq \mathopen{}\left\| \begin{pmatrix} \mathbb{I}- H & W_{12} Z\\ W_{12}^* Z^* & \mathbb{I}- H^* \end{pmatrix} \right\|\mathclose{}\\ &\leq \mathopen{}\left\| \begin{pmatrix} \mathbb{I}- H & 0\\ 0 & \mathbb{I}- H^* \end{pmatrix} \right\|\mathclose{} + \mathopen{}\left\| \begin{pmatrix} 0 & W_{12} Z\\ W_{12}^* Z^* & 0 \end{pmatrix} \right\|\mathclose{}\\ &\leq \|W_{12}\|^2 + \|W_{12}\|\\ &\leq 2 \|W_{12}\|\\ &= \|\Gamma_1 - \Gamma_2\|, \end{align}\] which is what we had set out to prove. ◻

Proof. The algorithm is \(\mathsf{GaussianTomo}(\Phi(Q) | 0^{n} \rangle, N_{{\mathrm{act}}})\) from 8 except we output the orthogonal matrix \(\boldsymbol{W}\) directly. The number of copies follows from 41 6, and 4 converts the covariance matrix error to orthogonal matrix error. ◻

Proof. Recall the relation between annihilation/creation operators and Majorana operators: \[a_j = \frac{\gamma_j - i\gamma_{j+n}}{2}, \quad a_j^\dagger = \frac{\gamma_j + i\gamma_{j+n}}{2}.\] Also recall that the covariance matrix entries are \[\Gamma_{j,k} = -i \mathop{\mathrm{tr}}(\gamma_j \gamma_k \rho) + i\delta_{jk}.\] Thus we can write the \(1\)-RDM as \[\begin{align} D_{jk} &= \mathop{\mathrm{tr}}(a_j^\dagger a_k \rho)\\ &= \frac{1}{4} \mathopen{}\left( \mathop{\mathrm{tr}}(\gamma_j \gamma_k \rho) + \mathop{\mathrm{tr}}(\gamma_{j+n} \gamma_{k+n} \rho) + i\mathop{\mathrm{tr}}(\gamma_{j+n} \gamma_k \rho) - i\mathop{\mathrm{tr}}(\gamma_j \gamma_{k+n} \rho) \right)\mathclose{}\\ &= \frac{1}{4} \mathopen{}\left( 2\delta_{jk} + i\Gamma_{j,k} + i\Gamma_{j+n,k+n} - \Gamma_{j+n,k} + \Gamma_{j,k+n} \right)\mathclose{}. \end{align}\] Equivalently in block matrix form, if we write \[\Gamma= \begin{pmatrix} \Gamma_{11} & \Gamma_{12}\\ \Gamma_{21} & \Gamma_{22} \end{pmatrix}\] with each block an \(n \times n\) real matrix, then \[D= \frac{\mathbb{I}+ \frac{1}{2}(\Gamma_{12} - \Gamma_{21}) + \frac{i}{2}(\Gamma_{11} + \Gamma_{22})}{2} = \frac{\mathbb{I}+ T(\Gamma)}{2}.\]

To compute the induced norm of \(T\), observe that we can write \[T(X) = i P X P^\dagger, \quad \text{where } P \mathrel{\vcenter{:}}= \frac{1}{\sqrt{2}} \begin{pmatrix} \mathbb{I}& i\mathbb{I} \end{pmatrix} \in \mathbb{C}^{n \times 2n}.\] Since \(PP^\dagger = \mathbb{I}\), \(\|P\| = 1\). Hence \(\|T(X)\| \leq \|P\|^2 \|X\| = \|X\|\). For the reverse direction, consider \(X = \mathop{\mathrm{diag}}(\mathbb{I}, \mathbb{I})\). Then \(T(X) = i\mathbb{I}\) attains \(\|T(X)\| = \|X\|\). ◻

Proof. Let \(\Gamma(\sigma)\) be the covariance matrix of \(\sigma\). By 29, we have \(D(\rho_1) = \frac{1}{2}(\mathbb{I}+ T(Q_1 \Gamma(\sigma) Q_1^\mathrm{T}))\) and \(D(\rho_2) = \frac{1}{2}(\mathbb{I}+ T(Q_2 \Gamma(\sigma) Q_2^\mathrm{T}))\). Set \(W \mathrel{\vcenter{:}}= Q_1^\mathrm{T}Q_2\). Using linearity of \(T\) and the fact that \(\|T\|_{\infty\to\infty} = 1\), we get \[\begin{align} \|D(\rho) - D(\rho')\| &= \frac{1}{2} \|T(Q_1 \Gamma(\sigma) Q_1^\mathrm{T}- Q_2 \Gamma(\sigma) Q_2^\mathrm{T})\| \notag\\ &\leq \frac{1}{2} \|T\|_{\infty\to\infty} \|Q_1 \Gamma(\sigma) Q_1^\mathrm{T}- Q_2 \Gamma(\sigma) Q_2^\mathrm{T}\| \notag\\ &= \frac{1}{2} \|\Gamma(\sigma) - W \Gamma(\sigma) W^\mathrm{T}\|\\ &= \frac{1}{2} \|\Gamma(\sigma)(\mathbb{I}- W^\mathrm{T}) + (\mathbb{I}- W) \Gamma(\sigma) W^\mathrm{T}\| \notag\\ &\leq \|\mathbb{I}- W\| = \|Q_1 - Q_2\|. \qedhere \end{align}\] ◻

Proof. Apply 5 with \(\sigma = \Phi_{\mathrm{pas}}(U) | 1_j \rangle\!\langle 1_j | \Phi_{\mathrm{pas}}(U)^\dagger\), \(Q_1 = Z\), and \(Q_2 = \mathbb{I}\). ◻

Proof. We condition on the event in 18 , which occurs for all \(j \in [n]\) except with probability at most \(n \eta_{{\mathrm{pas}}}\) by a union bound. Use 6 to assert that \[\begin{align} \|| \widehat{\boldsymbol{u}}_j \rangle\!\langle \widehat{\boldsymbol{u}}_j | - | u_j \rangle\!\langle u_j |\| &\leq 2\|\overline{\boldsymbol{D}}_j - | u_j \rangle\!\langle u_j |\|\\ &\leq 2(\|\overline{\boldsymbol{D}}_j - D_j\| + \|D_j - | u_j \rangle\!\langle u_j |\|)\\ &\leq 2(\delta_{{\mathrm{pas}}} + \delta_{{\mathrm{act}}}). \end{align}\] It is straightforward to convert this to a Euclidean distance between the vectors: \[\begin{align} \min_{\theta \in \mathbb{R}} \| | \widehat{\boldsymbol{u}}_j \rangle - e^{i\theta} | u_j \rangle \| &= \sqrt{2 - 2 |\langle \widehat{\boldsymbol{u}}_j | u_j \rangle|}\\ &\leq \sqrt{2(1 - |\langle \widehat{\boldsymbol{u}}_j | u_j \rangle|^2)}\\ &= \frac{\sqrt{2}}{2} \|| \widehat{\boldsymbol{u}}_j \rangle\!\langle \widehat{\boldsymbol{u}}_j | - | u_j \rangle\!\langle u_j |\|_1\\ &\leq \sqrt{2} \|| \widehat{\boldsymbol{u}}_j \rangle\!\langle \widehat{\boldsymbol{u}}_j | - | u_j \rangle\!\langle u_j |\|\\ &\leq 2\sqrt{2} (\delta_{{\mathrm{pas}}} + \delta_{{\mathrm{act}}}). \end{align}\] Hence there exists some diagonal real matrix \(\Theta\) such that \[\label{eq:unitary95fro95bound} \|\widehat{\boldsymbol{U}} - U e^{i\Theta}\| \leq \|\widehat{\boldsymbol{U}} - Ue^{i\Theta}\|_F \leq 2\sqrt{2n} (\delta_{{\mathrm{pas}}} + \delta_{{\mathrm{act}}}),\tag{20}\] from which 19 follows by another application of 6. ◻

Proof. The proof is identical to 23, using 6 as the base “up to column phases” algorithm. ◻

Proof. As alluded to above, this follows from a more general property of classical shadows. Recall that the \(\mathrm{U}(n)\)-shadow channel \(\mathcal{M}\) is an average over two basic operations: rotation by the group \(\Phi_{\mathrm{pas}}(\mathrm{U}(n))\) and projection into the number basis \(\{0, 1\}^n\). Both are block-diagonal in the number basis, which implies that \(\mathcal{M}\) admits the orthogonal decomposition \[\mathcal{M} = \bigoplus_{k=0}^n \mathcal{M}_k.\] In the case that \(\rho\) lies entirely in one particular number sector, say \(k\), then all postmeasurement states \(\boldsymbol{\sigma} = \Phi_{\mathrm{pas}}(\boldsymbol{V})^\dagger | \boldsymbol{b} \rangle\!\langle \boldsymbol{b} | \Phi_{\mathrm{pas}}(\boldsymbol{V})\) also lie in that subspace. Hence \(\mathcal{M}^{-1}(\boldsymbol{\sigma}) = \mathcal{M}_k^{-1}(\boldsymbol{\sigma})\) and we recover the familiar estimator from 13: \[\label{eq:rdm95shadow95decomp} \begin{align} \widehat{\boldsymbol{D}}_{ij} &\mathrel{\vcenter{:}}= \mathop{\mathrm{tr}}(a_i^\dagger a_j \mathcal{M}^{-1}(\boldsymbol{\sigma}))\\ &= \mathop{\mathrm{tr}}(a_i^\dagger a_j \mathcal{M}_\eta^{-1}(\boldsymbol{\sigma}))\\ &= [\boldsymbol{V}^\dagger ((n+1) \mathop{\mathrm{diag}}(\boldsymbol{b}) - k \mathbb{I}) \boldsymbol{V}]_{ij}. \end{align}\tag{21}\] Now consider arbitrary \(\rho\); each instance of \(\boldsymbol{\sigma}\) still lies in some \(|\boldsymbol{b}|\)-number sector, although \(|\boldsymbol{b}|\) itself is now a random variable. Therefore \(\mathcal{M}^{-1}(\boldsymbol{\sigma}) = \mathcal{M}_{|\boldsymbol{b}|}^{-1}(\boldsymbol{\sigma})\) so replacing \(k \to |\boldsymbol{b}|\) in 21 shows the claim. ◻

Proof. The norm bound \(B\) is the same as in 1 because it still holds that \(\|D\| \leq 1\) and \(\|\widehat{\boldsymbol{D}}\| = \|E(\boldsymbol{b})\| = \max\{n+1-|\boldsymbol{b}|, |\boldsymbol{b}|\} \leq n\).

For the variance, recall 4 except replacing \(\eta\) with \(|\boldsymbol{b}|\): \[E(\boldsymbol{b})^2 = (n + 1 - 2|\boldsymbol{b}|) E(\boldsymbol{b}) + |\boldsymbol{b}|(n+1-|\boldsymbol{b}|) \mathbb{I}.\] Then since \(\widehat{\boldsymbol{D}}^2 = \boldsymbol{V}^\dagger E(\boldsymbol{b})^2 \boldsymbol{V}\), we have \[\label{eq:sigma95bound95perturb} \begin{align} \sigma^2 &= \mathopen{}\left\| \sum_{\ell=1}^N \mathop{\mathrm{\mathbb{E}}}[\boldsymbol{X}_\ell^2] \right\|\mathclose{} = \frac{1}{N} \mathopen{}\left\| {\mathop{\mathrm{\mathbb{E}}}[\widehat{\boldsymbol{D}}^2] - D^2} \right\|\mathclose{}\\ &= \frac{1}{N} \mathopen{}\left\| {\mathop{\mathrm{\mathbb{E}}}[(n+1-2|\boldsymbol{b}|) \widehat{\boldsymbol{D}} + |\boldsymbol{b}|(n+1-|\boldsymbol{b}|) \mathbb{I}] - D^2} \right\|\mathclose{}\\ &= \frac{1}{N} \mathopen{}\left\| {(n+1) D- 2 \mathop{\mathrm{\mathbb{E}}}[|\boldsymbol{b}| \widehat{\boldsymbol{D}}] + \mathopen{}\left( (n+1) \mathop{\mathrm{\mathbb{E}}}|\boldsymbol{b}| -\mathop{\mathrm{\mathbb{E}}}|\boldsymbol{b}|^2 \right)\mathclose{} \mathbb{I}- D^2} \right\|\mathclose{}. \end{align}\tag{22}\] In contrast to the number-conserving case, \(|\boldsymbol{b}|\) is now a random variable, so we need to control its first two moments. These are simply the moments of the particle number operator \(\mathsf{Num}\mathrel{\vcenter{:}}= \sum_{j=1}^n a_j^\dagger a_j\) with respect to the fixed state \(| \psi \rangle\), as the random measurement bases \(\boldsymbol{V}\) commute with \(\mathsf{Num}\). We defer this calculation to 7 below, where we show that \[\begin{align} \langle \psi | \mathsf{Num} | \psi \rangle &\leq 1 + \|\beta\|_F^2, \tag{23}\\ \langle \psi | \mathsf{Num}^2 | \psi \rangle &\leq 2 + 7\|\beta\|_F^2 + \|\beta\|_F^4, \tag{24} \end{align}\] where \(\beta \in \mathbb{C}^{n \times n}\) is given by the representation of \(Z\) described in 33.

To bound the norm of \(\beta\) in terms of \(\delta_{{\mathrm{act}}} \geq \|Z - \mathbb{I}\|\), observe that \[\|\beta\| = \mathopen{}\left\| \begin{pmatrix} 0 & 0\\ \beta & 0 \end{pmatrix} \right\|\mathclose{} = \mathopen{}\left\| \begin{pmatrix} 0 & 0\\ 0 & \mathbb{I} \end{pmatrix} \begin{pmatrix} \alpha - \mathbb{I}& \beta^*\\ \beta & \alpha^* - \mathbb{I} \end{pmatrix} \begin{pmatrix} \mathbb{I}& 0\\ 0 & 0 \end{pmatrix} \right\|\mathclose{} \leq \| \Omega Z \Omega^\dagger - \mathbb{I}\|.\] Hence the Frobenius norm obeys \(\|\beta\|_F^2 \leq n \delta_{{\mathrm{act}}}^2\) and so \[\begin{align} \mathop{\mathrm{\mathbb{E}}}|\boldsymbol{b}| &= \langle \psi | \mathsf{Num} | \psi \rangle \leq 1 + n\delta_{{\mathrm{act}}}^2,\\ \mathop{\mathrm{\mathbb{E}}}|\boldsymbol{b}|^2 &= \langle \psi | \mathsf{Num}^2 | \psi \rangle \leq 2 + 7n\delta_{{\mathrm{act}}}^2 + n^2\delta_{{\mathrm{act}}}^4. \end{align}\] A crude triangle inequality suffices to get an estimate of 22 : \[N\sigma^2 \leq (n+2) + 2 \|{\mathop{\mathrm{\mathbb{E}}}[|\boldsymbol{b}| \widehat{\boldsymbol{D}}]}\| + (n+1) \mathop{\mathrm{\mathbb{E}}}|\boldsymbol{b}| + \mathop{\mathrm{\mathbb{E}}}|\boldsymbol{b}|^2.\] To handle the cross term we apply Cauchy–Schwarz for expectations: \[\begin{align} \|{\mathop{\mathrm{\mathbb{E}}}[|\boldsymbol{b}| \widehat{\boldsymbol{D}}]}\| &= \max_{v \in \mathbb{S}^{n-1}} \mathopen{}\left| {\mathop{\mathrm{\mathbb{E}}}[|\boldsymbol{b}| \langle v | \widehat{\boldsymbol{D}} | v \rangle]} \right|\mathclose{}\\ &\leq \max_{v \in \mathbb{S}^{n-1}} \mathopen{}\left( \sqrt{\mathop{\mathrm{\mathbb{E}}}|\boldsymbol{b}|^2} \sqrt{\mathop{\mathrm{\mathbb{E}}}[\langle v | \widehat{\boldsymbol{D}} | v \rangle^2]} \right)\mathclose{}\\ &\leq \sqrt{\mathop{\mathrm{\mathbb{E}}}|\boldsymbol{b}|^2} \sqrt{\mathop{\mathrm{\mathbb{E}}}\|\widehat{\boldsymbol{D}}\|^2}\\ &\leq n \sqrt{\mathop{\mathrm{\mathbb{E}}}|\boldsymbol{b}|^2}. \end{align}\] The advertised bound for \(\sigma^2\) follows from setting \(\delta_{{\mathrm{act}}} = \frac{c}{\sqrt{n}}\) and some elementary algebra. ◻

Proof. Follows from applying 34 to 8. ◻

Proof. First, we establish ?? : \[\begin{align} \langle \phi | b_j^\dagger b_k | \phi \rangle &= \sum_{p,q=1}^n \langle \phi | \mathopen{}\left( \alpha_{jp}^* a_p^\dagger + \beta_{jp} a_p \right)\mathclose{} \mathopen{}\left( \alpha_{kq} a_q + \beta_{kq}^* a_q^\dagger \right)\mathclose{} | \phi \rangle\\ &= \sum_{p,q=1}^n \langle \phi | \mathopen{}\left( \alpha_{jp}^* \alpha_{kq} a_p^\dagger a_q + \beta_{jp} \beta_{kq}^* a_p a_q^\dagger \right)\mathclose{} | \phi \rangle\\ &= \sum_{p,q=1}^n \mathopen{}\left( \alpha_{jp}^* \alpha_{kq} \langle \phi | a_p^\dagger a_q | \phi \rangle + \beta_{jp} \beta_{kq}^* (\delta_{pq} - \langle \phi | a_q^\dagger a_p | \phi \rangle) \right)\mathclose{}\\ &= [ \alpha^* C \alpha^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \beta^\dagger ]_{jk}, \end{align}\] where we used the fact that \(\langle \phi | a_p a_q | \phi \rangle = \langle \phi | a_p^\dagger a_q^\dagger | \phi \rangle = 0\) on the second line and the CAR \(a_p a_q^\dagger + a_q^\dagger a_p = \delta_{pq}\) on the third. Then the derivation for ?? is analogous, merely swapping the roles of \(\alpha\) and \(\beta\) acting on the right. Finally, ?? is a consequence of the identity \(\langle \phi | b_j b_k^\dagger | \phi \rangle = \delta_{jk} - \langle \phi | b_k^\dagger b_j | \phi \rangle\) (the \(b_j\)’s obey the same CAR because anticommutators are unitarily invariant), and ?? follows from \(\langle \phi | b_j b_k | \phi \rangle = \langle \phi | b_k^\dagger b_j^\dagger | \phi \rangle^*\) (also recall that \(C\) is Hermitian). ◻

Proof. We begin with the first moment. Define \(b_j\) as in 36. Using ?? , we get \[\begin{align} \langle \psi | \mathsf{Num} | \psi \rangle &= \sum_{j=1}^n \langle \psi | a_j^\dagger a_j | \psi \rangle\\ &= \sum_{j=1}^n \langle \phi | b_j^\dagger b_j | \phi \rangle\\ &= \mathop{\mathrm{tr}}(\alpha^* C \alpha^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \beta^\dagger). \end{align}\] We simplify this via the following facts. First, unitarity of \(\Omega Z \Omega^\dagger\) implies that \(\alpha^\mathrm{T}\alpha^* + \beta^\mathrm{T}\beta^* = \mathbb{I}\). Second, \(\mathop{\mathrm{tr}}(C) = 1\) because \(| \phi \rangle\) is a single-particle state. Finally, because both \(\beta^\mathrm{T}\beta^*\) and \(C\) are PSD, \(\mathop{\mathrm{tr}}(\beta^\mathrm{T}\beta^* C) \geq 0\). Altogether, these imply that \[\begin{align} \langle \psi | \mathsf{Num} | \psi \rangle &= \mathop{\mathrm{tr}}((\alpha^\mathrm{T}\alpha^* + \beta^\mathrm{T}\beta^*) C) - 2 \mathop{\mathrm{tr}}(\beta^\mathrm{T}\beta^* C) + \|\beta\|_F^2\\ &\leq 1 + \|\beta\|_F^2, \end{align}\] which is 25 .

Next we consider the second moment. For each term in \(\mathsf{Num}^2 = \sum_{j,k=1}^n a_j^\dagger a_j a_k^\dagger a_k\), we apply Wick’s theorem (35). Recall that the Pfaffian of a \(4 \times 4\) matrix \(S\) is \(\mathop{\mathrm{pf}}(S) = S_{12} S_{34} - S_{13} S_{24} + S_{14} S_{23}\). Hence with ?? , \[\begin{align} \langle \psi | a_j^\dagger a_j a_k^\dagger a_k | \psi \rangle &= \langle \phi | b_j^\dagger b_j b_k^\dagger b_k | \phi \rangle\\ &= \langle \phi | b_j^\dagger b_j | \phi \rangle \langle \phi | b_k^\dagger b_k | \phi \rangle - \langle \phi | b_j^\dagger b_k^\dagger | \phi \rangle \langle \phi | b_j b_k | \phi \rangle + \langle \phi | b_j^\dagger b_k | \phi \rangle \langle \phi | b_j b_k^\dagger | \phi \rangle\\ &= \langle \phi | b_j^\dagger b_j | \phi \rangle \langle \phi | b_k^\dagger b_k | \phi \rangle\\ &\hphantom{=~} - [ \alpha^* C \beta^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \alpha^\dagger ]_{jk} [\beta^* C \alpha^\mathrm{T}+ \alpha (\mathbb{I}- C^\mathrm{T}) \beta^\dagger]_{jk}\\ &\hphantom{=~} + [ \alpha^* C \alpha^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \beta^\dagger ]_{jk} [ \mathbb{I}- \alpha C^\mathrm{T}\alpha^\dagger - \beta^* (\mathbb{I}- C) \beta^\mathrm{T}]_{jk}. \end{align}\] Now sum each term from Wick’s theorem over \(j\) and \(k\). The first simply yields \(\langle \psi | \mathsf{Num} | \psi \rangle^2\). The second is \[\begin{align} &\hphantom{=~} \sum_{j,k=1}^n [ \alpha^* C \beta^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \alpha^\dagger ]_{jk} [\beta^* C \alpha^\mathrm{T}+ \alpha (\mathbb{I}- C^\mathrm{T}) \beta^\dagger]_{jk}\\ &= \sum_{j,k=1}^n [ \alpha^* C \beta^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \alpha^\dagger ]_{jk} [\beta C^* \alpha^\dagger + \alpha^* (\mathbb{I}- C^\dagger) \beta^\mathrm{T}]_{jk}^*\\ &= \sum_{j,k=1}^n [ \alpha^* C \beta^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \alpha^\dagger ]_{jk} [ \beta C^\mathrm{T}\alpha^\dagger + \alpha^* \beta^\mathrm{T}- \alpha^* C \beta^\mathrm{T}]_{jk}^*\\ &= -\sum_{j,k=1}^n [ \alpha^* C \beta^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \alpha^\dagger ]_{jk} [ \alpha^* C \beta^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \alpha^\dagger - \alpha^* \beta^\mathrm{T}- \beta \alpha^\dagger ]_{jk}^*\\ &= -\|\alpha^* C \beta^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \alpha^\dagger\|_F^2, \end{align}\] where on the last line we used that \(\alpha^* \beta^\mathrm{T}+ \beta \alpha^\dagger = 0\) due to unitarity of \(\Omega Z \Omega^\dagger\). The third term from Wick’s theorem is \[\begin{align} &\hphantom{=~} \sum_{j,k=1}^n [ \alpha^* C \alpha^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \beta^\dagger ]_{jk} [ \mathbb{I}- \alpha C^\mathrm{T}\alpha^\dagger - \beta^* (\mathbb{I}- C) \beta^\mathrm{T}]_{jk}\\ &= \sum_{j,k=1}^n [ \alpha^* C \alpha^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \beta^\dagger ]_{jk} [ \mathbb{I}- \alpha^* C^\dagger \alpha^\mathrm{T}- \beta (\mathbb{I}- C^*) \beta^\dagger ]_{jk}^*\\ &= \mathop{\mathrm{tr}}(\alpha^* C \alpha^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \beta^\dagger) - \|\alpha^* C \alpha^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \beta^\dagger\|_F^2\\ &= \langle \psi | \mathsf{Num} | \psi \rangle - \|\alpha^* C \alpha^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \beta^\dagger\|_F^2. \end{align}\] We collect these results to obtain \[\begin{align} \langle \psi | \mathsf{Num}^2 | \psi \rangle &= \langle \psi | \mathsf{Num} | \psi \rangle^2 + \|\alpha^* C \beta^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \alpha^\dagger\|_F^2 + \langle \psi | \mathsf{Num} | \psi \rangle - \|\alpha^* C \alpha^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \beta^\dagger\|_F^2 \notag\\ &\leq (1 + \|\beta\|_F^2)^2 + (1 + \|\beta\|_F^2) + \|\alpha^* C \beta^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \alpha^\dagger\|_F^2.\label{eq:num295intermediate95bound} \end{align}\tag{27}\] It remains to bound the last term. Note that the Frobenius norm has the property \(\|AB\|_F \leq \|A\|_F \|B\|\), which is stronger than mere submultiplicativity. Use this property to get \[\label{eq:fro95norm95aCb} \begin{align} \|\alpha^* C \beta^\mathrm{T}+ \beta (\mathbb{I}- C^\mathrm{T}) \alpha^\dagger\|_F &\leq \|\alpha^* C \beta^\mathrm{T}\|_F + \|\beta (\mathbb{I}- C^\mathrm{T}) \alpha^\dagger\|_F\\ &\leq \|\alpha\| \|C\| \|\beta\|_F + \|\beta\|_F \|\mathbb{I}- C\| \|\alpha\| \\ &\leq 2 \|\beta\|_F, \end{align}\tag{28}\] where we used the additional fact that \(\|\alpha\| \leq 1\) since \(\alpha\) is a block within a unitary matrix. Plug this bound into 27 to conclude the result of 25 . ◻

Proof. Begin by observing that the unitary we apply is \[\Phi(\widehat{Q}_{{\mathrm{act}}}^\mathrm{T}) \Phi(Q) \Phi_{\mathrm{pas}}(\widehat{U}^\dagger) = \Phi(Z) \Phi_{\mathrm{pas}}(U) \Phi_{\mathrm{pas}}(\widehat{U}^\dagger) = \Phi(Z) \Phi_{\mathrm{pas}}(W) \Phi_{\mathrm{pas}}(e^{i\theta} \mathbb{I}),\] where \(\theta \in [-\pi, \pi)\) and \(W \in \mathrm{U}(n)\) are as in 9 . The analysis is therefore identical to 26, up to conjugating \(X\) and \(Y\) by \(\Phi(Z)\). Let us re-define \(| \widetilde{\Psi}(\theta) \rangle \mathrel{\vcenter{:}}= \Phi(Z) \Phi_{\mathrm{pas}}(W) | \Psi(\theta) \rangle\) and write \([\alpha W]_{11} = re^{i\xi}\). We will show shortly that \[\label{eq:perturbed95cos} \begin{align} \langle \widetilde{\Psi}(\theta) | X | \widetilde{\Psi}(\theta) \rangle &= r \cos(\theta + \xi),\\ \langle \widetilde{\Psi}(\theta) | Y | \widetilde{\Psi}(\theta) \rangle &= r \sin(\theta + \xi). \end{align}\tag{29}\] The proof then follows 26 2 except that now \(|r - 1| \leq \|\alpha W - \mathbb{I}\| \leq \|W - \mathbb{I}\| + \|\alpha - \mathbb{I}\| \leq 2 \varepsilon\).

We conclude with the derivation of 29 . Embed \(Z\) into \(\mathrm{SO}(2n+2)\) such that it acts on the ancilla mode trivially. In the \((\alpha, \beta)\) representation, the Bogoliubov transformation is \[\Phi(Z)^\dagger a_1^\dagger a_{n+1}^\dagger \Phi(Z) = \sum_{m=1}^n (\alpha_{1m}^* a_m^\dagger + \beta_{1m} a_m) a_{n+1}^\dagger.\] The expectation of the \(a_m^\dagger a_{n+1}^\dagger\) terms was computed in 24: \[\langle \Psi(\theta) | \Phi_{\mathrm{pas}}(W)^\dagger a_m^\dagger a_{n+1}^\dagger \Phi_{\mathrm{pas}}(W) | \Psi(\theta) \rangle = \frac{1}{2} e^{-i\theta} W_{m1}^*.\] The \(a_m a_{n+1}^\dagger\) term follows an analogous calculation: \[\begin{align} \Phi_{\mathrm{pas}}(W)^\dagger a_m a_{n+1}^\dagger \Phi_{\mathrm{pas}}(W) &= \sum_{j,k=1}^{n+1} [W \oplus 1]_{mj} [W \oplus 1]_{n+1,k}^* a_j a_k^\dagger \end{align}\] and \[\begin{align} \langle \Psi(\theta) | a_j a_k^\dagger | \Psi(\theta) \rangle &= \frac{1}{2} \mathopen{}\left( \langle 0^{n+1} | a_j a_k^\dagger | 0^{n+1} \rangle + e^{i\theta} \langle 0^{n+1} | a_j a_k^\dagger | 1_1 1_{n+1} \rangle \right.\mathclose{}\\ &\hphantom{=~} \mathopen{}\left. + \, e^{-i\theta} \langle 1_1 1_{n+1} | a_j a_k^\dagger | 0^{n+1} \rangle + \langle 1_1 1_{n+1} | a_j a_k^\dagger | 1_1 1_{n+1} \rangle \right)\mathclose{}\\ &= \frac{1}{2} (\delta_{jk} - \langle 1_1 1_{n+1} | a_k^\dagger a_j | 1_1 1_{n+1} \rangle) = B_{jk}, \end{align}\] where we define the matrix \(B \mathrel{\vcenter{:}}= \frac{1}{2} \mathop{\mathrm{diag}}(0, 1, \ldots, 1, 0) \in \mathbb{R}^{(n+1) \times (n+1)}\). Hence \[\begin{align} \langle \Psi(\theta) | \Phi_{\mathrm{pas}}(W)^\dagger a_m a_{n+1}^\dagger \Phi_{\mathrm{pas}}(W) | \Psi(\theta) \rangle = [(W \oplus 1) B (W^\dagger \oplus 1)]_{m,n+1} = 0 \end{align}\] since the last column of \((W \oplus 1) B (W^\dagger \oplus 1)\) is \(0\). Altogether, we get \[\langle \Psi(\theta) | \Phi_{\mathrm{pas}}(W)^\dagger \Phi(Z)^\dagger a_1^\dagger a_{n+1}^\dagger \Phi(Z) \Phi_{\mathrm{pas}}(W) | \Psi(\theta) \rangle = \frac{1}{2} e^{-i\theta} [\alpha W]_{11}^*. \qedhere\] ◻

Proof. The algorithm begins by running \(\mathsf{GaussianTomo}\) on \(N_{{\mathrm{act}}}\) copies of the state \(\Phi(Q) | 0^{n} \rangle\). By 3, if we take \(N_{{\mathrm{act}}} = \mathopen{}\left\lceil \frac{32 n^2 \log(4n/\eta_{{\mathrm{act}}})}{\delta_{{\mathrm{act}}}^2} \right\rceil\mathclose{}\) then there exists a symplectic \(\boldsymbol{Q}_{{\mathrm{pas}}} \in \mathrm{O}(n) \cap \mathrm{Sp}(n, \mathbb{R})\) such that \(\|\widehat{\boldsymbol{Q}}_{{\mathrm{act}}} \boldsymbol{Q}_{{\mathrm{pas}}} - Q\| \leq \delta_{{\mathrm{act}}}\), except with probability \(\eta_{{\mathrm{act}}}\).

Define \(\boldsymbol{Z} \mathrel{\vcenter{:}}= \widehat{\boldsymbol{Q}}_{{\mathrm{act}}}^\mathrm{T}Q \boldsymbol{Q}_{{\mathrm{pas}}}^\mathrm{T}\) and let \(\boldsymbol{U}\) be the \(\mathrm{U}(n)\)-representation of \(\boldsymbol{Q}_{{\mathrm{pas}}}\). Condition on the success of the previous step. The circuit \(\boldsymbol{\mathcal{C}} \mathrel{\vcenter{:}}= \Phi(\widehat{\boldsymbol{Q}}_{{\mathrm{act}}}^\mathrm{T}) \Phi(Q)\) is equivalent to \(\Phi(\boldsymbol{Z}) \Phi_{\mathrm{pas}}(\boldsymbol{U})\), which we input into \(\mathsf{PassiveTomo}\). By convention, \(N_{{\mathrm{pas}}}\) is the number of copies per state of the form \(| \psi_j \rangle \mathrel{\vcenter{:}}= \boldsymbol{\mathcal{C}} | 1_j \rangle\) and \(| \widetilde{\psi}_j \rangle \mathrel{\vcenter{:}}= \boldsymbol{\mathcal{C}} \Phi_{\mathrm{pas}}(F^\dagger) | 1_j \rangle\) prepared by \(\mathsf{PassiveTomo}\), for a total of \(2n N_{{\mathrm{pas}}}\) queries to \(\Phi(Q)\). The error analysis for this subroutine is as follows. Set \(\delta_{{\mathrm{act}}} = \frac{c}{\sqrt{n}}\) for some small \(c < 1\) to be determined later. Then we can use 5, which implies that \(N_{{\mathrm{pas}}} = \mathopen{}\left\lceil \frac{48n \log(2n/\eta_{{\mathrm{pas}}})}{\delta_{{\mathrm{pas}}}^2} \right\rceil\mathclose{}\) copies suffices to learn an RDM to error \(\delta_{{\mathrm{pas}}}\) in operator norm, except with probability \(\eta_{{\mathrm{pas}}}\). We set \(\delta_{{\mathrm{pas}}} = \frac{c}{\sqrt{n}}\) as well, allowing us to use 31 to find a unitary \(\widehat{\boldsymbol{U}}\) such that \(\mathop{\mathrm{\mathsf{dist}_{ph}}}(\widehat{\boldsymbol{U}}, \boldsymbol{U}) \leq 200\sqrt{2}c\) except with probability \(2n\eta_{{\mathrm{pas}}}\).

The final piece of \(\mathsf{PassiveTomo}\) is the \(\mathrm{U}(1)\) phase estimation. This returns a phase \(\widehat{\boldsymbol{\theta}} \in (-\pi, \pi]\) such that, by 37, the unitary \(\boldsymbol{U}^\sharp \mathrel{\vcenter{:}}= e^{i\widehat{\boldsymbol{\theta}}} \widehat{\boldsymbol{U}}\) obeys \(\|\boldsymbol{U}^\sharp - \boldsymbol{U}\| \leq 1800\sqrt{2}c\). Conditioned on all prior steps succeeding, this holds with probability \(1 - 2\eta_{{\mathrm{ph}}}\) if we make \(2N_{{\mathrm{ph}}} = 2\mathopen{}\left\lceil \frac{(6 + 4\sqrt{2}) \log(2/\eta_{{\mathrm{ph}}})}{80000 c^2} \right\rceil\mathclose{}\) queries to \(\Phi(Q)\) (26). The unconditional success probability is therefore at least \(1 - \eta_{\mathrm{act}}- 2n\eta_{\mathrm{pas}}- 2\eta_{\mathrm{ph}}\) by a union bound, to achieve an error of \[\|\widehat{\boldsymbol{Q}}_{\mathrm{act}}\widehat{\boldsymbol{Q}}_{\mathrm{pas}}- Q\| \leq \frac{c}{\sqrt{n}} + 1800\sqrt{2}c.\] Choosing \(c = \frac{\varepsilon}{2600}\) is more than enough to bound this by \(\varepsilon\), and choosing \(\eta_{\mathrm{act}}= \frac{\delta}{3}\), \(\eta_{\mathrm{pas}}= \frac{\delta}{6n}\), and \(\eta_{\mathrm{ph}}= \frac{\delta}{6}\) bounds the failure probability by \(\delta\). The resulting query complexity is \[\begin{align} N_{\mathrm{act}}+ 2n N_{\mathrm{pas}}+ 2 N_{\mathrm{ph}}&= \mathopen{}\left\lceil \frac{32 n^2 \log(4n/\eta_{{\mathrm{act}}})}{\delta_{{\mathrm{act}}}^2} \right\rceil \mathclose{}+ 2n \mathopen{}\left\lceil \frac{48n \log(2n/\eta_{{\mathrm{pas}}})}{\delta_{{\mathrm{pas}}}^2} \right\rceil \mathclose{}+ 2 \mathopen{}\left\lceil \frac{(6 + 4\sqrt{2}) \log(2/\eta_{{\mathrm{ph}}})}{80000 c^2} \right\rceil \mathclose{}\notag\\ &\leq \mathopen{}\left\lceil \frac{C_1 n^3 \log(K n/\delta)}{\varepsilon^2} \right\rceil \mathclose{}+ 2n \mathopen{}\left\lceil \frac{C_2 n^2 \log(K n^2/\delta)}{\varepsilon^2} \right\rceil \mathclose{}+ 2 \mathopen{}\left\lceil \frac{C_3 \log(K/\delta)}{\varepsilon^2} \right\rceil\mathclose{}, \label{eq:base95queries} \end{align}\tag{30}\] where one can take \(C_1 = 2.2 \times 10^8\), \(C_2 = 3.3 \times 10^8\), \(C_3 = 1000\), and \(K = 12\). ◻

Proof. As 7 indicates, \(\mathsf{ActiveTomo}(\Phi(Q), \frac{1}{10}, \delta)\) makes \(\mathcal{O}(n^3 \log(n/\delta))\) queries. Hence by 21 the bootstrapped process with error \(\varepsilon/n\) makes a total of \(\mathcal{O}(n^4 \log(n/\delta) / \varepsilon)\) queries. This error in the operator norm of the \(\mathrm{O}(2n)\)-representation is chosen such that the diamond distance error is at most \(\varepsilon\), per 10. See 28 for the gate complexity argument (note that both passive and active FLOs use \(\mathcal{O}(n^2)\) gates).

For the classical cost, we have that the circuit for \(\Phi(\boldsymbol{V}_t^\dagger)\) can be determined in \(\mathcal{O}(n^3)\) time (this only needs to be calculated once per iteration). For each iteration \(t\),

1. \(\mathsf{GaussianTomo}\) and computing the normal form costs \(\mathcal{O}(N_{\mathrm{act}}n^\omega + n^3) = \mathcal{O}(n^{\omega+3} \log(n/\delta_t))\) operations (42);

2. Determining the circuit for \(\Phi(\widehat{\boldsymbol{Q}}_{{\mathrm{act}}})\) costs \(\mathcal{O}(n^3)\) operations;

3. \(\mathsf{PassiveTomo}\) costs \(\mathcal{O}(n^4 \log(n/\delta_t))\) operations (from 28, adjusted to account for the fact that the perturbed RDM estimates are rank-\(\mathcal{O}(n)\) rather than rank-\(1\));

4. Forming \(\widehat{\boldsymbol{Q}}_{{\mathrm{act}}} \widehat{\boldsymbol{Q}}_{{\mathrm{pas}}}\) costs \(\mathcal{O}(n^\omega)\) operations.

The \(\mathsf{GaussianTomo}\) step asymptotically dominates, so the total time complexity is \[\mathcal{O}\mathopen{}\left( \sum_{t=0}^T n^{\omega+3} \log(n/\delta_t) \right)\mathclose{} = \mathcal{O}\mathopen{}\left( n^{\omega+3} (\log(n/\delta) T + T^2) \right)\mathclose{} = \mathcal{O}\mathopen{}\left( n^{\omega+3} \log^2(n/{\min\{\varepsilon, \delta\}}) \right)\mathclose{}. \qedhere\] ◻

Proof. The formula can be derived from any of the aforementioned papers on matchgate shadows; we will follow [38] due to its relatively compact presentation. There they show that, for any quadratic Majorana observable \(-i \gamma_j \gamma_k\), with \(j \neq k\), \[\mathop{\mathrm{tr}}(-i \gamma_j \gamma_k \mathcal{M}^{-1}(\Phi(\boldsymbol{R})^\dagger | \boldsymbol{b} \rangle\!\langle \boldsymbol{b} | \Phi(\boldsymbol{R}))) = (2n-1) \mathop{\mathrm{pf}}((\boldsymbol{R}^\mathrm{T}J(\boldsymbol{b}) \boldsymbol{R})[(j, k)]),\] where \(A[(j, k)]\) denotes the \((j, k)\)-principal submatrix of a matrix \(A\). For a \(2 \times 2\) skew-symmetric matrix, the Pfaffian is simply the upper-right entry, so the right-hand side is \((2n-1) [\boldsymbol{R}^\mathrm{T}J(\boldsymbol{b}) \boldsymbol{R}]_{jk} = \widehat{\boldsymbol{\Gamma}}_{jk}\). The left-hand side is \(\mathop{\mathrm{tr}}(-i \gamma_j \gamma_k \rho) = \Gamma_{jk}\) in expectation, by construction of the classical shadows. ◻

Proof. Since \(\boldsymbol{R}\) and \(J(\boldsymbol{b})\) are both orthogonal, so too is \(\boldsymbol{R}^\mathrm{T}J(\boldsymbol{b}) \boldsymbol{R}\). Recall also that any covariance matrix obeys \(\|\Gamma\| \leq 1\). Thus for any \(\ell \in [N]\), \[\|\boldsymbol{X}_\ell\| \leq \frac{\|\widehat{\boldsymbol{\Gamma}}_\ell\| + \|\Gamma\|}{N} \leq \frac{(2n-1) + 1}{N} = \frac{2n}{N}.\] For the variance, observe that \(J(\boldsymbol{b})^2 = -\mathbb{I}\), so \[\sigma^2 = \mathopen{}\left\| \frac{1}{N^2} \sum_{\ell=1}^N (\mathop{\mathrm{\mathbb{E}}}[\widehat{\boldsymbol{\Gamma}}_\ell^2] - \Gamma^2) \right\|\mathclose{} = \frac{1}{N} \mathopen{}\left\| -(2n-1)^2 \mathbb{I}- \Gamma^2 \right\|\mathclose{} \leq \frac{(2n-1)^2 + 1}{N}.\] ◻

Proof. Let \(\widehat{\boldsymbol{\Gamma}}_1, \ldots, \widehat{\boldsymbol{\Gamma}}_N\) as in 8 and set \(\overline{\boldsymbol{\Gamma}} \mathrel{\vcenter{:}}= \frac{1}{N} \sum_{\ell=1}^N \widehat{\boldsymbol{\Gamma}}_\ell\). We recall the matrix Bernstein inequality (8), where note that the matrices have linear dimension \(2n\) and the \(B\) and \(\sigma^2\) parameters are given by 8: \[\Pr\mathopen{}\left( \|\overline{\boldsymbol{\Gamma}} - \Gamma\| \geq \varepsilon \right)\mathclose{} \leq 4n \exp\mathopen{}\left( \frac{-N \varepsilon^2 / 2}{(2n-1)^2 + 1 + 2n\varepsilon/3} \right)\mathclose{}.\] For \(n \geq 1\) and \(\varepsilon < 1\), taking \(N\) as in 32 suffices to bound this probability by \(\delta\). ◻

Proof. First we check the runtime of 8. Any FLO requires \(\mathcal{O}(n^2)\) elementary gates, and random instances can be constructed in \(\mathcal{O}(n^2)\) time [61]. Each repetition is dominated by the matrix multiplication of \(\boldsymbol{R}^\mathrm{T}J(\boldsymbol{b}) \boldsymbol{R}\). The algorithm concludes by computing the normal form of \(\overline{\boldsymbol{\Gamma}}\), which takes \(\mathcal{O}(n^3)\) time (5).

For the copy complexity, let \(\Gamma\), \(\overline{\boldsymbol{\Gamma}}\), and \(\boldsymbol{\Gamma}^\star\) be the covariance matrix of \(| \psi \rangle\), the unrounded estimate, and the rounded estimate corresponding to Gaussian \(| \widehat{\boldsymbol{\psi}} \rangle\), respectively. By 6 15, we have \[\mathop{\mathrm{\mathsf{dist}_{tr}}}(| \widehat{\boldsymbol{\psi}} \rangle, | \psi \rangle) \leq \frac{1}{4} \|\boldsymbol{\Gamma}^\star - \Gamma\|_F \leq \frac{3}{4} \|\overline{\boldsymbol{\Gamma}} - \Gamma\|_F \leq \frac{3\sqrt{2n}}{4} \|\overline{\boldsymbol{\Gamma}} - \Gamma\|.\] Per 41, if we take \(N = \mathopen{}\left\lceil \frac{9n^3 \log(4n/\delta)}{\varepsilon^2} \right\rceil\mathclose{}\) then ?? holds. ◻

Proof. Write the vacuum state as \[| 0^{2n} \rangle\!\langle 0^{2n} | = \prod_{j=1}^n \mathopen{}\left( \frac{\mathbb{I}- i \gamma_j \gamma_{j + n}}{2} \right)\mathclose{} \mathopen{}\left( \frac{\mathbb{I}- i \gamma_{j+2n} \gamma_{j + 3n}}{2} \right)\mathclose{}.\] Graphically, this is corresponds to a perfect matching on \([4n]\), as is 33 . Let \(\pi \in \mathcal{S}_{4n}\) be the permutation that maps between the two perfect matchings via \(j+n \leftrightarrow j+2n\). This consists of \(n\) swaps, so if \(P_\pi \in \mathrm{O}(4n)\) is the matrix representation of \(\pi\), then the corresponding unitary is \(\Phi(P_\pi)\) with \(\det(P_\pi) = (-1)^n\). To ensure even parity for all \(n\), we can additionally flip every other matching by another introducing another permutation \(\sigma \in \mathcal{S}_{4n}\) that performs \(j \leftrightarrow j+2n\) if and only if \(j\) is odd. This is equivalent to the staggered sign appearing in 33 since Majoranas anticommute. Overall, this implies that \(\Phi(P_\sigma P_\pi) | 0^{2n} \rangle = | \mathrm{fEPR} \rangle\) where \(\det(P_\sigma P_\pi) = (-1)^{2n} = 1\). ◻

Proof. Let \(1 \leq j < k \leq 2n\). Use the fact that \(\widetilde{Q}\) acts trivially on the auxiliary modes and that Majorana monomials are trace-orthogonal to get \[\Gamma_{j,k+2n} = -i \langle \mathrm{fEPR}(Q) | \gamma_j \gamma_{k+2n} | \mathrm{fEPR}(Q) \rangle = \frac{(-1)^{k}}{2^{2n}} \mathop{\mathrm{tr}}\mathopen{}\left( \gamma_j \Phi(\widetilde{Q}) \gamma_k \Phi(\widetilde{Q})^\dagger \right)\mathclose{} = (-1)^k Q_{jk}.\] Meanwhile the diagonal blocks of \(\Gamma\) vanish because \(| \mathrm{fEPR}(Q) \rangle\) is pure Gaussian, so \(\Gamma\) must be orthogonal. ◻

Proof. Use 41 to get an \((\varepsilon/2)\)-estimate \(\overline{\boldsymbol{\Gamma}}\) of the covariance matrix \(\Gamma\) of \(| \mathrm{fEPR}(Q) \rangle\). This uses \(N = \mathopen{}\left\lceil \frac{128 n^2 \log(8n/\delta)}{\varepsilon^2} \right\rceil\mathclose{}\) copies. By 6, if we extract the top-right block of \(\overline{\boldsymbol{\Gamma}}\), round it to a nearby orthogonal matrix (e.g., by taking the SVD), and right-multiply by \(S\), then the solution satisfies 34 . ◻

Proof. Run \(\mathsf{Bootstrap}(\mathop{\mathrm{\mathcal{A}}}; \frac{\varepsilon}{n}, \delta)\) as in 2, where \(\mathop{\mathrm{\mathcal{A}}}\) is the algorithm described in 9. By 21, the base cost for constant error \(\frac{1}{10}\) and failure probability \(\delta_t = \frac{\delta}{2^{T+1-t}}\) is \(\mathcal{O}(n^2 \log(n/\delta_t))\) queries, so the entire procedure uses \(\mathcal{O}(n^3 \log(n/\delta) / \varepsilon)\) queries. That all operations are parity-conserving Gaussian follows from the fact that the unitary which prepares \(| \mathrm{fEPR} \rangle\) is an \(\mathrm{SO}(4n)\) FLO (43). ◻

Proof (of 22). [28] begin by writing the output of each state tomography as \[| \widehat{\boldsymbol{u}}_j \rangle = e^{i\boldsymbol{\alpha}_j} \sqrt{1 - \boldsymbol{\varepsilon}_j} | u \rangle + \sqrt{\boldsymbol{\varepsilon}_j} | \boldsymbol{w} \rangle\] where \(\boldsymbol{\varepsilon}_j \leq \varepsilon^2\) with probability at least \(1 - \frac{\delta}{2n}\) and \(| \boldsymbol{w} \rangle\) is Haar-random on the subspace orthogonal to \(| u_j \rangle\). By 19, we can guarantee this with \(N = \mathopen{}\left\lceil \frac{384(11n + 5\log(4n/\delta)}{\varepsilon^2} \right\rceil\mathclose{}\) copies of \(\Phi_{\mathrm{pas}}(U) | 1_j \rangle\). They then re-express this as \[\widehat{\boldsymbol{U}} - U \boldsymbol{A} = U \boldsymbol{A} \boldsymbol{\Delta} + \boldsymbol{W} \boldsymbol{E},\] where \(\boldsymbol{W}\) is the random matrix with \(| \boldsymbol{w}_j \rangle\) as its columns, and we have defined \(\boldsymbol{A} \mathrel{\vcenter{:}}= \mathop{\mathrm{diag}}(e^{i\boldsymbol{\alpha}_1}, \ldots, e^{i\boldsymbol{\alpha}_n})\), \(\boldsymbol{\Delta} \mathrel{\vcenter{:}}= \mathop{\mathrm{diag}}(\sqrt{1 - \boldsymbol{\varepsilon}_1}, \ldots, \sqrt{1 - \boldsymbol{\varepsilon}_n}) - \mathbb{I}\), and \(\boldsymbol{E} \mathrel{\vcenter{:}}= \mathop{\mathrm{diag}}(\sqrt{\boldsymbol{\varepsilon}_1}, \ldots, \sqrt{\boldsymbol{\varepsilon}_n})\). By a union bound, both \(\|\boldsymbol{\Delta}\|\) and \(\|\boldsymbol{E}\|\) are at most \(\varepsilon\) except with probability \(\delta/2\). This implies that \[\label{eq:bound95with95W} \min_{\Theta \in \mathop{\mathrm{diag}}(\mathbb{R}^n)} \|\widehat{\boldsymbol{U}} - U e^{i\Theta}\| \leq (1 + \|\boldsymbol{W}\|) \varepsilon \quad \text{with probability at least } 1 - \frac{\delta}{2}.\tag{35}\] Using techniques from random matrix theory, [28] argue that the norm of \(\boldsymbol{W}\) is bounded by some unspecified constant with high constant probability, say \(\geq 0.98\). To boost this probability also to \(\geq 1 - \frac{\delta}{2}\) they use a standard median-of-means trick, repeating the column tomography process \(\mathcal{O}(\log(1/\delta))\) times.

We provide an alternative proof of the statement here which gets the \(\delta\) failure probability directly. Observe that the columns of \(\boldsymbol{W}\) are:

1. Independent,

2. Subgaussian with Orlicz \(\psi_2\)-norm \(\| | \boldsymbol{w}_j \rangle \|_{\psi_2} \leq \frac{1}{\sqrt{n-1}}\), and

3. Isotropic on average: \(\mathop{\mathrm{\mathbb{E}}}[\boldsymbol{W} \boldsymbol{W}^\dagger] = \mathbb{I}\).

The first point is by construction. The second is because each \(| \boldsymbol{w}_j \rangle\) is uniform on the sphere orthogonal to \(| u_j \rangle\), hence subgaussian; we show at the end how to derive the constant in the \(\psi_2\)-norm. The third follows from the fact that \(\mathop{\mathrm{\mathbb{E}}}| \boldsymbol{w}_j \rangle\!\langle \boldsymbol{w}_j |\) is equal to the normalized projector onto the subspace orthogonal to \(| u_j \rangle\): \[\begin{align} \mathop{\mathrm{\mathbb{E}}}[\boldsymbol{W} \boldsymbol{W}^\dagger] = \mathop{\mathrm{\mathbb{E}}}\sum_{j=1}^n | \boldsymbol{w}_j \rangle\!\langle \boldsymbol{w}_j | = \sum_{j=1}^n \frac{\mathbb{I}- | u_j \rangle\!\langle u_j |}{n - 1} = \mathbb{I}. \end{align}\]

Now we use random matrix theory to bound the norm of \(\boldsymbol{W}\). The columns of \(\boldsymbol{W}\) are not exactly isotropic, but only isotropic on averge, so we need to use a non-isotropic concentration bound appearing in [68]: for every \(t \geq 0\), \[\label{eq:non-isotropic95bound} \Pr\mathopen{}\left( \mathopen{}\left\| \frac{1}{n} \boldsymbol{W} \boldsymbol{W}^\dagger - \frac{1}{n} \mathbb{I}\right\|\mathclose{} \leq \max\{\gamma, \gamma^2\} \right)\mathclose{} \geq 1 - 2\exp\mathopen{}\left( -\frac{c_1 t^2}{K^4} \right)\mathclose{} \quad \text{where } \gamma = K^2 \sqrt{\frac{\log 9}{c_1}} + \frac{t}{\sqrt{n}},\tag{36}\] \(K = \max_{j \in [n]} \| | \boldsymbol{w}_j \rangle \|_{\psi_2}\), and one can choose \(c_1 = \frac{1}{128e^2}\). This implies that with the same probability, \[\|\boldsymbol{W}\| \leq \sqrt{1 + n \max\{\gamma, \gamma^2\}} \leq 1 + \frac{1}{2} n \gamma,\] assuming that \(n\) is sufficiently large enough so that \(\gamma \leq 1\). Then, it suffices to set \(t = K^2 \sqrt{\frac{\log(4/\delta)}{c_1}}\) to get \[\|\boldsymbol{W}\| \leq 1 + \frac{1}{2} n K^2 \mathopen{}\left( \sqrt{\frac{\log 9}{c_1}} + \sqrt{\frac{\log(4/\delta)}{c_1 n}} \right)\mathclose{}\] except with probability \(\delta/2\). Assuming that the final failure probability \(\delta\) (by a union bound with the event in 35 ) is no less than \(e^{-5n}\),¹³ we can conclude that \[\min_{\Theta \in \mathop{\mathrm{diag}}(\mathbb{R}^n)} \|\widehat{\boldsymbol{U}} - U e^{i\Theta}\| \leq 120 \varepsilon \quad \text{with probability at least } 1 - \delta.\] Rescaling \(\varepsilon\) implies that the constant \(C\) in [line:phaseless95N] of 3 is no larger than \(5.6 \times 10^6\).

It remains to establish the constant in [item:psi295norm]. It is a standard fact that Haar-random unit vectors \(| \boldsymbol{w} \rangle\) in \(\mathbb{C}^d\) are subgaussian with \(\|| \boldsymbol{w} \rangle\|_{\psi_2} \lesssim \frac{1}{\sqrt{d}}\) [69]; we make this constant explicit. First, the subgaussian (or Orlicz \(\psi_2\)-)norm of a random scalar variable \(\boldsymbol{X}\) is defined as \[\|\boldsymbol{X}\|_{\psi_2} \mathrel{\vcenter{:}}= \inf\{ b > 0 : \mathop{\mathrm{\mathbb{E}}}e^{\boldsymbol{X}^2 / b^2} \leq 2 \}.\] The generalization to random vectors is then \(\|| \boldsymbol{w} \rangle\|_{\psi_2} \mathrel{\vcenter{:}}= \sup_{| v \rangle} \| \langle v | \boldsymbol{w} \rangle \|_{\psi_2}\). Denoting \(\boldsymbol{X} = |\langle v | \boldsymbol{w} \rangle|\), we expand the exponential: \[\mathop{\mathrm{\mathbb{E}}}e^{\boldsymbol{X}^2 / b^2} = 1 + \sum_{p=1}^\infty \frac{1}{p!} \frac{\mathop{\mathrm{\mathbb{E}}}[\boldsymbol{X}^{2p}]}{b^{2p}}.\] The moments of \(|\langle v | \boldsymbol{w} \rangle|^{2p}\) are well-known; for example, using [70] gets \[\mathop{\mathrm{\mathbb{E}}}[|\langle v | \boldsymbol{w} \rangle|^{2p}] = \frac{1}{\binom{p + d - 1}{p}}.\] Using a computer algebra system, we find that \[\sum_{p=1}^\infty \frac{1}{p! \binom{p + d - 1}{p} b^{2p}} = e^{b^2/2} b^{2(d-1)} ((d-1)! - \Gamma(d, b^{-2})),\] where \(\Gamma(s, x) = \int_x^\infty t^{s-1} e^{-t} \, dt\) is the incomplete Gamma function; all we use is that \(\Gamma(s, x) \geq 0\) for real arguments. Write \(b = \sqrt{\frac{c}{d}}\) for some constant \(c > 0\) to be determined. Using Stirling’s approximation, \[\begin{align} \sum_{p=1}^\infty \frac{1}{p! \binom{p + d - 1}{p} b^{2p}} &\leq \exp\mathopen{}\left( \frac{c}{2d} \right)\mathclose{} \mathopen{}\left( \frac{c}{d} \right)^{d-1}\mathclose{} \sqrt{2\pi(d-1)} \mathopen{}\left( \frac{d-1}{e} \right)^{d-1}\mathclose{} \exp\mathopen{}\left( \frac{1}{12(d-1)} \right)\mathclose{}\\ &\leq \exp\mathopen{}\left( \frac{6c + 1}{12(d-1)} \right)\mathclose{} \mathopen{}\left( \frac{e}{c} \right)^{-(d-1)}\mathclose{} \sqrt{2\pi(d-1)}. \end{align}\] We can choose \(c = 1\) for simplicity; for all \(d \geq 2\) this bound is at most \(\sqrt{2}e^{-41/24} < 0.26\). Hence \(\mathop{\mathrm{\mathbb{E}}}e^{\boldsymbol{X}^2 / b^2} < 1.26 < 2\) for \(b = \sqrt{\frac{1}{d}}\), making this is a valid bound on the subgaussian norm. We apply this to \(| \boldsymbol{w}_j \rangle\) with \(d = n - 1\). ◻

Proof (of 33). Define the following “vector-of-operators” notation: \[\vec{A} \mathrel{\vcenter{:}}= \begin{pmatrix} \vec{A}_1\\ \vec{A}_2 \end{pmatrix} \quad \text{where } \vec{A}_1 \mathrel{\vcenter{:}}= \begin{pmatrix} a_1\\ \vdots\\ a_n \end{pmatrix} \text{ and } \vec{A}_2 \mathrel{\vcenter{:}}= \begin{pmatrix} a_1^\dagger\\ \vdots\\ a_n^\dagger \end{pmatrix},\] \[\vec{\gamma} \mathrel{\vcenter{:}}= \begin{pmatrix} \vec{\gamma}_1\\ \vec{\gamma}_2 \end{pmatrix} \quad \text{where } \vec{\gamma}_1 \mathrel{\vcenter{:}}= \begin{pmatrix} \gamma_1\\ \vdots\\ \gamma_n \end{pmatrix} \text{ and } \vec{\gamma}_2 \mathrel{\vcenter{:}}= \begin{pmatrix} \gamma_{n+1}\\ \vdots\\ \gamma_{2n} \end{pmatrix}.\] These two are related via \[\vec{A} = \frac{1}{\sqrt{2}} \Omega \vec{\gamma} \quad \text{where } \Omega = \frac{1}{\sqrt{2}} \begin{pmatrix} \mathbb{I}& i\mathbb{I}\\ \mathbb{I}& -i\mathbb{I} \end{pmatrix}.\] We also use the notation that operators transform elementwise, e.g., \[\Phi(Z)^\dagger \vec{\gamma} \Phi(Z) = \begin{pmatrix} \Phi(Z)^\dagger \gamma_1 \Phi(Z)\\ \vdots\\ \Phi(Z)^\dagger \gamma_{2n} \Phi(Z) \end{pmatrix} = Z \vec{\gamma}.\] Write \(Z\) in \(n \times n\) blocks: \[Z = \begin{pmatrix} Z_{11} & Z_{12}\\ Z_{21} & Z_{22} \end{pmatrix}.\] Then we straightforwardly compute: \[\begin{align} \Phi(Z)^\dagger \vec{A} \Phi(Z) &= \frac{1}{\sqrt{2}} \Phi(Z)^\dagger (\Omega \vec{\gamma}) \Phi(Z)\\ &= \frac{1}{2} \Phi(Z)^\dagger \begin{pmatrix} \vec{\gamma}_1 + i\vec{\gamma}_2\\ \vec{\gamma}_1 - i\vec{\gamma}_2 \end{pmatrix} \Phi(Z)\\ &= \frac{1}{2} \begin{pmatrix} Z_{11} \vec{\gamma}_1 + Z_{12} \vec{\gamma}_2 + iZ_{21} \vec{\gamma}_1 + iZ_{22} \vec{\gamma}_2\\ Z_{11} \vec{\gamma}_1 + Z_{12} \vec{\gamma}_2 - iZ_{21} \vec{\gamma}_1 - iZ_{22} \vec{\gamma}_2 \end{pmatrix}\\ &= \frac{1}{2} \begin{pmatrix} (Z_{11} + iZ_{21}) (\vec{A}_1 + \vec{A}_2) -i (Z_{12} + iZ_{22})(\vec{A}_1 - \vec{A}_2)\\ (Z_{11} - iZ_{21}) (\vec{A}_1 + \vec{A}_2) -i (Z_{12} - iZ_{22})(\vec{A}_1 - \vec{A}_2) \end{pmatrix}\\ &= \frac{1}{2} \begin{pmatrix} [Z_{11} + Z_{22} - i(Z_{12} - Z_{21})] \vec{A}_1 + [Z_{11} - Z_{22} + i(Z_{12} + Z_{21})] \vec{A}_2\\ [Z_{11} - Z_{22} - i(Z_{12} + Z_{21})] \vec{A}_1 + [Z_{11} + Z_{22} + i(Z_{12} - Z_{21})] \vec{A}_2 \end{pmatrix}\\ &= \begin{pmatrix} \alpha & \beta^*\\ \beta & \alpha^* \end{pmatrix} \vec{A}. \end{align}\] where we have defined \(\alpha \mathrel{\vcenter{:}}= \frac{1}{2} [Z_{11} + Z_{22} - i(Z_{12} - Z_{21})]\) and \(\beta \mathrel{\vcenter{:}}= \frac{1}{2} [Z_{11} - Z_{22} - i(Z_{12} + Z_{21})]\) in the final line. In particular, one checks that \[\Omega Z \Omega^\dagger = \begin{pmatrix} \alpha & \beta^*\\ \beta & \alpha^* \end{pmatrix}. \qedhere\] ◻