10-qubit entanglement and parallel logic operations with a superconducting circuit

TL;DR

This paper demonstrates the creation of 10-qubit entangled states and parallel quantum operations in a superconducting circuit, marking the largest entanglement in solid-state systems and advancing scalable quantum computing.

Contribution

It reports the first production and tomography of 10-qubit GHZ states in a superconducting circuit with parallel operations and high fidelity.

Findings

Achieved 10-qubit GHZ state with 66.8% fidelity

Demonstrated parallel qubit-qubit interactions

Largest entanglement in solid-state quantum systems

Abstract

Here we report on the production and tomography of genuinely entangled Greenberger-Horne-Zeilinger states with up to 10 qubits connecting to a bus resonator in a superconducting circuit, where the resonator-mediated qubit-qubit interactions are used to controllably entangle multiple qubits and to operate on different pairs of qubits in parallel. The resulting 10-qubit density matrix is unambiguously probed, with a fidelity of $0.668 \pm 0.025$. Our results demonstrate the largest entanglement created so far in solid-state architectures, and pave the way to large-scale quantum computation.

Tables

Table 1. Table S1: Qubit characteristics. ω j 0 superscript subscript 𝜔 𝑗 0 \omega_{j}^{0} is the maximum resonant frequency of Q j subscript 𝑄 𝑗 Q_{j} , i.e., the resonant frequency at the qubit’s sweetpoint where dephasing is minimum. ω I subscript 𝜔 𝐼 \omega_{I} is the interaction frequency for qubits to entangle, which is about 2 ​ π × 140 2 𝜋 140 2\pi\times 140  MHz below the resonator frequency ω B subscript 𝜔 𝐵 \omega_{B} . Qubit coherence parameters including the lifetime T 1 , j subscript 𝑇 1 𝑗 T_{1,j} , the Ramsey Gaussian dephasing time T 2 , j ∗ superscript subscript 𝑇 2 𝑗 ∗ T_{2,j}^{\ast} , and the spin echo Gaussian dephasing time T 2 , j SE superscript subscript 𝑇 2 𝑗 SE T_{2,j}^{\textrm{SE}} of Q j subscript 𝑄 𝑗 Q_{j} are measured at ω I subscript 𝜔 𝐼 \omega_{I} where other qubits are far detuned. ω j subscript 𝜔 𝑗 \omega_{j} is the idle frequency of Q j subscript 𝑄 𝑗 Q_{j} , where the X / 2 𝑋 2 X/2 gate for state preparation and the tomographic pulses for qubit measurement are applied. g j subscript 𝑔 𝑗 g_{j} is the coupling strength between Q j subscript 𝑄 𝑗 Q_{j} and the resonator B 𝐵 B , defined in the interaction Hamiltonian ℏ ​ g j ​ ( σ j + + σ j − ) ​ ( a † + a ) Planck-constant-over-2-pi subscript 𝑔 𝑗 superscript subscript 𝜎 𝑗 superscript subscript 𝜎 𝑗 superscript 𝑎 † 𝑎 \hbar g_{j}(\sigma_{j}^{+}+\sigma_{j}^{-})(a^{\dagger}+a) , where σ j + superscript subscript 𝜎 𝑗 \sigma_{j}^{+} ( σ j − superscript subscript 𝜎 𝑗 \sigma_{j}^{-} ) is the raising (lowering) operator of Q j subscript 𝑄 𝑗 Q_{j} and a † superscript 𝑎 † a^{\dagger} and a 𝑎 a are field operators of resonator B 𝐵 B . Fidelity of the X / 2 𝑋 2 X/2 gate on Q j subscript 𝑄 𝑗 Q_{j} is characterized by RB at its idle frequency ω j subscript 𝜔 𝑗 \omega_{j} ; we also simultaneously run RBs on 8 qubits to characterize the X / 2 𝑋 2 X/2 gate fidelities for these 8 qubits. Each qubit is measured through its own readout resonator that connects to a common circumferential transmission line, and ω j r superscript subscript 𝜔 𝑗 𝑟 \omega_{j}^{r} is the resonant frequency of Q j subscript 𝑄 𝑗 Q_{j} ’s readout resonator R j subscript 𝑅 𝑗 R_{j} , with the Q j subscript 𝑄 𝑗 Q_{j} - R j subscript 𝑅 𝑗 R_{j} coupling strength noted as g j r superscript subscript 𝑔 𝑗 𝑟 g_{j}^{r} . δ ​ ω j m 𝛿 superscript subscript 𝜔 𝑗 𝑚 \delta\omega_{j}^{m} is the AC stark shift of Q j subscript 𝑄 𝑗 Q_{j} during its measurement, i.e., the qubit frequency shifts downwards from ω j subscript 𝜔 𝑗 \omega_{j} as its readout resonator R j subscript 𝑅 𝑗 R_{j} becomes populated. The readout pulse is 1  μ 𝜇 \mu s-long, and carries 10 tones targeting the 10 readout resonators. n j r superscript subscript 𝑛 𝑗 𝑟 n_{j}^{r} ( = | δ ​ ω j m / 2 ​ χ j | absent 𝛿 superscript subscript 𝜔 𝑗 𝑚 2 subscript 𝜒 𝑗 =|\delta\omega_{j}^{m}/2\chi_{j}| , where χ j subscript 𝜒 𝑗 \chi_{j} is the dispersive frequency shift of the readout resonator when qubit changes state from | 0 ⟩ ket 0 |0\rangle to | 1 ⟩ ket 1 |1\rangle ) is the maximum average photon number in Q j subscript 𝑄 𝑗 Q_{j} ’s readout resonator R j subscript 𝑅 𝑗 R_{j} as excited by the readout pulse for qubit measurement, which decays at the end of the readout pulse with a leakage rate of κ j r superscript subscript 𝜅 𝑗 𝑟 \kappa_{j}^{r} for R j subscript 𝑅 𝑗 R_{j} returning to the vacuum state and being ready for the next cycle. F 0 , j subscript 𝐹 0 𝑗 F_{0,j} ( F 1 , j subscript 𝐹 1 𝑗 F_{1,j} ) is the probability of correctly reading out Q j subscript 𝑄 𝑗 Q_{j} in | 0 ⟩ ket 0 |0\rangle ( | 1 ⟩ ket 1 |1\rangle ) when it is prepared in | 0 ⟩ ket 0 |0\rangle ( | 1 ⟩ ket 1 |1\rangle ), which are used to correct the measured qubit probabilities  [ 43 ] . Q 3 subscript 𝑄 3 Q_{3} ’s readout fidelities apply to the entangling experiments with N ≤ 9 𝑁 9 N\leq 9 . We use two concatenated X / 2 𝑋 2 X/2 gates to prepare the qubit in | 1 ⟩ ket 1 |1\rangle . Values in parenthesis, if available, are one standard deviation.

		$Q_1$	$Q_2$	$Q_3$	$Q_4$	$Q_5$	$Q_6$	$Q_7$	$Q_8$	$Q_9$	$Q_{10}$
${\omega }_j^0/2\pi$ (GHz)		$5.782$	$5.831$	$5.828$	$5.780$	$5.760$	$5.863$	$5.780$	6.004	5.893	5.930
${\omega }_I/2\pi$ (GHz)		$\approx 5.655$
$T_{1,j}$ ($\mu$s)		27.2	24.4	10.9	15.0	19.2	23.7	13.8	11.8	17.1	22.0
$T_{2,j}^{\ast }$ ($\mu$s)		2.9	2.8	2.8	2.2	2.6	1.8	1.1	2.1	1.7	4.4
$T_{2,j}^{\text{SE}}$ ($\mu$s)		11.8	10.6	10.0	10.8	11.7	8.9	8.0	8.0	7.9	11.8
${\omega }_j/2\pi$ (GHz)		5.080	5.467	5.657	5.042	5.179	5.605	4.960	5.260	5.146	5.560
$g_j/2\pi$ (MHz)		14.2	20.5	19.9	20.2	15.2	19.9	19.6	18.9	19.8	16.3
$X/2$ fidelity		0.9985(2)	0.9992(1)	0.9984(1)	0.9987(2)	0.9991(1)	0.9964(5)	0.9987(1)	0.9980(3)	0.9988(3)	0.9989(1)
Simultaneous
$X/2$ fidelity		0.9978(7)	0.9980(2)	0.9953(5)	0.9955(5)	0.9985(2)	0.9962(3)	–	–	0.9979(2)	0.9931(12)
${\omega }_j^r/2\pi$ (GHz)		6.509	6.541	6.615	6.614	6.635	6.694	6.691	6.794	6.809	6.891
$g_j^r/2\pi$ (MHz)		41.3	39.9	40.6	38.2	38.5	40.4	41.8	40.9	40.2	38.7
$\delta {\omega }_j^m/2\pi$ (MHz)		31.1	32.7	21.1	46.5	9.0	45.1	22.5	19.5	26.0	70.2
$n_j^r$		92	59	31	180	30	81	93	74	103	203
$1/{\kappa }_j^r$ (ns)		291	275	272	348	223	284	248	266	299	242
$F_{0,j}$		0.921	0.955	0.982	0.974	0.962	0.988	0.950	0.970	0.961	0.971
$F_{1,j}$		0.867	0.915	0.904	0.928	0.927	0.917	0.922	0.880	0.894	0.934

Table 2. Table S2: Nearest-neighbor X ​ X 𝑋 𝑋 XX -type crosstalk couplings.

		$Q_1$-$Q_2$	$Q_2$-$Q_3$	$Q_3$-$Q_4$	$Q_4$-$Q_5$	$Q_5$-$Q_6$
${\lambda }_{j,j^{\prime }}^c/2\pi$ (MHz)		1.7	2.6	2.3	2.2	0.2
		$Q_6$-$Q_7$	$Q_7$-$Q_8$	$Q_8$-$Q_9$	$Q_9$-$Q_{10}$	$Q_{10}$-$Q_1$
${\lambda }_{j,j^{\prime }}^c/2\pi$ (MHz)		2.2	2.3	2.1	2.1	0.06

Equations

$$H/\hslash ={\omega }_B{a}^{+}a+\sum _{j=1}^{10}\left[{\omega }_j|1_j⟩⟨1_j|+g_j({\sigma }_j^{+}a+{\sigma }_j^{-}{a}^{+})\right],$$

$$\begin{gathered} H_1/\hslash =\sum _{j\in \{k,l,m\}}{\lambda }_j\left({\sigma }_j^{-}{\sigma }_{j^{\prime }}^{+}+{\sigma }_j^{+}{\sigma }_{j^{\prime }}^{-}\right) \\ +\sum _{j\in \{k,l,m\}}\left[\frac{g_j^2}{{\Delta }_j}|1_j⟩⟨1_j|+\frac{g_{j^{\prime }}^2}{{\Delta }_j}|1_{j^{\prime }}⟩⟨1_{j^{\prime }}|\right], \end{gathered}$$

$$\left[\begin{array}{cccccccccc}1& -0.023& -0.001& 0.008& 0.028& 0.020& 0.009& 0.011& 0.003& -0.008\\ -0.018& 1& 0.050& 0.024& 0.008& 0.003& 0.000& -0.002& 0.005& -0.016\\ -0.021& -0.081& 1& 0.054& 0.022& 0.012& 0.003& 0.001& -0.004& -0.017\\ 0.017& 0.055& 0.080& 1& -0.024& -0.012& -0.003& 0& 0.003& 0.013\\ 0.016& 0.019& 0.009& -0.003& 1& -0.015& -0.002& 0.005& 0.009& 0.014\\ 0.001& -0.002& -0.004& -0.005& -0.025& 1& 0.009& 0.022& 0.013& 0.004\\ 0.001& 0& -0.003& -0.006& -0.028& -0.046& 1& 0.078& 0.035& 0.008\\ -0.004& 0& 0.001& 0.003& 0.013& 0.020& 0.025& 1& -0.029& -0.006\\ -0.012& -0.006& -0.002& 0& 0.008& 0.015& 0.016& 0.065& 1& -0.014\\ 0.002& 0.011& 0.015& 0.015& 0.029& 0.023& 0.010& 0.008& -0.011& 1\end{array}\right].$$

$$P_k^j=⟨k|U^j\rho {\left(U^j\right)}^{†}|k⟩=\sum _{l,m=1}^{2^N}{U}_{kl}^j{\left(U_{km}^j\right)}^{\ast }{\rho }_{lm},$$

Full text

Entanglement is one of the most counter-intuitive features of quantum mechanics. The creation of increasingly large number of maximally entangled quantum bits (qubits) is central for measurement-based quantum computation [1], quantum error correction [2, 3], quantum simulation [4], and foundational studies of nonlocality [5, 6] and quantum-to-classical transition [7]. A significant experimental challenge for engineering multiqubit entanglement [8, 9, 10] has been noise control [11, 12]. With solid-state platforms, the largest number of entangled qubits reported so far is five [10], and further scaling up would be difficult as constrained by the qubit coherence and the employed sequential-gate method.

Superconducting circuits are a promising solid-state platform for quantum state manipulation and quantum computing owing to the microfabrication technology scalability, individual qubit addressability, and ever-increasing qubit coherence time [13]. The past decade has witnessed significant progresses in quantum information processing and entanglement engineering with superconducting qubits: preparation of three- and four-qubit entangled states [14, 15, 16, 17], demonstration of elementary quantum algorithms [18, 19], realization of three-qubit Toffoli gates and quantum error correction [20, 24, 23, 22, 21]. In particular, a recent experiment has achieved a two-qubit controlled-phase gate with a fidelity above 99 percent with a superconducting quantum processor [10], where five transmon qubits with nearest-neighbor coupling are arranged in a linear array. Based on this gate, a 5-qubit Greenberger-Horne-Zeilinger (GHZ) state was produced step by step; the number of entangled qubits is increased by one at a time. With a similar architecture consisting of 9 qubits, digitized Trotter steps were used to emulate the adiabatic change of the system Hamiltonian that encodes a computational problem [25], where the digital evolution into a GHZ state with a fidelity of 0.55 was demonstrated for a 4-qubit system.

In this letter we demonstrate a versatile superconducting quantum processor featuring high connectivity with programmable qubit-qubit couplings mediated by a bus resonator, and experimentally produce GHZ states with up to 10 qubits using this quantum processor. The resonator-induced qubit-qubit couplings result in a phase shift that is quadratically proportional to the total qubit excitation number, evolving the participating qubits from an initially product state to the GHZ state after a single collective interaction, irrespective of the number of the entangled qubits [26]. We characterize the multipartite entanglement by quantum state tomography achieved by synchronized local manipulations and detections of the entangled qubits, and measure a fidelity of $0.668\pm 0.025$ for the 10-qubit GHZ state, which confirms the genuine tenpartite entanglement [27] with 6.7 standard deviations ($\sigma$). We also implement parallel entangling operations mediated by the resonator, simultaneously generating three Einstein-Podolsky-Rosen (EPR) pairs; this feature was previously suggested in the context of ion traps [28] and quantum dots coupled to an optical cavity [29], but experimental demonstrations are still lacking.

The superconducting quantum processor is illustrated in Fig. 1(a), which is constructed as 10 transmon qubits ($Q_{j}$ for $j$ = 1 to 10), with resonant frequencies $\omega_{j}/2\pi$ tunable from 5 to 6 GHz, symmetrically coupled to a central resonator ($B$), whose resonant frequency is fixed at $\omega_{B}/2\pi\approx$ 5.795 GHz. Measured qubit-resonator ($Q_{j}$-$B$) coupling strengths $g_{j}/2\pi$ range from 14 to 20 MHz (see Supplemental Material [30] for details on device, operation, and readout) [42]. The central resonator serves as a multipurpose actuator, enabling controlled long-range logic operations, scalable multiqubit entanglement, and quantum state transfer. In the rotating-wave approximation and ignoring the crosstalks between qubits (see Supplemental Material [30]), the Hamiltonian of the system is given by

[TABLE]

where $\sigma_{j}^{+}$ ($\sigma_{j}^{-}$) is the raising (lowering) operator of $Q_{j}$ and $a^{+}$ ($a$) is the creation (annihilation) operator of $B$.

The qubit-qubit coupling can be realized through the superexchange (SE) interaction [44] mediated by the bus resonator $B$ [45, 46, 47, 48]. With multiplexing we can further arrange multiple qubit pairs at different frequencies to turn on the intra-pair SE interactions simultaneously. To illustrate this feature, we consider three qubit pairs, $Q_{k}$-$Q_{k^{\prime}}$, $Q_{l}$-$Q_{l^{\prime}}$, and $Q_{m}$-$Q_{m^{\prime}}$, detuned from resonator $B$ by $\Delta_{j}$ ($\equiv\omega_{j}-\omega_{B}$, and $\omega_{j}=\omega_{j^{\prime}}$) for $j=k$, $l$, and $m$, respectively, while all other qubits are far detuned and can be neglected for now. In the dispersive regime and when the resonator $B$ is initially in the ground state, it will remain so throughout the procedure and the effective Hamiltonian for the qubit pairs is

[TABLE]

where $\lambda_{j}=\frac{g_{j}g_{j^{\prime}}}{\Delta_{j}}$, $|\Delta_{j}|\gg g_{j},\,g_{j}^{\prime}$, and $|\Delta_{j_{1}}-\Delta_{j_{2}}|\gg\lambda_{j_{1}},\,\lambda_{j_{1}}^{\prime},\,\lambda_{j_{2}},\,\lambda_{j_{2}}^{\prime}$ for $j_{1},\,j_{2}\in\{k,l,m\}$ and $j_{1}\neq j_{2}$. With this setting, the resonator $B$ is simultaneously used for three intra-pair SE processes; the inter-pair couplings are effectively switched off due to large detunings between different pairs.

With the fast Z control on each qubit, coupling between any two qubits can be dynamically turned on and off by matching (intra-pair) and detuning (inter-pair), respectively, their frequencies, i.e., we can reconfigure the coupling structure in-situ without modifying the physical wiring of the circuit. For example, by arranging $\Delta_{k}$, $\Delta_{l}$, and $\Delta_{m}$ in Eq. (2) at three distinct frequencies, we create three qubit pairs ($Q_{2}$-$Q_{9}$, $Q_{3}$-$Q_{8}$, and $Q_{5}$-$Q_{6}$) featuring programmable intra-pair SE interactions with negligible inter-pair crosstalks, enabling parallel couplings as demonstrated in Fig. 1(b). According to the probability evolutions shown in Fig. 1(b), a characteristic gate time, $t_{\sqrt{\text{iSWAP}}}$, for each qubit pair can be be identified.

Operating multiple pairs in parallel naturally produces multiple EPR pairs [45, 46]. As the pulse sequence shows in Fig. 2(a), three EPR pairs are produced after the completion of all three SE-$\sqrt{\text{iSWAP}}$ gates [48], with the 6-qubit quantum state tomography measuring an overall state fidelity of $0.904\pm 0.018$. The inferred density matrix $\rho$ is validated by satisfying the physical constraints of Hermitian, unit trace, and positive semi-definite. We further perform partial trace on $\rho$ to obtain three 2-qubit reduced density matrices, each corresponding to a EPR pair with a fidelity above 0.93 (Fig. 2(b)).

Remarkably, our architecture allows high-efficiency generation of multiqubit GHZ states. In contrast to the previous approach where GHZ states are generated by a series of controlled-NOT (CNOT) gates [10], here all the qubits connected to the bus resonator can be entangled with a single collective qubit-resonator interaction. In the theoretical proposal [26, 49], $N$ qubits are assumed to be equally coupled to the resonator and are detuned from the resonator by the same amount $\Delta$ that is much larger than the qubit-resonator coupling. When all qubits are initialized in the same equal superpositions of $0\rangle$ and $|1\rangle$, e.g., $\left(\ket{0}-i\ket{1}\right)/\sqrt{2}$, the SE interaction does not induce any energy exchange between qubits; instead, it produces a dynamic phase that nonlinearly depends upon the collective qubit excitation number $k$ as $k(N+1-k)\theta$, where $\theta$ is determined by the effective qubit-qubit coupling strength and the interaction time. With the choice $\theta=\pi/2$, this gives rise to the GHZ state $\left(|+_{1},+_{2},...,+_{N}\rangle+i|-_{1},-_{2},...,-_{N}\rangle\right)/\sqrt{2}$, where $|\pm_{j}\rangle=\left(\ket{0_{j}}\pm i^{N}\ket{1_{j}}\right)/\sqrt{2}$ [26].

Here we apply this proposal to our experiment. We find that, though the qubit-resonator couplings are not uniform and unwanted crosstalk couplings exist in our circuit, we can optimize each qubit’s detuning and the overall interaction time to achieve GHZ states with high fidelities as guided by numerical simulation. The pulse sequence is shown in Fig. 3(a). We start with initializing the chosen $N$ qubits in $\left(\ket{0}-i\ket{1}\right)/\sqrt{2}$ by applying $\pi/2$ pulses at their respective idle frequencies, following which we bias them to nearby $\Delta/2\pi\approx-140$ MHz for an optimized duration of approximately twice $t_{\sqrt{\text{iSWAP}}}$. The phase of each qubit’s XY drive is calibrated according to the rotating frame with respect to $\Delta$, ensuring that all $N$ qubits are in the same initial state just before their SE interactions are switched on [16, 30]. After the optimized interaction time, these qubits approximately evolve to the GHZ state $|\Psi_{1}\rangle=\left(|+_{1},+_{2},...,+_{N}\rangle+e^{i\varphi}|-_{1},-_{2},...,-_{N}\rangle\right)/\sqrt{2}$, where $\varphi$ may not be equal to $\pi/2$ as in the ideal case with uniform qubit-qubit interactions; however, this phase variation does not affect entanglement. Later on we bias these $N$ qubits back to their idle frequencies; during the process a dynamical phase $\phi_{j}$ is accumulated between $|0\rangle$ and $|1\rangle$ of $Q_{j}$. Re-defining $|\pm_{j}\rangle=\left(\ket{0_{j}}\pm i^{N}e^{i\phi_{j}}\ket{1_{j}}\right)/\sqrt{2}$ ensures that the above-mentioned formulation of $|\Psi_{1}\rangle$ remains invariant, which is equivalent to a $z$-axis rotation of the $x$-$y$-$z$ reference frame, i.e., $x\rightarrow x^{\prime}$ and $y\rightarrow y^{\prime}$. Tracking the new axes is important for characterization of the produced GHZ states.

Tomography of the produced states requires individually measuring the qubits in bases formed by the eigenvectors of the Pauli operators $X$, $Y$, and $Z$, respectively. Measurement in the $Z$ basis can be directly performed. For each state preparation and measurement event, we record the 0 or 1 outcomes of each qubit and do so for $N$ qubits simultaneously; repeating the state preparation and measurement event thousands of times we count $2^{N}$ probabilities of {$P_{00...0}$, $P_{00...1}$, …., $P_{11...1}$}. Measurement in the $X$ ($Y$) basis is achieved by inserting a Pauli $Y$ ($X$) rotation on each qubit before readout. All directly measured qubit occupation probabilities are corrected for elimination of the measurement errors [43]. The $3^{N}$ tomographic operations and the $2^{N}$ probabilities for each operation allow us to reconstruct all elements of the density matrix $\rho$ (see Supplemental Material [30] for various aspects of our tomography technique including measurement stability, reliability with reduced sampling size, and pre-processing for minimizing the computational cost). The resulting 10-qubit GHZ density matrix is partially illustrated in Fig. 3(b), with a fidelity of $0.668\pm 0.025$, and the $N$-qubit GHZ fidelity as function of $N$ is plotted in Fig. 3(b) inset. The achieved fidelities are well above the threshold for genuine multipartite entanglement [27].

The full tomography technique, though general and accurate, is costly when $N$ is large. The produced GHZ states can also be characterized by a shortcut, since the ideal GHZ density matrix consists of only four non-zero elements in a suitably chosen basis. To do so, we apply to each qubit a $\pi/2$ rotation around its $y^{\prime}$ or $x^{\prime}$ axis, transforming $|\Psi_{1}\rangle$ to $|\Psi_{2}\rangle=\left(|00...0\rangle+e^{i\varphi}|11...1\rangle\right)/\sqrt{2}$ (here and below we omit the subscripts of the qubit index for clarity). The diagonal elements $\rho_{00...0}$ and $\rho_{11...1}$ can be directly measured; the off-diagonal elements $\rho_{00...0,\,11...1}$ and $\rho_{11...1,\,00...0}$ can be obtained by measuring the system parity, defined as the expectation value of the operator $P(\gamma)=\otimes_{j=1}^{N}(\cos{\gamma}Y^{\prime}_{j}+\sin{\gamma}X^{\prime}_{j})$, which is given by $\langle P(\gamma)\rangle=2\left|\rho_{00...0,\,11...1}\right|\cos(N\gamma+\varphi)$ for $|\Psi_{2}\rangle$ [9]. Polarization along the axis defined by $\cos\gamma Y^{\prime}+\sin⁡\gamma X^{\prime}$ can be measured after applying to each qubit a rotation by an angle $\gamma$ around the $z^{\prime}$ axis [16]. The oscillation patterns of the measured parity as functions of $\gamma$ confirm the existence of coherence between the states $|00...0\rangle$ and $|11...1\rangle$ (Fig. 3(c)). The fidelity of the $N$-qubit GHZ state $|\Psi_{2}\rangle$ can be estimated using the four non-zero elements, which is $0.660\pm 0.020$ for $N=10$. This value agrees with that of the GHZ state $|\Psi_{1}\rangle$ obtained by full state tomography.

A key advantage of the present protocol for generating GHZ states is its high scalability as demonstrated in Fig. 3(b). If limited by decoherence, the achieved fidelity based on the sequential-CNOT approach, $F_{\text{N,C}}$, scales approximately as $F_{\text{N,C}}\propto F_{2,C}^{N^{2}/2}$ at large $N$ (see the red dashed line in Fig. 3(b) inset), while that based on our protocol scales as $F_{\text{N}}\propto F_{2}^{N}$ (blue dashed line). Here $F_{2,C}$ ($F_{2}$) is quoted as the decoherence-limited fidelity of the CNOT gate (present protocol) involving two qubits. The falling of the experimental data (blue dots) below the scaling line when $N\geq 6$ is due to the inhomogeneity of $g_{j}$ and the crosstalk couplings. One can see that, even with the two-qubit gate fidelity above 0.99 as demonstrated in two recent experiments [10, 51], the coherence performance of the devices does not allow generation of 10-qubit GHZ state with fidelity above the genuine entanglement threshold using the sequential-CNOT approach.

In summary, our experiment demonstrates the viability of the multiqubit-resonator-bus architecture with essential functions including high-efficiency entanglement generation and parallel logic operations. We deterministically generate the 10-qubit GHZ state, the largest multiqubit entanglement ever created in solid-state systems, which is verified by quantum state tomography for the first time as well. In addition, our approach allows instant in situ rewiring of the qubits, featuring all-to-all connectivity that is critical in a recent proposal [52]. These unique features show the great potential of the demonstrated approach for scalable quantum information processing.

This work was supported by the National Basic Research Program of China (Grants No. 2014CB921201 and No. 2014CB921401), the National Natural Science Foundations of China (Grants No. 11434008, No. 11374054, No. 11574380, No. 11374344, and No. 11404386), the Fundamental Research Funds for the Central Universities of China (Grant No. 2016XZZX002-01), and the National Key Research and Development Program of China (Grant No. 2016YFA0301802). S. H. was supported by National Science Foundation (Grant No. PHY-1314861). Devices were made at the Nanofabrication Facilities at Institute of Physics in Beijing, University of Science and Technology of China in Hefei, and National Center for Nanoscience and Technology in Beijing.

1 1. General information of the experimental setup

Figure S1 provides an overview of the experimental setup, with individual components being addressed below in detail. In this section we also discuss some characteristics of the microwave and Z bias crosstalks of our device, and introduce the numerical simulation which takes into account the device characteristics.

1.1 1.1. Device

Our sample is a superconducting circuit consisting of 10 Xmon qubits interconnected by a bus resonator, fabricated with two steps of aluminum deposition as described elsewhere [31]. There is no crossover layer to connect segments of grounding planes on the circuit chip. Instead, aluminum bonding wires are manually applied as many as possible to reduce the impact of parasitic slotline modes. The bus resonator $B$ is a half-wavelength coplanar waveguide resonator with a resonant frequency fixed at $\omega_{B}/2\pi\approx 5.795$ GHz, which is measured with all coupling qubits staying in $|0\rangle$ at their respective idle frequencies (see below). The resonator has 10 side arms, each is capacitively coupled to an Xmon qubit with the coupling strength $g_{j}$ listed in Tab. S1. The Xmon qubit is a variant of the transmon qubit, each with an individual flux line for dynamically tuning its frequency and a microwave drive ($Q_{2}$, $Q_{4}$, and $Q_{6}$ share other qubits’ microwave lines in this experiment) for controllably exciting its $|0\rangle\leftrightarrow|1\rangle$ transition. The Xmon qubit reaches its maximum resonant frequency $\omega_{j}^{0}$ at the sweetpoint, where it is insensitive to flux noise and exhibits the longest phase coherence time. $\omega_{j}^{0}$ for $j=1$ to 10 in our device are around or slightly above $\omega_{B}$, whose values are roughly estimated through the flux-biased spectroscopy measurement. In this experiment all qubits are initialized to the ground state $|0\rangle$ at their respective idle frequencies $\omega_{j}/2\pi$ that spread in the range from 4.96 to 5.66 GHz, corresponding to 840 to 140 MHz below $\omega_{B}/2\pi$ (Tab. S1), where single-qubit rotations and the qubit-state measurement are performed. For initialization we simply idle all qubits for more than 200 $\mu$s [1, 32, 33, 35, 34]; For the entangling operation we dynamically bias all target qubits from their idle frequencies $\omega_{j}$ to the interaction frequency $\omega_{I}$ for a specified interaction period, following which we bias all these qubits back to their $\omega_{j}$ for measurement. Qubit coherence performance at $\omega_{I}$ can be found in Tab. S1 and in Fig. S2.

1.2 1.2. XY control

Our instrument has 7 (expendable to more) independent XY signal channels controlled by digital analog converters (DACs): 3 channels are selected to output two tones per channel and the rest 4 channels are programmed to output a single tone per channel, for a total of XY controls with 10 tones targeting 10 qubits. The 10 tones are generated with 10 sidebands mixed with a continuous microwave whose carry frequency is $\omega_{c}/2\pi=5.324$ GHz. Microwave leakage is minimized with this standard mixing method as calibrated by the room temperature electronics.

In our setup, $Q_{1}$ and $Q_{2}$ share $Q_{1}$’s on-chip XY line connecting to the 1st two-tone XY channel, $Q_{3}$ and $Q_{4}$ share $Q_{3}$’s on-chip XY line connecting to the 2nd two-tone XY channel, and $Q_{6}$ and $Q_{7}$ share $Q_{7}$’s on-chip XY line connecting to the 3rd two-tone XY channel. For the two-tone XY signals that are supposed to simultaneously act on the two qubits, we sequentially place two rotation pulses, each with a single tone targeting one of the two qubits. Each $\pi$-rotation ($\pi/2$-rotation) pulse has a length of 60 (30) ns and a full width half maximum of 30 (15) ns, designed following the derivative reduction by adiabatic gate (DRAG) control theory [50]. We perform randomized benchmark (RB) on the two qubits simultaneously, and verify that overlapping in time the two rotation pulses through the two-tone XY channel is not a problem except for $Q_{6}$ and $Q_{7}$, which might be due to the large sideband used for $Q_{7}$.

The idle frequencies $\omega_{j}$ of the 10 qubits are detuned from each other to minimize the microwave crosstalk during single-qubit rotations. For two qubits sharing the same XY line, the cross-resonance interaction reported elsewhere [36] is not a major factor due to the large detuning. For each qubit, we optimize the quadrature correction term with DRAG coefficient $\alpha$ to minimize leakage to higher levels [37], yielding the $X/2$ ($\pi/2$ rotation around the $x$ axis) gate fidelities no less than 0.998 for all 10 qubits as verified by RB on each qubit (Tab. S1 and in Fig. S2). We also select eight qubits and simultaneously perform RBs on them (see pulse sequences in Fig. S3(a)), finding that the $X/2$ gate fidelities remain reasonably high, no less than 0.993 ($Q_{7}$ is not selected due to the above-mentioned large sideband issue; $Q_{8}$’s DAC control has a smaller memory, with a maximum sequence length only half of the others’). We also perform RB simultaneously with shorter pulse sequences (smaller $m$-Number of gates) on $Q_{8}$ and $Q_{9}$, two neighboring qubits with the idle frequencies being very close, and find that the $X/2$ gate fidelities remain no less than 0.997.

We note that further optimization of the single-qubit gates are possible by shortening the gate time and introducing a slight detuning to the XY pulse to minimize the phase error [37]. We carry out the optimization procedure on $Q_{5}$ and obtain fidelities of the typical single-qubit gates, $X$, $Y$, $\pm X/2$, and $\pm Y/2$, all above 0.999 (Fig. S3(b)).

We note that although single-qubit gates are optimized at the qubits’ idle frequencies $\omega_{j}$, the gate performance may slightly degrade after a big square pulse used to tune the qubit frequency due to the finite rise-up time on the order of a few to a few tens of nanoseconds for an ideally sharp step-edge. In our experiment, after biasing the qubits from the interaction frequency $\omega_{I}$ back to their idle frequencies $\omega_{j}$, we wait for another 10 ns before applying single-qubit gates. We also use the GHZ tomography data to benchmark the gate fidelities of our $\pi/2$ rotations. For $N=7$, the density matrix of $|\Psi_{1}\rangle$ with four major elements in the $|\pm\rangle$ basis has a fidelity of $0.796\pm 0.021$. After applying a $\pi/2$ rotation to each qubit we transform $|\Psi_{1}\rangle$ to $|\Psi_{2}\rangle$ in the $|0\rangle$ and $|1\rangle$ basis, whose density matrix is measured to exhibit a fidelity of $0.771\pm 0.022$ (each 7-qubit full tomography is done within 40 minutes). Therefore we estimate that our $\pi/2$ rotation after the big square pulse for entanglement has an average fidelity around $(0.771/0.796)^{1/7}>0.995$, while a more detailed numerical simulation suggests that the average $\pi/2$ gate fidelity is $>0.993$.

1.3 1.3. Z control

Our instrument has 10 (expendable to more) independent Z signal channels controlled by DACs, which give us the full capability of simultaneously tuning all 10 qubits’ resonant frequencies. To correct for the finite rise-up time of an ideally sharp edge of a square pulse, we generate a step-edge output from the DAC and capture the waveform with a high-speed sampling oscilloscope. The measured response of the step-edge gives two time constants describing the response of the room-temperature wirings, based on which we use de-convolution to correct for desired step-edge pulses. Imperfection due to the cryogenic wirings are partially compensated using the qubit’s transition frequency response influenced by a step-edge Z bias as a caliber [38].

The effective Z bias of a qubit due to a unitary bias applied to other qubits’ Z lines is calibrated, which yields the Z-crosstalk matrix $\tilde{M}_{Z}$ as

[TABLE]

With the Z biases applied to the 10 qubits written in a column format as $\tilde{Z}_{\textrm{applied}}$ and the actual Z biases sensed by these 10 qubits written as $\tilde{Z}_{\textrm{actual}}$, we have the mapping relation of $\tilde{Z}_{\textrm{actual}}=\tilde{M}_{Z}\cdot\tilde{Z}_{\textrm{applied}}$. The Z crosstalks reach maximum at about 8% between two neighboring qubits. We note that the Z crosstalks may not contribute to the GHZ state errors as we iteratively fine-tune the Z bias of each qubit within a small range for optimal GHZ entanglement.

1.4 1.4. Qubit readout

Besides the above-mentioned 7 XY signal channels for the qubit control, our instrument has an XY signal channel that can output a readout pulse with multiple tones achieved by sideband mixing; this readout pulse is captured by a room-temperature ADC, which simultaneously demodulates the multiple (up to 20) tones and returns a pair of $I$ and $Q$ values for each tone. An impedance-transformed JPA operating at 20 mK is used before the ADC to enhance the signal-to-noise ratio. The signal-line impedance of the JPA continuously varies, in a manner of the Klopfenstein taper, so that the environmental characteristic impedance changes from 50 to 15 $\Omega$, which gives a JPA bandwidth of more than 200 MHz centered around 6.72 GHz. The fabrication procedure can be found elsewhere [31]. The JPA can be switched “ON” and “OFF” by turning on and off, respectively, an appropriate pump tone that is about twice the signal frequency. The signal transmission spectra with the JPA in the states “ON” (red line) and “OFF” (blue line) are displayed in Fig. S4, where the 10 dips correspond to the 10 readout resonators. The amplification band of the JPA, identified by the vertical difference between the red and blue lines, is tunable with a DC bias applied to the JPA.

The readout pulse is 1 $\mu$s-long, with the input tones and the power at each tone optimized for high-fidelity readout. The $j$-th tone of the readout pulse, where $j$ is up to 10 in this experiment, populates $Q_{j}$’s readout resonator $R_{j}$ with an average photon number of $n_{j}^{r}$ in 1 $\mu$s, which dispersively interacts with $Q_{j}$ with the coupling strength $g_{j}^{r}$ and shifts $Q_{j}$’s frequency downwards by an amount of $\delta\omega_{j}^{m}$. Reversely, the qubit state affects the state of its readout resonator, which is encoded in the $I$-$Q$ values at the tone $j$ of the transmitted readout pulse. At the end of 1 $\mu$s, photons in $Q_{j}$’s readout resonator leak into the circumferential transmission line at the rate of $\kappa_{j}^{r}$, and the readout resonator returns to the ground state before the next sequence cycle starts.

Repeated readout signals amplified by the JPA are demodulated at room temperature, yielding the $I$-$Q$ points at each tone on the complex plane forming two blobs to differentiate the states $|0\rangle$ and $|1\rangle$ of each qubit (see Tab. S1 and in Fig. S2). The probabilities of correctly reading out each qubit in $|0\rangle$ and $|1\rangle$ are listed in Tab. S1.

We note that the readout resonators of $Q_{3}$ and $Q_{4}$ are very close in frequency, and so are those of $Q_{6}$ and $Q_{7}$. We carefully choose the readout tones and powers to minimize the readout crosstalk if $Q_{3}$ and $Q_{4}$ are both being measured, which has a slight side-effect that the readout visibility of $Q_{3}$ drops a little bit compared with the case when $Q_{4}$ is not being measured (Fig. S4). Nevertheless, our readout choice for minimizing the crosstalk is fully verified by preparing various product states of $Q_{3}$, $Q_{4}$, and $Q_{5}$ with high-fidelity single-qubit gates and performing single- and two-qubit state tomography, with the fidelities of all reconstructed density matrices being around or above 0.99.

1.5 1.5. $XX$-type crosstalk coupling

Due to insufficient crossover bonding wires to tie the ground segments on-chip, we experience unwanted microwave crosstalk coupling between nearest-neighbor qubits. The crosstalk coupling is calibrated by measuring the qubit-qubit energy swap process around the interaction frequency $\omega_{I}$.

To understand the crosstalk coupling, we measure in detail the energy swap process of $Q_{8}$ and $Q_{9}$ as a function of the qubit detuning from the resonator $\omega_{B}$, with the result shown in Fig. S5(a) ($Q_{8}$ and $Q_{9}$ are chosen since they are the nearest-neighbor qubits with two highest sweetpoint frequencies $\omega_{8}^{0}$ and $\omega_{9}^{0}$): We excite $Q_{8}$ to $\ket{1}$ and then detune both $Q_{8}$ and $Q_{9}$ simultaneously to the same detuning $\Delta$ from the resonator, with $\Delta$ being varied; by subsequently monitoring the $|1\rangle$-state population of $Q_{8}$, $P_{1}$, as a function of the interaction time, we obtain the energy-swap dynamics of the system at various detunings. The qubit-qubit interaction strength can be inferred from the oscillation period of $P_{1}$.

We find that, in addition to the resonator mediated SE coupling $\lambda$ (see Eq. (2) of the main text) which changes sign across $\omega_{B}$, a direct $XX$-type coupling with a magnitude of $\approx 2\pi\times 2.1$ MHz, named as $\lambda_{8,9}^{c}$ for $Q_{8}$-$Q_{9}$, must be taken into account to explain our experimental data. Figure S5(b) shows the Fourier transform of the data in Fig. S5(a) along the $y$ axis, based on which the net qubit-qubit interaction as a function of $\Delta$ can be inferred, as shown in Fig. S5(d). Overall the experimental data agree well with the numerical simulation taking into account the extra $XX$-type coupling $\lambda^{c}$ (Figs. S5(c) and (d)). The $XX$-type crosstalk couplings $\lambda_{j,j^{\prime}}^{c}$ between other nearest-neighbor qubits ($Q_{j}$ and $Q_{j^{\prime}}$) are roughly estimated as the differences between the measured qubit-qubit coupling strengths and theoretical SE interaction strengths, as given in Tab. S2.

1.6 1.6. Numerical simulation

The numerical simulation in Fig. S5(c) is based on the Hamiltonian shown in Eq. (1) of the main text with the additional nearest-neighbor crosstalk coupling term $\sum_{j,j^{\prime}}{\lambda_{j,j^{\prime}}^{c}\left(\sigma_{j}^{-}\sigma_{j^{\prime}}^{+}+\sigma_{j}^{+}\sigma_{j^{\prime}}^{-}\right)}$. The decoherence impact, if considered, is included using the Lindblad master equation taking into account a Markovian environment to avoid numerical complexity. Two characteristic decay times, the energy relaxation time $T_{1,j}$ and the pure dephasing time $T_{\varphi,j}$, are used for qubit $j$. However, the non-Markovian $1/f$ character of the phase noise prevents us from directly using the $T_{2,j}^{\ast}$ values listed in Tab. S1. Meanwhile, $T_{2,j}^{\ast}$s are measured when the qubits are detuned and uncoupled, where each qubit frequency depends on the flux in its own transmon loop. On the opposite, during the resonant coupling as done in Fig. S5, the two qubits of the pair form a new system with eigenenergies that depend very weakly on each qubit flux. In other words, the eigenenergies at an anticrossing are always flat and are much less sensitive to noise, and the two hybridized levels form a kind of decoherence-free subspace [40]. Consequently, we use empirical $T_{\varphi,j}$ values ($\approx 10T_{2,j}^{\ast}$) to capture the decoherence impact if necessary, which ensures a good agreement between the numerical results and the experimental data.

During the numerical optimization of the parameters for generating GHZ states, effects of the $XX$-type nearest-neighbor couplings as listed in Tab. S2 are investigated. It is found that the introduction of the crosstalk couplings lowers the GHZ fidelities for $N>6$ but actually raises the fidelity for $N=10$ (Fig. S6).

2 2. GHZ entanglement generation

2.1 2.1. Multiqubit GHZ phase calibration at $\omega_{I}$

For the multiqubit GHZ entanglement, it is critical that the phase of each qubit’s XY drive is calibrated according to the rotating frame at $\omega_{I}$, after taking into account the extra dynamical phases accumulated during the frequency adjustment of all qubits. Here we follow the approach as done previously [16]: For simplicity we consider the product state of two qubits as $\frac{1}{\sqrt{2}}\left(|0\rangle+|1\rangle\right)\otimes\frac{1}{\sqrt{2}}\left(|0\rangle+e^{i\varphi}|1\rangle\right)$, where the extra $\varphi$ on the second qubit is seen right after the two qubits are placed on-resonance at $\omega_{I}$; we intend to find a way to adjust $\varphi$ to be zero. With the interaction Hamiltonian as $\hbar\lambda(|01\rangle\langle 10|+|10\rangle\langle 01|)$ (see Eq. (2) of the main text), the amplitudes of $|01\rangle$ and $|10\rangle$ then oscillate in time, as described (in the rotating frame) by the unitary transformation $U_{\textrm{int}}=\left[\begin{smallmatrix}1&0&0&0\\ 0&\cos(\lambda t)&-i\sin(\lambda t)&0\\ 0&-i\sin(\lambda t)&\cos(\lambda t)&0\\ 0&0&0&1\\ \end{smallmatrix}\right]$. At $t=\pi/4|\lambda|$, where $\lambda$ is negative, the two-qubit state evolves to $|00\rangle/2+(e^{i\varphi}+i)/(2\sqrt{2})|01\rangle+(ie^{i\varphi}+1)/(2\sqrt{2})|10\rangle+e^{i\phi}|11\rangle/2$, which gives equal probabilities for the four two-qubit computational states with $\varphi=0$.

Experimentally we choose $Q_{3}$ as the reference and adjust the phase of the other qubit’s microwave; we perform the check pairwisely, with the data after all phase calibrations shown in Fig. S7.

2.2 2.2. Pulse optimization

Experimentally we scan over each qubit’s detuning and the overall interaction time to optimize the GHZ state fidelity. Varying these parameters and checking the resulting GHZ density matrix is the most straightforward way, but performing the multiqubit tomography is time consuming when the qubit number becomes large. Alternatively we only repeatedly measure the reduced density matrix on a selected number of qubits, according which we optimize the pulse parameters.

For example, in order to optimize $Q_{1}$’s detuning, we perform the single-qubit tomography on $Q_{1}$ while ignoring outcomes of all other qubits, which, in the ideal case, would yield the one-qubit density matrix as $\rho_{1,\textrm{ideal}}=\left(\begin{smallmatrix}0.5&0\\ 0&0.5\end{smallmatrix}\right)$. We compare the experimental matrix with $\rho_{1,\textrm{ideal}}$, and find the best detuning at the minimum of the norm of the difference of the two matrices, as shown in Fig. S8. Other qubits’ parameters are optimized in a similar way.

3 3. 10-qubit tomography

All directly measured qubit occupation probabilities are corrected for elimination of the measurement errors before any further processing [43]. With all the corrected probability data we perform unconstrained linear inversion to obtain an initial guess of the density matrix, which might be unphysical, e.g., with small negative eigenvalues. Then we extract a Hermitian, unit trace, and positive semi-definite density matrix that is closest in distance to this initial guess [41]. We note that the state fidelity of the inferred matrix with respect to the ideal GHZ matrix typically drops by 0.03$\sim$0.04 (in absolute value) after it is validated.

3.1 3.1. Effect of reduced sample size

The $N$-qubit tomography takes $3^{N}$ tomographic operations, and for each operation $2^{N}$ occupation probabilities of the $N$-qubit computational states are measured. Once a GHZ state is generated, a tomographic operation is appended, following which the single-shot measurement yields a binary outcome for each qubit and for $N$ qubits simultaneously; running the whole pulse sequence, which includes the state preparation, the tomographic operation, and the measurement, once is one sampling event. Therefore a sufficiently large sample size, i.e., repeating the same pulse sequence many times for many synchronized binary outcomes of all $N$ qubits, is necessary to precisely count all $2^{N}$ probabilities, which would significantly slow down the measurement. For the multiqubit tomography, we maintain a fixed sample size of 3000, which is only about 3 times $2^{N}$ when $N=10$; even in this case the full tomography measurement takes about 40 hours if uninterrupted, and we have to constantly monitor our measurement to ensure that the system performance is reasonably stable (Fig. S2). We note that for $N=10$ a huge number ($3^{10}$) of operations are involved, resulting in a set of over-constrained equations to solve for the system’s density matrix, which may overcome the shortage of an insufficient sample size.

Furthermore, we carry out a test by performing the $N$-qubit tomography with a variable sample size for $N=6$ to 9. The results show that, with the sample size around $3\times 2^{N}$, the reconstructed density matrix has a fidelity very close to that with a sufficiently large sample size (Fig. S9).

3.2 3.2. Reducing the computation complexity

Here we use the unitary matrix $U^{j}$ to describe the $j$-th tomographic operation on the $N$-qubit system whose density matrix is $\rho$, where both $U^{j}$ and $\rho$ are of size $2^{N}\times 2^{N}$. The measured probability of the $k$-th computational state after the tomographic operation is therefore

[TABLE]

where $k=1$ to $2^{N}$ indexes the $N$-qubit computational states. Vectorization of the matrix $\rho$ by stacking its columns into a single column vector $\tilde{\rho}$, we have $\tilde{P}^{j}=\tilde{U}^{j}\tilde{\rho}$, where $\tilde{P}^{j}$ is the vector format of $\{P_{1}^{j},\,P_{2}^{j},\,\ldots,\,P_{N}^{j}\}$ and $\tilde{U}^{j}$ is a $2^{N}\times 4^{N}$ matrix replacing the summation terms of $U_{kl}^{j}$ and $(U_{km}^{j})^{\ast}$ in Eq. (S2). Stacking all tomographic operation matrices $\tilde{U^{j}}$ and all measured probabilities $\tilde{P}^{j}$ for $j=1$ to $3^{N}$ into $\tilde{U}$ and $\tilde{P}$, respectively, we obtain the linear equations of $\tilde{U}\tilde{\rho}=\tilde{P}$, which is used to solve for $\tilde{\rho}$ given $\tilde{P}$ and $\tilde{U}$.

$\tilde{U}$ is a column full rank matrix of size $6^{N}\times 4^{N}$. When $N$ approaches 10, it becomes extremely difficult to fully load $\tilde{U}$ into a computer’s memory and solve for $\tilde{\rho}$. Fortunately, only a small fraction of $\tilde{U}$’s elements are non-zero, so that we can use sparse matrix for storage and employ the pseudo-inverse method. With $\tilde{U}^{\dagger}$ as the Hermitian conjugate of $\tilde{U}$, we have $\tilde{U}^{\dagger}\tilde{U}\tilde{\rho}=\tilde{U}^{\dagger}\tilde{P}$, where $\tilde{U}^{\dagger}\tilde{U}$ is a symmetric and positive definite matrix of size $4^{N}\times 4^{N}$. $\tilde{U}^{\dagger}\tilde{U}$ is not only smaller in size, but also more sparse than $\tilde{U}$, which greatly reduce the complexity when solving the equations.

Here we quote the time complexity to quantitatively describe the advantage of using $\tilde{U}^{\dagger}\tilde{U}$. For a general full rank matrix of size $4^{N}\times 4^{N}$, the time complexity involved in computing the inverse operation is $O\left(\left(4^{N}\right)^{3}\right)$. For comparison, $\tilde{U}^{\dagger}\tilde{U}$ has non-zero elements only at the indices of $\left[(k-1)\times 2^{N}+k+l,\,1+l\right]$ and $\left[1+l,\,(k-1)\times 2^{N}+k+l\right]$, where $k=1,\,2,\,\ldots,\,2^{N}$ and $l$ can be any non-negative integers for the indices to be valid. The number of non-zero elements in each row of $\tilde{U}^{\dagger}\tilde{U}$ is less than $2^{N}$, and thus the number of column elementary operations needed for each row is less than $2^{N}$ during matrix inversion; the total number of operations for the $4^{N}$ rows is less than $2^{N}\times 2^{N}\times 4^{N}$. We conclude that the time complexity of solving the inverse matrix of $\tilde{U}^{\dagger}\tilde{U}$ is $O\left(\left(4^{N}\right)^{2}\right)$.

3.3 3.3. 10-qubit $\rho$ in the $|0\rangle$ and $|1\rangle$ basis

The GHZ density matrix shown in Fig. 3 of the main text has four major elements in the $|\pm\rangle$ basis, while our measurement is in the $|0\rangle$ and $|1\rangle$ basis. Here we show the partial matrix elements for the 10-qubit GHZ density matrix in the $|0\rangle$ and $|1\rangle$ basis. It is seen that all matrix elements of $\rho$ are no higher than 0.003 in amplitude (Fig. S10).