Verifying Multipartite Entangled GHZ States via Multiple Quantum Coherences

TL;DR

This paper introduces a scalable, noise-robust method for verifying multipartite entanglement in quantum systems, demonstrated on a 20-qubit superconducting device with high fidelity for an 18-qubit GHZ state.

Contribution

The paper develops a new entanglement metric based on multiple quantum coherences and demonstrates its effectiveness on a large superconducting quantum device.

Findings

Achieved 0.5165 fidelity for an 18-qubit GHZ state

The metric is robust to noise and easy to implement

Verified multipartite entanglement across 18 qubits

Abstract

The ability to generate and verify multipartite entanglement is an important benchmark for near-term quantum devices devices. We develop a scalable entanglement metric based on multiple quantum coherences, and demonstrate experimentally on a 20-qubit superconducting device - the IBM Q System One. We report a state fidelity of 0.5165$\pm$0.0036 for an 18-qubit GHZ state, indicating multipartite entanglement across all 18 qubits. Our entanglement metric is robust to noise and only requires measuring the population in the ground state; it can be readily applied to other quantum devices to verify multipartite entanglement.

Tables

Table 1. Table 1: Fidelity Comparison for small GHZ states

Method	2q GHZ	3q GHZ	4q GHZ	5q GHZ
QST	0.96	0.93	0.87	0.85
MQC	0.98	0.94	0.87	0.86

Table 2. Table 2: Qubit parameters on IBM System One The qubit frequency, T 1 subscript 𝑇 1 T_{1} , T 2 , echo subscript 𝑇 2 echo T_{2,\text{echo}} , readout fidelity are presented.

Qubit	Frequency (GHz)	$T_1$ ($\mu s$)	$T_{2,\text{echo}}$ ($\mu s$)	Readout Fidelity
Q0	4.666	88.1	76.6	98.1
Q1	4.760	69.0	75.7	96.4
Q2	4.609	58.3	65.4	97.2
Q3	5.031	60.9	73.0	79.7
Q4	4.657	69.1	78.1	96.6
Q5	4.752	74.4	71.9	95.9
Q6	4.829	60.2	65.8	98.1
Q7	4.698	80.7	79.5	96.4
Q8	4.893	64.0	75.7	96.5
Q9	4.731	63.3	70.7	93.0
Q10	4.840	59.1	62.9	96.6
Q11	4.755	64.1	56.3	97.8
Q12	4.621	85.4	87.2	96.6
Q13	4.859	69.4	83.2	93.6
Q14	4.394	101.6	86.6	93.5
Q15	4.693	76.1	74.3	98.1
Q16	4.512	70.3	80.1	95.0
Q17	4.719	66.4	79.2	97.8
Q18	4.321	73.6	80.7	93.0
Q19	4.593	83.3	85.5	97.6
Median	4.708	69.2	76.2	96.6

Table 3. Table 3: Physical qubits used on IBM System One for state tomography and MQC experiments

$N$	Physical qubits used
2	[5, 10]
3	[5, 10, 6]
4	[5, 10, 6, 11]
5	[5, 10, 6, 11, 0]
11	[5, 10, 6, 11, 0, 12, 7, 15, 1, 8, 13]
12	[5, 10, 6, 11, 0, 12, 7, 15, 1, 8, 13, 16]
13	[5, 10, 6, 11, 0, 12, 7, 15, 1, 8, 13, 16, 2]
14	[5, 10, 6, 11, 0, 12, 7, 15, 1, 8, 13, 16, 2, 9]
15	[5, 10, 6, 11, 0, 12, 7, 15, 1, 8, 13, 16, 2, 9, 17]
16	[5, 10, 6, 11, 0, 12, 7, 15, 1, 8, 13, 16, 2, 9, 17, 4]
17	[5, 10, 6, 11, 0, 12, 7, 15, 1, 8, 13, 16, 2, 9, 17, 4, 14]
18	[5, 10, 6, 11, 0, 12, 7, 15, 1, 8, 13, 16, 2, 9, 17, 4, 14, 3]
19	[5, 10, 6, 11, 0, 12, 7, 15, 1, 8, 13, 16, 2, 9, 17, 4, 14, 3, 18]
20	[5, 10, 6, 11, 0, 12, 7, 15, 1, 8, 13, 16, 2, 9, 17, 4, 14, 3, 18, 19]

Equations

$$S_{\varphi }=|⟨\mathrm{000...00}|U_{\text{GHZ}}^{†}{U}_{\varphi }{U}_{\text{GHZ}}|\mathrm{000...00}⟩{|}^2=\text{Tr}({\rho }_{\varphi }\rho )$$

$$S_{\varphi }^{\text{ideal}}=\frac{1}{2}(1+\cos (N\varphi ))$$

$$I_q={\mathcal{N}}^{-1}|\sum _{\varphi }{e}^{iq\varphi }{S}_{\varphi }|$$

$$2\sqrt{I_N}\le F\le \sqrt{I_0/2}+\sqrt{I_N}$$

$$|0⟩$$

$${|0⟩}_0$$

$$|A{v}_{\text{cal}}-v_{\text{mea}}{|}^2$$

$$e^{-i\frac{\varphi }{2}{\sum }_j{\sigma }_z^j}{\rho }_q{e}^{i\frac{\varphi }{2}{\sum }_j{\sigma }_z^j}=e^{-iq\varphi }{\rho }_q,\text{commutator}{\sum }_j{\sigma }_z^j/2{\rho }_q=q{\rho }_q$$

$$S_{\varphi }$$

$$=\sum _q{e}^{-iq\varphi }\text{Tr}({\rho }_q{\rho }_{-q})=\sum _q{e}^{-iq\varphi }{I}_q$$

$${\rho }_0^{\text{GHZ}}$$

$${\rho }_N^{\text{GHZ}}$$

$${\rho }_{-N}^{\text{GHZ}}$$

$$I_0^{\text{GHZ}}=\text{Tr}({\rho }_0^{\text{GHZ}}{\rho }_0^{\text{GHZ}})=\frac{1}{2},I_N^{\text{GHZ}}=\text{Tr}({\rho }_N^{\text{GHZ}}{\rho }_{-N}^{\text{GHZ}})=\frac{1}{4},I_{-N}^{\text{GHZ}}=I_N^{\text{GHZ}}$$

$$2\sqrt{I_N}\le F\le \sqrt{I_0/2}+\sqrt{I_N}$$

$$\text{Tr}(\rho {\rho }^{\text{GHZ}})=\sum _{q=-N}^N\text{Tr}({\rho }_q{\rho }_{-q}^{\text{GHZ}})$$

$$\le \sum _{q=-N}^N\sqrt{\text{Tr}({\rho }_q{\rho }_{-q})\text{Tr}({\rho }_q^{\text{GHZ}}{\rho }_{-q}^{\text{GHZ}})}=\sqrt{I_0/2}+\sqrt{I_N}$$

$$F$$

$$\ge 2(\text{Tr}({\rho }_N{\rho }_{-N}^{\text{GHZ}})+\text{Tr}({\rho }_{-N}{\rho }_N^{\text{GHZ}}))$$

$$|{\psi }_j⟩={\alpha }_j|\mathrm{00...0000}⟩+{\beta }_j|\mathrm{11...1111}⟩+\cdots$$

$$\frac{1}{2}\sum _j{w}_j(|{\alpha }_j{|}^2+|{\beta }_j{|}^2)\ge \frac{1}{2}\sum _j{w}_j({\alpha }_j{\beta }_j^{\ast }+{\alpha }_j^{\ast }{\beta }_j)$$

$$\text{Tr}({\rho }_0{\rho }_0^{\text{GHZ}})\ge \text{Tr}({\rho }_N{\rho }_{-N}^{\text{GHZ}})+\text{Tr}({\rho }_{-N}{\rho }_N^{\text{GHZ}})$$

$$F=\frac{1}{2}(P_{\mathrm{000..00}}+P_{\mathrm{111..11}})+\sqrt{I_N}$$

Full text

Universal quantum computers promise to solve many problems that are intractable classically Shor (1996); Bravyi et al. (2018), but achieving fault tolerance will require a number of resources that are unavailable today. Until we can implement error correction, quantum systems will be beset with a certain amount of noise. Understanding how to best benchmark these near-term quantum devices is an active question Preskill (2018). Traditionally, the field has relied on local metrics such as one- and two-qubit gate fidelities since these are experimentally feasible even with full tomographic methods Knill et al. (2008); Magesan et al. (2012a, b); Sheldon et al. (2016a). However, it has become increasingly clear that such local metrics do not capture the full intricacies of a multi-qubit device. Therefore, a number of multi-qubit metrics such as direct fidelity estimation Flammia and Liu (2011); da Silva et al. (2011), three qubit simultaneous randomized benchmarking (RB) McKay et al. (2017a), direct RB Proctor et al. (2018), and quantum volume Cross et al. (2018) have been proposed and measured. Another powerful multi-qubit metric is that of entanglement, specifically, measuring the largest possible multipartite entangled state on a device Wang et al. (2018a); Mooney et al. (2019). Not only is the ability to generate entanglement indicative of high fidelity gate operations and qubit coherence, entangled states are the cornerstone of quantum speedups and they can be direct resources for quantum computing Shor (1997); Raussendorf and Briegel (2001). Multipartite entanglement in Greenberger-Horne-Zeilinger (GHZ) states have been demonstrated with 10 superconducting circuits Song et al. (2017), 14 trapped ions Monz et al. (2011), and 18 photons Wang et al. (2018b). Recently, multipartite entanglement in a 12-qubit linear graph state was verified in a superconducting qubit architecture Gong et al. (2019).

Here we generate and verify an 18-qubit entangled GHZ state on a 20-qubit superconducting device. Our entanglement metric is inspired by quantum sensing Degen et al. (2017) and can be used to directly bound the state fidelity. The device is comprised of 20 fixed frequency transmon qubits, and implements two-qubit gates based on cross-resonance driving Chow et al. (2011); Rigetti and Devoret (2010); Sheldon et al. (2016b). The device layout and the two qubit errors are shown in FIG. 1.

We verify the generation of GHZ states by measuring multiple quantum coherences (MQC) Baum et al. (1985), a tool traditionally used in solid state NMR and more recently in trapped ions to study many-body correlations and quantum information scrambling Wei et al. (2018); Gärttner et al. (2017). The experimental method to measure MQC has a strong overlap with quantum sensing and entanglement assisted metrology Cappellaro et al. (2005). In the prototypical quantum sensing circuit shown in FIG. 2A, a GHZ state is used to sense static magnetic fields with Heisenberg-limited sensitivity; it works by taking advantage of an ideal GHZ state’s amplified sensitivity to phase rotations of each of the individual qubits in the entangled state Bollinger et al. (1996); Leibfried et al. (2004); Giovannetti et al. (2004); Degen et al. (2017). If each qubit has a phase rotation of $\phi$, then the $N$-qubit GHZ state rotates collectively by $N\phi$. By observing how sensitive a nonideal GHZ state responds to rotations, we can deduce how much entanglement is present in the state. The quantum circuit for measuring MQC is illustrated in FIG. 2B, and it can be described in four steps:

Starting from the $N$-qubit ground state: $\ket{\text{GS}}=\ket{00..000}$, apply a Hadamard gate on qubit 0 followed by a sequence of CX gates. Ideally this brings the system into the GHZ state: $\ket{\text{GHZ}}=\frac{1}{\sqrt{2}}(\ket{000..00}+\ket{111..11})$ 2. 2.

Apply a collective rotation given by the unitary $U_{\phi}$ on all qubits. This amounts to adding a phase $N\phi$ to the GHZ state: $\frac{1}{\sqrt{2}}(\ket{000..00}+e^{-iN\phi}\ket{111..11})$ 3. 3.

Disentangle the GHZ state by performing the CX gate sequence in reverse order. The amplified phase is mapped onto qubit 0: $\frac{1}{\sqrt{2}}(\ket{0}+e^{-iN\phi}\ket{1})\otimes\ket{00..00}$ 4. 4.

Read out the amplified phase by measuring the probability of the system returning to its initial state: $\ket{\text{GS}}$

The measured signal of this protocol is given by

[TABLE]

where $\rho=U_{\text{GHZ}}\outerproduct{\text{GS}}{\text{GS}}U_{\text{GHZ}}^{\dagger}$, $U_{\text{GHZ}}=U_{\text{CX}}H_{0}$, and $\rho_{\phi}=U_{\phi}\rho U_{\phi}^{\dagger}$. If our controls are perfect and there is no decoherence, Eq. (1) reduces to

[TABLE]

which can also be obtained by measuring the state $\ket{0}$ on qubit 0 in the final step. Since $S_{\phi}$ comes from the overlap between a rotated and unrotated density matrix, in the final step of the protocol the probability of all qubits being in the zero state must be measured. The constant term in $S^{\text{ideal}}$ comes from the diagonal elements of the GHZ density matrix, whereas the oscillating term comes from the off-diagonal corner elements. Any difference between $S_{\phi}$ and $S_{\phi}^{\text{ideal}}$ is an indication that our GHZ state is imperfect. To quantify the state fidelity we focus on the MQC amplitudes, defined as the discrete Fourier transform of $S_{\phi}$:

[TABLE]

where $\mathcal{N}$ is a normalization factor. The $N$-qubit GHZ state fidelity defined as $F=\matrixelement{\text{GHZ}}{\rho}{\text{GHZ}}$ can be bounded by

[TABLE]

For a perfect GHZ state we have $I_{0}=2I_{N}=1/2$, and all other $I_{q}$ being zero. We can also directly obtain the state fidelity as $F=\frac{1}{2}(P_{000..00}+P_{111..11})+\sqrt{I_{N}}$, where $P_{000..00}$ and $P_{111..11}$ are the populations of $\ket{000..00}$ and $\ket{111..11}$ in the density matrix. A discussion on MQC amplitudes and proof of Eq. (2) are given in SM. For a $N$-qubit state to have multipartite GHZ entanglement, it needs to have a minimal fidelity of 0.5 Gühne and Tóth (2009); Gühne and Seevinck (2010).

While this method departs from parity oscillation measurements commonly used in trapped ions to verify GHZ entanglement Sackett et al. (2000); Leibfried et al. (2005); Monz et al. (2011), MQC offers two main benefits: robustness to noise and scalability in readout correction. Parity oscillations measure the coherence between the $\ket{000..00}$ and $\ket{111..11}$ states. It works by looking at the oscillations in the expectation value $\expectationvalue{ZZZ...ZZ}{\text{GHZ}_{\phi}}$ as a function of $\phi$, where $\text{GHZ}_{\phi}=\otimes^{N}_{j}e^{\frac{i\pi}{4}(\cos\phi\sigma_{x}^{j}+\sin\phi\sigma_{y}^{j})}\ket{\text{GHZ}}$ Monz et al. (2011). The amplitude of the parity oscillations gives coherence $C=|\rho_{000..00,111..11}|+|\rho_{111..11,000..00}|$, which is related to the fidelity via $F=\frac{1}{2}(P_{000..00}+P_{111..11}+C)$ Monz et al. (2011). The coherence is related to MQC amplitudes via $C=2\sqrt{I_{N}}$. At first our entanglement metric appears disadvantageous compared to parity oscillations since it takes twice the circuit length. However our experiments can be made robust against noise. Just as a Hahn echo refocuses low frequency noise and reduces dephasing Hahn (1950), adding a $\pi$-pulse after making the GHZ state can dramatically improve the measured fidelity. The quantum circuit for refocused MQC is illustrated in FIG. 2C. We find experimentally the 20-qubit state fidelity to increase by nearly $11\%$ by adding the refocusing $\pi$-pulse, as shown in FIG. 4B. In addition, our entanglement metric only requires measurement of the initial state.

In addition to accommodating dynamical decoupling techniques, the MQC method is also less sensitive to readout errors. We point out that parity oscillations in GHZ states have been measured previously on the IBM Q 16-qubit device Cruz et al. (2018) with average readout error of $7\%$, but multipartite entanglement cannot be established beyond five qubits. Aside from control imperfections and decoherence, readout errors limit our ability to measure the entangled state, even though the state itself can be highly entangled. Since readout errors are independent from entanglement we can calibrate them out of the measurement. To mitigate measurement errors, we experimentally construct a $2^{N}$ by $2^{N}$ calibration matrix, $A$, where each row vector corresponds to the measured outcome probabilities of a prepared basis state. In the nominal case of no readout error, $A$ is an identity matrix. With readout error, we can correct the measured counts $v_{\text{mea}}$ by minimizing

[TABLE]

under the constraint $\sum_{j}v_{\text{cal},j}=1$ and $v_{\text{cal},j}\geq 0$. Here $v_{\text{cal}}$ is the calibrated counts of $v_{\text{mea}}$. Eq. (3) can be recasted into a convex optimization problem and solved by quadratic programming using packages such as CVXOPT Anderson et al. (2019). The overhead for measuring $A$ and minimizing Eq. (3) increases exponentially with $N$. We modify this calibration procedure to have a scalable way to perform readout correction using two key features of the MQC method. For one, we only need to measure the probability that the sate is in $\ket{000..00}$. Two, we expect that the MQC output for imperfections in the GHZ state to result in low excitation states. This is not entirely unexpected since for an ideal GHZ state, the MQC output has only two distinct states: $\ket{000..00}$ and $\ket{100..00}$. With an imperfect GHZ state and readout errors we expect the output counts to spread out but stay within the low excitation manifolds. Combining these two features we can significantly reduce the overhead for readout calibration by truncating $A$ to only correct for states with significant weights. For example, we can reasonably conclude that the final measurement will not include states such as $\ket{111..11}$. We justify this truncation in SM, where we present the corrected data as a function of the number of states kept in the calibration. A rapid convergence is observed with just 256 states. While keeping such a low number of states may not be sufficient to correct for the entire output vector, it is sufficient to accurately correct for the $\ket{000..00}$ state. In addition, the largest 256 states all have similar excitation distributions centered around three excitations independent of $N$, as shown in SM. This suggests the readout calibration based on the truncated $A$, call it $A_{t}$, is scalable for measuring MQC amplitudes in GHZ states. The $A_{t}$ matrix cannot correct readout errors from parity oscillation experiments, since the output counts will be distributed across all eigenstates. An alternative to scalable readout calibration is to approximate the full $A$ matrix as a tensor product of $A$ matrices of each qubit Song et al. (2017); Gong et al. (2019), this approach is valid when there is little to no readout cross-talk between qubits.

We have experimentally generated GHZ states for $N=11$ to $N=20$, and measured $S_{\phi}$ in Eq. (1) and extracted the corresponding MQC amplitudes. The data for 12, 14, 16, and 18-qubit GHZ states are shown in FIG. 3, and the circuit used in the 18-qubit experiment is shown in FIG. 2D. For each $N$, we measure $S_{\phi}$ for $\phi=\frac{\pi j}{N+1}$, where $j=0,1,2,\cdots,2N+1$ so the highest frequency detectable is $N+1$. The result of each experiment is averaged over 16384 shots, and the errorbar corresponds to one standard error obtained from eight experiments. There is considerable difference between the results with and without readout correction. Here $A_{t}$ is constructed using 256 basis states, and each basis state measurement is averaged over 4096 shots.

From the experimentally extracted MQC amplitudes shown in the bottom row of FIG. 3, we see one peak located at $q=0$ and two peaks at $q=\pm N$, characteristic of $N$-qubit GHZ states. Peak amplitudes become lower with increasing $N$, indicating larger $N$-qubit GHZ states have lower fidelities. Using Eq. (2) we extract the upper and lower bounds on state fidelities with readout calibrations as a function of $N$, as shown in FIG. 4A. For $N=11$ to $N=17$ the fidelity lower bound is clearly higher than the 0.5 threshold for multipartite entanglement. For $N=18$ the lower bound is $0.5006\pm 0.0067$, in this case we measure $P_{000..00}$ and $P_{111..11}$ for the GHZ state in addition to MQC amplitudes to obtain the state fidelity of $F=0.5165\pm 0.0036$, confirming that the 18-qubit GHZ state is multipartite entangled. We have not been able to establish multipartite entanglement with 19 and 20-qubit GHZ states. Without applying the aforementioned readout calibration, the highest number of multipartite entangled qubits we can measure is $N=14$ with a fidelity lower bound of $0.5406\pm 0.0037$. We compare fidelities extracted from our method with that from quantum state tomography (QST) Smolin et al. (2012) for small GHZ states, the results are summarized in Table 1. While the MQC method appears to give slightly higher state fidelity, we expect these results are within the errors of the tomography experiments. Each experiment is averaged over 16384 shots and readout corrected. Our method to experimentally quantify multipartite entanglement for GHZ states can be applied to other states that are locally equivalent to GHZ states, such as star graph and complete graph states Hein et al. (2004). The only difference is in the rotation step. For star graph states, instead of applying $U_{\phi}$ on all qubits, apply $U_{\phi}$ on the central qubit and $HU_{\phi}H$ on the rest. For complete graph states, apply $e^{-\frac{i\pi}{4}\sigma_{x}}U_{\phi}e^{\frac{i\pi}{4}\sigma_{x}}$ on all qubits.

There are several experimental limitations to entangling large GHZ states in our device. First, the circuit depth required to generate a $N$-qubit GHZ state scales as $\mathcal{O}(N)$; in contrast with the linear graph state, where only two steps are needed independent of $N$ Gong et al. (2019); Mooney et al. (2019). This makes GHZ states particularly fragile to decoherence. We choose an entangling path on our device that takes the least amount of time to complete; the physical qubits involved for each $N$ are listed in SM. Second, there are instances where CX gates, implemented by cross-resonance driving, are run on adjacent qubit pairs. Simultaneous adjacent CX gates will have lower fidelity than individual CX gates due to always-on $ZZ$ interactions and cross-driving effects between neighboring qubits. Third, since the CX gates are applied sequentially, there will inevitably be free evolution on idle qubits leading to unitary errors. In addition, due to pulse alignment restrictions in the software, the entangling and the disentangling operations take different times to complete, making the $\pi$-pulse not as effective as it can be. This might explain why the measured $S_{\phi}$ appears phase-shifted.

We use a basic noise model built from device parameters to simulate the 20-qubit MQC experiments using Qiskit Aleksandrowicz et al. (2019). The simulation models one- and two-qubit gate errors as a depolarizing error plus a thermal relaxation error such that the total error equals the error measured experimentally from randomized benchmarking; the details of the basic noise model can be found in Qiskit tutorial Aleksandrowicz et al. (2019). To compare with experiments, we turn off readout error in the simulation and average over 2048 shots. The simulation shows higher fidelities than the experiment and does not appear to capture the effects of the refocusing $\pi$ pulse, as shown in FIG. 4C. This suggests that the experimental system has slow drifts which can be refocused by $\pi$-pulses. Interestingly, it has been demonstrated that dynamical decoupling is remarkably effective at extending the lifetime of GHZ states Kaufmann et al. (2017).

We demonstrate in this work an experimentally scalable entanglement metric based on multiple quantum coherences, and applied it to verify 18-qubit multipartite GHZ entanglement. Our experiments show encouraging results in the ability to entangle and disentangle highly correlated many-body states in near term quantum devices. We are exploring new variations of CX gates which can be applied simultaneously on adjacent qubit pairs while canceling $ZZ$ errors Takita et al. (2017). This should improve gate fidelity in the entangling and disentangling steps. The lifetimes of MQC amplitudes should be measured and compared to those of parity oscillations, which were reported to decrease as $N^{2}$ in trapped ions Monz et al. (2011) and $N$ in superconducting qubits Ozaeta and McMahon (2019). It will be interesting to extend MQC to other entangled states such as the W-state and study their entanglement properties. Lastly, the newly developed error mitigation techniques Temme et al. (2017); Kandala et al. (2019) may give us insights to the maximum GHZ fidelity achievable in our device in the limit of zero noise.

During the preparation of this manuscript, we became aware of recent experiments demonstrating 18-qubit multipartite GHZ entanglement in a tunable-frequency transmon device Song et al. (2019) and 20-qubit multipartite GHZ entanglement in a Rydberg atoms array Omran et al. (2019).

SUPPLEMENTARY MATERIAL

.1 Multiple Quantum Coherences

Consider writing the density matrix as $\rho=\sum_{m,m^{\prime}}\rho_{m,m^{\prime}}\outerproduct{m}{m^{\prime}}$, where the basis states satisfy $\sum_{j}\sigma_{z}^{j}/2\ket{m}=m\ket{m}$. We can expand the density matrix as $\rho=\sum_{q}\rho_{q}$, where $\rho_{q}=\sum_{m}\rho_{m,m-q}\outerproduct{m}{m-q}$. It can be shown that $\rho_{q}$ satisfies the following:

[TABLE]

Since $\rho$ is hermitian we also have $\rho_{q}^{\dagger}=\rho_{-q}$. Each $\rho_{q}$ occupies a different part of the density matrix, and it obeys the orthogonality condition $\text{Tr}(\rho_{q}\rho_{p})=\delta_{q,-p}\text{Tr}(\rho_{q}\rho_{-q})$. While each $\rho_{q}$ is not directly observable, the trace $I_{q}=\text{Tr}(\rho_{q}\rho_{-q})$ is. $I_{q}$ is the multiple quantum coherence amplitude, and it can be found by Fourier transforming the overlap signal $S_{\phi}=\text{Tr}(\rho_{\phi}\rho)$, where $\rho_{\phi}=e^{-i\frac{\phi}{2}\sum_{j}\sigma_{z}^{j}}\rho e^{i\frac{\phi}{2}\sum_{j}\sigma_{z}^{j}}$. For a general time-dependent density matrix, measuring $S_{\phi}$ requires the ability to implement time-reversed evolution. Upon expanding $\rho$ inside $S_{\phi}$ and using the first relation in Eq. (4) we have

[TABLE]

where in the second to last step we used the orthogonality condition. Fourier transforming gives $I_{q}=\mathcal{N}^{-1}\sum_{\phi}e^{iq\phi}S_{\phi}$, where $\mathcal{N}$ is a normalization constant depending on the number of $\phi$ used in the experiments. If the maximum coherence order to be measured is $q_{\text{max}}$, we need at least $2q_{\text{max}}$ experiments. The angle $\phi$ can be chosen as $\phi=\frac{\pi j}{q_{\text{max}}}$ where $j=0,1,2,\cdots 2q_{\text{max}}-1$ and $\mathcal{N}=2q_{\text{max}}$. While our discussion assumes unitary evolution, it has been shown that in certain types of decoherence the experimental method for measuring multiple quantum coherence is still valid Gärttner et al. (2018).

For an ideal $N$-qubit GHZ state, the nonzero elements in the density matrix resides only in the four corners. Therefore only three components arise in the expansion: $\rho^{\text{GHZ}}=\rho_{0}^{\text{GHZ}}+\rho_{N}^{\text{GHZ}}+\rho_{-N}^{\text{GHZ}}$. Explicitly they are given by

[TABLE]

the corresponding multiple quantum amplitudes are given by

[TABLE]

Multiple quantum coherence amplitudes are symmetric: $I_{q}=I_{-q}$.

.2 Fidelity Bounds from MQC amplitudes

The state fidelity, given by $F=\matrixelement{\text{GHZ}}{\rho}{\text{GHZ}}=\text{Tr}(\rho\rho^{\text{GHZ}})$, can be bounded by:

[TABLE]

The upper bound on $F$ follows from the Cauchy-Schwarz inequality:

[TABLE]

To prove the lower bound we first notice that

[TABLE]

This can be proved by writing the density matrix as $\rho=\sum_{j}w_{j}\outerproduct{\psi_{j}}{\psi_{j}}$, where $w_{j}\geq 0$ and $\sum_{j}w_{j}=1$. The state vectors $\ket{\psi_{j}}$ need not to be orthogonal, we can in general expand $\ket{\psi_{j}}$ as

[TABLE]

Since $|\alpha_{j}-\beta_{j}|^{2}\geq 0$, upon expanding we have $|\alpha_{j}|^{2}+|\beta_{j}|^{2}\geq\alpha_{j}\beta_{j}^{*}+\alpha_{j}^{*}\beta_{j}$. It then follows that

[TABLE]

which is the same as

[TABLE]

thereby proving Eq. (.2). One can go one step further by noting that $\rho_{N}=\kappa\rho_{N}^{\text{GHZ}}$, where $\kappa$ is a complex constant. Using $\text{Tr}(\rho_{N}\rho_{-N}^{\text{GHZ}})=\kappa I_{N}^{\text{GHZ}}=\frac{1}{\kappa^{*}}I_{N}$ we can show that $|\kappa|=2\sqrt{I_{N}}$. Notice $\kappa$ can always be made real by appropriately rotating the density matrix $\rho$. Substituting $\rho_{N}=2\sqrt{I_{N}}\rho_{N}^{\text{GHZ}}$ into Eq. (.2) gives $2\sqrt{I_{N}}$ as the lower bound of $F$. In addition to the fidelity bounds, we can also obtain the GHZ state fidelity $F$:

[TABLE]

where $P_{000..00}=\expectationvalue{\rho}{000...00}$ and $P_{111..11}=\expectationvalue{\rho}{111...11}$ are the probabilities of finding all zeroes and all ones in the state $\rho$.

.3 Readout Calibration

We construct a truncated readout calibration matrix $A_{t}$ according to the largest 256 states by weight in the output of each experiments. In this section we justify using only 256 states as opposed to $2^{N}$ states. Our metric to verify GHZ entanglement requires MQC amplitudes $I_{0}$ and $I_{N}$. In FIG. 5A we plot the corrected $I_{0}$ and $I_{N}$ as a function of the number of states used in $A_{t}$; we see a rapid convergence after just 32 states. Interestingly, $I_{0}$ decreases as we add more states into the readout calibration, while $I_{N}$ is relatively unchanged.

In FIG. 5B we compare all counts and the largest 256 counts from all experiments for each $N$. The counts are grouped according to excitation number (number of ones) and divided by the total number of counts. For small $N$, there is little difference between all counts and the largest 256 counts, indicating most of the weight in the output are contained in the largest 256 counts. As $N$ increases however, the total counts begins to spread out to higher excitation numbers, while the largest 256 counts are still localized in the low excitation numbers. This discrepancy however does not affect the calibrated values of the all zeroes count, as demonstrated by the convergence of MQC amplitudes.

In FIG. 5C we examine the distribution of the largest 245 counts as a function of excitation number for each $N$. Similar distributions centered around three excitation are observed for all $N$.

.4 Device Parameters and Qubits involved in making GHZ states

In Table 2 we show the typical qubit parameters for the device. In Table 3 we list the physical qubits involved in the state tomography experiments and the $N$-qubit MQC experiments on the device. The number labeling the physical qubit is shown in FIG. 1 in the main text.