Observation of Entangled States of a Fully Controlled 20-Qubit System

We generate and characterise entangled states of a register of 20 individually controlled qubits, where each qubit is encoded into the electronic state of a trapped atomic ion. Entanglement is generated amongst the qubits during the out-of-equilibrium dynamics of an Ising-type Hamiltonian, engineered via laser fields. Since the qubit-qubit interactions decay with distance, entanglement is generated at early times predominantly between neighbouring groups of qubits. We characterise entanglement between these groups by designing and applying witnesses for genuine multipartite entanglement. Our results show that, during the dynamical evolution, all neighbouring qubit pairs, triplets, most quadruplets, and some quintuplets simultaneously develop genuine multipartite entanglement. Witnessing genuine multipartite entanglement in larger groups of qubits in our system remains an open challenge.


I. INTRODUCTION
The ability to generate quantum entanglement [1, 2] between large numbers of spatially-separated and individually-controllable quantum systems -such as qubits -is of fundamental importance to a broad range of current research endeavours, including studies of nonlocality [3], quantum computing [4], quantum simulation [5], quantum communication [6,7], and quantum metrology [8][9][10][11]. For example, in order for quantum computers and simulators to go beyond the capabilities of conventional computers, large amounts of entanglement (or other quantum correlations) must be generated between their components [12].
As such, there is an ongoing effort to generate and characterise entangled states of increasing numbers of qubits, in systems which permit preparation of arbitrary initial states, the control of interactions between constituent particles, and readout of individual sites. In such systems, the largest number of qubits entangled to date is 14, achieved in a trapped-ion system [13], followed by 10 entangled superconducting qubits [14], and 10 entangled photonic qubits [15].
Since every qubit added to an experimental system doubles the Hilbert space dimension in which the collective quantum state is described, the task of characterisation of an unknown state in the laboratory can soon become a significant challenge. Indeed, all generated entangled states of more than 6 qubits to date have been of a highly symmetric form, such as Greenberger-Horne- * These authors contributed equally to this work. † Current address: ARC Centre for Engineered Quantum Systems, School of Physics, The University of Sydney, 2006 NSW, Australia ‡ ben.lanyon@uibk.ac.at Zeilinger (GHZ) or W states, for which efficient characterisation techniques exist [16]. How to generate and detect more complex multiqubit entangled states remains an open challenge.
In this paper, we report on the deterministic generation of complex entangled states of 20 trapped-ion qubits and their partial characterisation via custom-built witnesses for genuine multipartite entanglement (GME). Our states are complex in the sense that they are generated during quench dynamics of an engineered manybody Hamiltonian and their exact description requires specifying a number of parameters that grows exponentially in the number of qubits involved. Each qubit in our system can be, and is in this work, individually manipulated and read out, as required for universal quantum computation and quantum simulation. This paper is structured as follows. Section II presents the experimental system and explains how the 20-qubit quantum states are generated and measured in the laboratory. Section III presents results of basic properties measured for the generated 20-qubit states, using established methods. Section IV introduces and applies analytically derived GME witnesses to reveal genuine tripartite entanglement in all groups of 3 neighbouring qubits. To go beyond tripartite correlations, we then turn to more computationally demanding witnesses in Sec. V, which enable GME to be detected in groups of up to 5 neighbouring qubits. Finally, we discuss our results and possible future directions in Sec. VI.

II. EXPERIMENTAL SETUP
Our register of N = 20 qubits is realized using a 1D string of 40 Ca + ions confined in a linear Paul trap, with axial (radial) centre-of-mass vibrational frequency of 220 kHz (2.712 MHz) [17]. A qubit is encoded into two long-lived states of the outer valence electron in each ion.
Under the influence of laser-induced forces that offresonantly drive all 40 transverse normal vibrational modes of the ion string, the interactions between the qubits are well described by an "XY " model in a dominant transverse field [19][20][21], with Hamiltonian Here J ij is an N × N qubit-qubit coupling matrix, σ + i (σ − i ) is the qubit raising (lowering) operator for qubit i, B is the transverse field strength (B max{|J ij |}) and σ z j ≡ Z j is the Pauli Z matrix for qubit j with eigenvectors satisfying Z |0 = − |0 and Z |1 = |1 . Interactions reduce approximately with a power law J ij ∝ 1/|i − j| α with qubit separation number |i − j|, where in this work α ≈ 1.1.
The ground state of H XY has all qubits in the state |0 . The excited states are split into m uncoupled and nondegenerate manifolds. Each manifold contains an integer number of qubit excitations (qubits in the state |1 ), with the mth manifold containing states with m qubit excitations. In previous work, we have shown that an initial state consisting of a single localised qubit excitation coherently disperses in the system, distributing quantum correlations as it propagates [22]. Here, we study the entanglement generated during the time evolution of the initial 20-qubit Néel-ordered product state |ψ(t = 0) = |1, 0, 1, ... under H XY . That is, we study the state in the laboratory that is ideally described by |ψ(t) = exp(−iH XY t) |ψ(0) . The initial Néel state |ψ(0) contains localised qubit excitations at every other site, which should ideally coherently disperse in the subsequent dynamics, entangling groups of neighbouring qubits, as we have previously shown for neighbouring pairs with a string of up to 14 qubits [19]. While Ref. [19] presented scalable tomography techniques [23][24][25], here we study multipartite entanglement dynamics in the system.
The initial state is prepared as follows. Standard Doppler cooling, optical pumping, and resolved-sideband cooling prepare the initial qubit state |0, 0, 0, ..., 0 and all 40 transverse vibrational string modes into the motional ground state [19,20]. Next, a combination of qubit-resonant laser beams that illuminate all ions simultaneously and off-resonant single-ion-focused laser beams flip every second qubit, preparing the state |ψ(0) . The interactions (laser-induced forces) simulating H XY are then turned on.
After the desired evolution time t, the interactions are turned off and the state [ideally |ψ(t) ] is measured via qubit-state-dependent resonance fluorescence, using a single-ion-resolved electron multiplying charge coupled device (EMCCD) camera. Specifically, detecting a fluo- The experimental sequence proceeds as follows. a. Standard Doppler cooling, optical pumping, and resolved-sideband cooling prepare the initial qubit state |0, 0, 0, ..., 0 and all 40 transverse vibrational string modes close (< 1 phonon per mode) to the motional ground state [19,20]. An image of the 20 ions in our trap during state-dependent fluorescence detection of the state |0, 0, 0, ..., 0 is shown. The string length is 108 µm. b. A combination of single-ion-focused and global laser beams prepares the initial state |ψ(t = 0) = |1, 0, 1, ... . c. A bichromatic light field is applied to |ψ(t = 0) , subsequently inducing qubit-qubit interactions (not shown) as described in the text [19][20][21]. d. After any desired evolution time, the interactions are turned off, leaving a nonclassical state of qubits (well approximated by) |ψ(t) . e. A combination of single-ion-focused and global laser beams is used to rotate the basis of individual qubits, determining the desired measurement basis in the next step. An example laser pattern is shown. f. The standard state-dependent resonant fluoresce technique is used to determine the state of each qubit (see Methods in Ref. [19] for more details). An example outcome, imaged on an EMCCD camera, is shown schematically. The ions are then cooled again and initialised, ready for the sequence to be repeated. rescing (nonfluorescing) ion corresponds to the measurement outcome |0 (|1 ). Such a measurement corresponds to projecting each of the 20 qubits into the eigenstates of the Pauli Z operator, for which there are 2 20 possible outcomes each corresponding to a 20-qubit projective measurement outcome. After repeated state preparations and measurements in the "Z basis" any single-qubit expectation value ( Z i ), 2-qubit correlator ( Z i Z j ), or indeed any other n-qubit "Z-type" correlator can be estimated between any qubits (up to n = 20). Performing single-qubit operations, with a single-ion-focused laser, before the aforementioned measurement process enables projective measurement of any qubit in any desired basis, and therefore the construction of any multiqubit correlation function. That is, in this work full local control over the individual qubits is available and necessary for state preparation and analysis. A conceptual schematic of our experimental protocol is presented in Fig. 1.
One approach to studying the entanglement properties of an N -qubit system is to perform full quantum state tomography to estimate the N -qubit density matrix and then develop and apply entanglement measures to that matrix. While this is technically possible for our 20-qubit system (i.e., the required measurements can each be performed in principle), it is practically not feasible as, e.g., billions of measurement bases are required. In general, the number of measurement bases required for full state tomography grows exponentially in N as 3 N . Several of us have recently shown that matrix product state (MPS) tomography can provide a pure-state estimate of states generated in quantum systems with finite-range interactions, using a number of measurements (and all other resources) that scales efficiently (polynomially) with the system size [19]. However, MPS tomography failed to produce a useful pure-state description in our present 20-qubit system, probably due to errors in preparation of the initial state that lead to mixed states and the longrange nature of the interactions present in our 20-qubit Hamiltonian.
A more favourable approach to detecting and characterising entanglement in N -qubit systems is to develop entanglement witnesses that are not a function of every element of the density matrix and can be directly measured in the laboratory with a practical number of measurements.

III. INITIAL RESULTS: MAGNETISATION AND ENTANGLEMENT
As the first experimental step, we prepare the timeevolved state of our system [ideally described by |ψ(t) ] and measure each qubit in the Z basis. The dynamical evolution of Z i (proportional to the probability of finding a qubit excitation at site i) shows how the multiple excitations in the initial state disperse, and then partially refocus at later times, in close agreement with predictions from the exact model (see Figs. 2 a and 2 b). One sees in the data and theory (Figs. 2 a and 2 b) that the qubit excitations at the ends of the string disperse more slowly than those in the centre of the string (the green color representing "delocalised" excitations appears at later times for qubits farther from the centre of the string). The qubit-qubit interactions (J ij ) are not homogeneous across the string, leading to slower evolution for qubits farther away from the centre. J ij for spin pairs near the string centre is approximately 25% larger than for spin pairs located at either end of the string. The physical origin of this effect is a combination of the Gaussian intensity profile of the laser fields involved in generating the interactions (the laser beam is centred on the middle ions and therefore weaker on the outer ions) and the intrinsic interaction inhomogeneity across the string set by the ion-string vibrational mode frequencies, for our chosen trapping parameters. Both effects are taken into account in the theoretical model presented in Fig. 2  As the final experimental step, we perform the set of measurements that would be sufficient to reconstruct (via full quantum state tomography) the density matrices of all neighbouring k = 3 qubits (qubit triplets), during the simulator dynamics. This set consists of 3 k = 27 measurement bases (with 1000 measurements performed per basis), corresponding to all possible combinations of choosing three Pauli operators. We carry out a simple scheme (choice of measured Pauli operators) that allows measurements on all 18 neighbouring qubit triplets (out of the 20-qubit string) to be performed in parallel, requiring a total of only twenty-seven 20-qubit measurement bases. From this data set, we could reconstruct the density matrices of all single qubits, neighbouring qubit pairs, and neighbouring qubit triplets. Generalising this approach to arbitrary k, all N −k+1 groups of neighbouring k-qubit density matrices in an N -qubit string can be fully characterised by measuring in 3 k bases (independent of the number of qubits N ). For fixed k, this measurement approach is clearly efficient (constant overhead) in the system size N . We nonetheless stop at k = 3, as the number of measurement bases for k = 4 is already quite demanding and, as we show, k = 3 is already sufficient to observe genuine multipartite entanglement in groups of up to five qubits.
We reconstruct the density matrices of all neighbouring qubit pairs from the experimental data, via the standard maximum likelihood method, which finds the most likely physical density matrix to have produced the data. For each of the reconstructed 2-qubit states we evaluate the genuine multipartite negativity N g , an established mea-sure for GME [26,27]. A positive value of N g for a given k-qubit state implies the existence of genuine k-partite entanglement in this state, since N g vanishes for all biseparable states. For two qubits, N g is directly related to the logarithmic negativity [28]. More details on N g are given in Sec. V. From the results one sees that all neighbouring qubit pairs become entangled during the time evolution of the system, as is shown for t = 2 ms in Fig. 2 c. Error bars, on properties calculated using the tomographically reconstructed density matrices, are derived from the finite number of measurements (1000 for each global basis) used to estimate expectation values and calculated using the standard Monte Carlo method [29].
Naturally, one may wonder if entanglement extends beyond qubit pairs, for instance, in the form of bipartite entanglement between distant qubits or in terms of genuine multipartite entanglement between groups of more than two (adjacent) qubits. In fact, one may even be tempted to ask, is multipartite entanglement not implied if every neighbouring qubit pair is entangled? The answer to this question is simply no. One can indeed have states that feature entanglement in every 2-qubit reduction, yet still feature only bipartite entanglement 1 . Nonetheless, in theory, it is often possible to detect GME purely from inspection of the reduced density matrices of overlapping groups of qubits [16]. This is the basis for our first approach to detecting GME, presented in the next section, where we derive GME witnesses purely from neighbouring two-body observables.
In general, the task of determining if and how an Nqubit quantum system is entangled is highly nontrivial. For an arbitrary known mixed state (e.g., reconstructed via full state tomography), the problem is at least NP hard and even the best known relaxations are semidefinite programs (SDPs) that are not feasible beyond 5 qubits with our computers [26,30]. However, if the density matrix is close to a given pure target state |ψ T , then targeted witnesses can be constructed to detect its entanglement without resorting to full tomography [16]. There, the canonical ansatz would be to estimate the fidelity to the ideal target state Tr(ρ |ψ T ψ T |) and check whether it is above the maximum possible fidelity of a biseparable state. That is, a corresponding witness would be W = β1 − |ψ T ψ T |, where β := min A|A ||Tr A (|ψ T ψ T |)|| ∞ , see, e.g., Ref. [16,Sec. 3.6]. While this witness could in principle be successful in detecting GME if the experimental state is indeed very close to the intended pure state, it suffers from poor noise resistance 2 and the 1 Take, for example, the k-qubit state with density operator ρ = 1 2 ( where k is assumed to be even and |φ + i,j = 1 √ 2 (|0 i |0 j + |1 i |1 j ) with the second index being understood as modulo k. The state ρ is a convex combination of k 2 -separable states, yet the reduced state ρ i,i+1 = 1 2 (|φ + φ + | + 1 4 1) of every pair of neighbouring qubits is entangled. 2 For qubits β ≤ 1 2 and thus the best possible noise resistance is task of determining the state fidelity for arbitrary pure states still requires a number of measurement settings that scales exponentially in qubit number [31]. For example, if the state is "well conditioned" [31], i.e., if only few Pauli-expectation values are of a significant size and all others vanish, one could estimate the fidelity via randomised measurements [31] with effort that scales efficiently in qubit number. Although our states are not well conditioned, in Ref. [19] we implemented this randomised measurement strategy to obtain a fidelity estimate for a 14-qubit version of the states presented here, at one time step. That experiment required preparing 5×10 5 sequential copies of the state and involved over 5 hours of data taking (with periodic recalibration of experimental parameters). Numerical simulations of the randomised measurement technique applied to the 20qubit states considered in this work show that more than 3 times the number of copies, and therefore impractical measurement time, would be required to yield accurate fidelity estimations. As such, we aim to develop novel, and more time-efficient, approaches to characterising entanglement in our 20-qubit system. As a first step in this direction, we focus on the dynamics of subsystem entanglement percolating through the system and are able to make statements about the entanglement using measurements completed in a few tens of minutes in our system.

IV. GME WITNESSES BASED ON 2-QUBIT OBSERVABLES
In this section, we construct analytical witnesses for GME based on 2-qubit observables and use them to detect GME in groups of up to three neighbouring qubits (k = 3) within the 20-qubit (N = 20) register. Recall that a (multipartite) pure state |ψ is called biseparable if there exists a bipartition A|B such that |ψ = |φ A |χ B for some |φ A and |χ B , and is called genuinely multipartite entangled otherwise. Mixed states are GME if their density operators cannot be written as convex combinations of biseparable pure states. For more details, see Appendix A.II.
Following the observation in the previous section of strong entanglement between neighbouring qubits, the first type of GME witnesses we consider is based on average fidelities of the 2-qubit density matrices with Bell states. As such, only expectation values of pairs of Pauli operators, on k-qubit subsets of choice, are required. Linear combinations of these expectation values are then evaluated and compared to their respective thresholds for 50% white noise. While this may seem strong from a bipartite intuition, many multipartite states in fact have asymptotically perfect white noise resistance, i.e., approaching 100% polynomially in system dimension, which would inevitably be missed by such simple witness constructions. biseparable k-qubit states. Surpassing a k-qubit biseparability threshold then detects genuine k-partite entanglement.
We now present a short technical summary of our method, and refer the reader to Appendix A for more details. The main quantity of interest for detecting k-qubit GME in this section is the k-qubit symmetric average Bell fidelityF Bell , which we define as (k−2)! , and the subscripts i and j denote operators acting nontrivially only on the ith and jth qubits, i.e., i is chosen unitarily equivalent to the usual triple of Pauli operators X, Y , and Z, although the unitary U i ∈ SU (2) may be chosen differently for each qubit (for each i). This ensures thatF (k) Bell can be written as a linear combination (the absolute values can be replaced by appropriate sign changes) of pairs of Pauli operators. As we show in detail in Appendix A, any quantum state of k qubits for which is genuinely k-partite entangled for any choice of U 1 , . . . , U k . For example, for k = 3 and k = 4 one can detect GME forF Bell > 0.5, is a well-known result.
If the underlying N -partite quantum state is known (the U i are known), one could directly measure all of the 3b k 2-qubit correlators Õ iÕj appearing in Eq. (2) for optimally chosen {U i }. However, when the optimal local measurements are unknown, one strategy is to measure the 6k k-qubit basis settings corresponding to the set where O (i) AB = A i i =j B j , and perform the optimisation when evaluating the corresponding results. From the 2 k outcomes of each of these simultaneous measurements of k qubits one can obtain the expectation values of all pairwise combinations of Pauli operators.
For our purposes, we exploit the fact that the results from the twenty-seven 20-qubit measurement bases already taken in the laboratory are also sufficient to calculate all the expectation values appearing in the witnesses F (k) Bell for k = 2 and k = 3. That is, they contain as a subset, all the 2-qubit observables required to calcu-lateF (2) Bell andF (3) Bell , without knowledge of the states. The results, for increasing system interaction times and for the optimisation of the U i limited to the X-Y plane for each qubit, are shown in Fig. 3. First, the witness for bipartite entanglement (F (2) Bell > 0.5) reaffirms that all qubits are (bipartite) entangled with their direct neighbours throughout the interaction time (from time 0.5 to 3.5 ms). Second, genuine tripartite entanglement between neighbouring qubit triples builds up more slowly and is initially detected at time 1.5 ms. At time 2 ms, most triples of neighbouring qubits are genuinely tripartite entangled, before the GME gradually disappears again at later times.
The experimental uncertainties (error bars) in Fig The small deviations between the theoretical and measured dynamics ofF Fig. 3 are due to experimental imperfections, which we discuss in the next section.
A number of interesting observations, based on analyses beyond those presented in Fig. 3, are now made. First, up to the evolution time presented in Fig. 3, entanglement between any 2-qubits spaced farther apart than direct neighbours was never detected-in agreement with the ideal theoretical model. For instance, on these timescales, qubit 1 does not directly become entangled with qubit 3 alone, but qubits 1, 2, and 3 do become genuinely tripartite entangled with each other. In fact, the absence of next-nearest neighbour pairwise entanglement was necessary to detect 3-qubit GME. Specifically, we found that a GME witness based only on the entanglement between direct neighbours (see Appendix A.II.2) is not able to verify genuine tripartite entanglement, and it was only possible to do so once the (separable) correlations between non-neighbouring qubits (e.g., qubits 1 and 3) are also taken into account.
Second, although there are states for which the witness of Eqs. (2) and (3) could be used to detect GME between more than 3 parties 3 , it is not sensitive enough to do so for the states presented here in our setup (neither for the theoretical predictions nor for the experimental data).
To address the question of whether genuine multipartite quantum correlations occur in groups of more than Bell > 1 12 (3 + √ 15) ≈ 0.573), respectively. The error bars for each qubit pair or triple represent 1 standard deviation of the mean in each direction. It can be seen in a that all qubits immediately become entangled with their direct neighbours and remain entangled throughout, whereas genuine tripartite entanglement is detected in time step 1.0 ms for the first time. At time step 2.0 ms, the witnessF (3) Bell indicates that all neighbouring qubit triples are genuinely tripartite entangled simultaneously (although the witness is less than 1 standard deviation above the threshold for two of these triples).
3 qubits, in our setup, we hence turn to more computationally demanding procedures, which we present in the next section.
The observed and predicted entanglement peak amplitude and dynamics for qubits near the centre and those near the ends (Fig. 3) are markedly different. We attribute those differences to the interaction inhomogeneity across the qubit string and boundary effects.

V. GME WITNESSES BASED ON NUMERICAL SEARCH
In this section, we present and apply a method that employs a numerical search to find k-qubit witnesses for GME. This search is computationally intensive: an optimisation is performed that takes computational resources that increase exponentially with k. Finding a GME witness operator for mixed 5-qubit states is already at the practical limit of our available computers and algorithms. Nonetheless, we find witnesses that succeed in detecting GME in groups of up to 5 qubits in our 20-qubit experimental system. In the following, we give a brief overview of the new witnesses and defer to Appendix B for a more detailed discussion of the technical aspects.
We make use of the genuine multipartite negativity (N g ), an established measure for GME [26,27]. A positive value of N g for a given k-qubit state implies the existence of genuine k-partite entanglement in this state, since N g vanishes for all biseparable states. The N g can be calculated given knowledge of the density matrix [26,27]. However, we have not performed a tomographically complete set of measurements for more than k = 3 qubits (we do not have the density matrices for the state in the lab, for k > 3). Our approach, to detect GME in any given group i, of k qubits, is to find a k-qubit witness operator Q (k) i whose expectation value provides a lower bound on the k-qubit N g , and which can be written as a function of the set of measurements that were carried out in our experiment. We now provide more details on this approach.
We perform a search to find a k-qubit operator Q (k) i , subject to two important constraints. First, we search for an operator which both maximises the following inequality, for a specific k-qubit state of interest ρ k i , and satisfies the inequality for all possible k-qubit states. We call Q (k) i a quantitative entanglement witness (QEW) because it provides a lower bound on N g . It is straightforward to constrain the search in this way and also to verify that any given Q (k) i is a QEW. When searching for the optimal witness, we use a theoretical model for the time-evolved k-qubit state for ρ k i . Second, we include the additional constraint that the Q (k) i can be written as a linear function of the k-qubit measurement operators (projectors) that were done, involving qubit group i. Specifically, we where the projectors P (k) s, α correspond to the marginal distributions of the twenty-seven 20-qubit projective measurement settings carried out in the lab (Sec. III), and c (k) i; s, α denote some coefficients. Here, s and α label, respectively, the qubit outcome and the local basis of the measurements, see also Appendix B.II.
The search for the optimal witnesses operator Q (k) i is carried out using a semidefinite program. The run time of the SDP is polynomial in the dimension of the Hilbert space [32,33] but the dimension of our Hilbert space naturally increases exponentially with the number of qubits k. This makes the optimisation demanding already for medium numbers of qubits: Our available computational resources are not sufficient to determine optimal witnesses for states of more than 5 qubits.
The Q (k) i which satisfies Eqs. (5) and (6), and maximises the left-hand side of Eq. (5), determines an optimal witness tailored to the target state (from a theoretical model) and the available measurements. Once this optimal Q (k) i is found we can calculate its expecta-tion value from the outcomes of the measurement done on the state in the laboratory. A witness expectation value (S (k) i ) larger than zero then detects k-qubit GME (N g > 0), for the ith group of k qubits.
The experimental results for k = 3 presented in Fig. 4 a show that all neighbouring qubit triplets soon develop GME during the dynamics, to within many standard deviations, reaching a maximum at t = 2 ms. Furthermore, for the times 2, 2.5, and 3 ms, Figs. 4 b and 4 c show that GME is detected in the majority of all neighbouring groups of 4 and 5 qubits, to within at least 1 standard deviation of experimental uncertainty. Figure 4 compares the witness results obtained from the data with those derived from two theoretical models. The first "pure" model employs the perfect pure 20-qubit time-evolved state [|ψ(t) ] and uses exact knowledge of k-qubit density matrices to optimise and apply the witnesses. The witness expectation value for the pure model yields S . Although the pure model succeeds in qualitatively describing the multipartite entanglement dynamics, the witness expectation values from the data are generally offset to lower values. A more sophisticated "mixed" model is able to explain part of this offset, which includes known imperfections in preparing the initial Néel-ordered state. Specifically, out of 1000 attempts to generate the Néel state, we observe the correct output state 829 times. In the remaining 171 cases, 146 correspond to single qubit flip errors and the rest to errors with two or more qubit flips. We attribute these errors to uncontrolled fluctuations in laser intensity and frequency, and model them as leading to the preparation of a statistical mixture of those different logical initial states, with corresponding weights. The witnesses used for the data were obtained by a search involving the mixed-model reduced states for the targets. We attribute the remaining small differences between data and theory, in Fig. 4, to additional mixing processes that occur in the laserinduced qubit-qubit interactions.
The mixed theory predictions in Fig. 4 include error bars due to the use of a finite number of numerically simulated measurements per measurement basis (1000). Error bars indicate 1 standard deviation of the mean, which is estimated as in Sec. IV, and show that the fluctuations in the data are largely consistent with those expected from such statistical noise. We conclude that, in order to witness 4-and 5-partite GME with greater statistical significance in future work, we could benefit from taking more measurements. However, it will be challenging to ensure that the experimental configuration remains stable over the longer time required to take such additional measurements.
The sizes of the error bars on both data and mixed theory points, in Fig. 4, increase with increasing k. This can be understood as follows: there are more measurement outcomes available in the data for the k = 3 witness calculations than for larger k. Amongst the twenty-seven 20-qubit measurement bases, 3-qubit measurements are repeated (duplicated) more often in the measurement pattern than 4-qubit or 5-qubit measurements, leading to better statistics.

VI. DISCUSSION AND CONCLUSION
We have experimentally generated and detected the presence of entanglement in a register of 20 qubits. In particular, we detected the dynamical evolution of genuine multipartite entanglement in the system following a quench, and developed new characterisation techniques to do so. While we cannot say that 20-qubit GME was generated, we can say that every qubit simultaneously became genuine multipartite entangled with a least two of its neighbours and, in most cases, three and four of its neighbours.
Our experimental apparatus represents the largest joint system of individually controllable subsystems to date where the presence of entanglement has been demonstrated. Each qubit can be individually controlled and qubit-qubit interactions can be turned on and off as desired (and tuned to have various forms). As such, our system has the capability to perform universal quantum simulation and quantum computation.
Confirming GME beyond groups of 5 qubits, even for the ideal states, is currently beyond our available classical computational resources and algorithms. A possible approach to overcome that problem is to exploit symmetries in the system and initial state, to reduce the size of the search space for witnesses. Another is to tune the experimental system Hamiltonian into regimes where more symmetries are apparent or approximated, e.g., infinite or nearest-neighbour-only qubit-qubit interaction ranges.
Finally, witnesses based on average Bell-state fidelities are straightforward to use and measure in the lab. As with all such witnesses, they detect entanglement without the need to carry out state tomography and can be evaluated with only a few measurements. This can be important for the detection of weakly entangled states, where estimates based on state tomography are known to overestimate entanglement [34]. Our witnesses based on brute-force numerical searching have the advantage of placing the least constraints on the form of the state in the lab and the measurements that should be taken: this can be important in the case of unknown local rotations of qubits during the dynamics. As such, our witnesses should find application beyond the present trapped-ion setting. In this section of the appendix, we introduce a method for the detection of genuine multipartite entanglement in N -qubit systems. This method is based on 2-qubit observables and does not require full state tomography. In particular, our detection criteria can be phrased as biseparability thresholds for average Bell fidelities, i.e., expectation values of linear combinations of pairs of Pauli operators. At the heart of this method lies the anticommutativity theorem (ACT) from Refs. [35,36], which we use to provide bounds on the average Bell fidelities. Although our approach is not able to detect all types of GME in multipartite systems, its advantage lies in providing linear entanglement witnesses that can be practically evaluated with only a few measurements. More specifically, our approach does not require obtaining a good estimate of the N -partite correlation tensor [37, 38] with 3 N components, but instead only needs at most 6N measurements of strings of N local Pauli operators to test for N -partite GME. As we discuss, the linearity of the witness also makes it amenable to a simple treatment of the potentially correlated statistical errors arising from deriving expectation values of bipartite observables from simultaneous measurements of N qubits.
Following the brief description of these results in Sec. IV of the main text, we now present more detailed derivations of the quantities and bounds that we consider. In Appendix A.I we briefly define and motivate the basic quantities of interest, before we construct our GME witnesses in Appendix A.II.

A.I. Framework
In this section, we explain the basic quantities and notions of interest, i.e., the anticommutativity theorem of Ref. [36] and the average Bell fidelities to establish a basis for the more detailed discussion of GME that is to follow in Appendix A.II.

A.I.1. The Anticommutativity Theorem
Let us consider a set {A n } n=1,2,...,k of self-adjoint, normalized, anticommuting operators on a Hilbert space H with dim(H) = d, i.e., Tr{A m , A n } + = Tr A m A n + A n A m = 2dδ mn (A.1) for all m, n = 1, . . . , k. The anticommutativity theorem [35,36] then states that for all states ρ ∈ L(H), A simple example for the applicability of this theorem is the set of single-qubit Pauli operators {X, Y, Z}. Since all of these operators anticommute and square to the identity, the ACT then simply requires that In other words, for single-qubit Pauli operators, the ACT is equivalent to demanding that Bloch vectors are (sub)normalized, i.e., positivity of the density operator ρ.
A less trivial example of the ACT arises for 2 qubits. Consider the set of operators where the shorthand notation for N -qubit operators is and O ∈ {X, Y, Z}. We can sort the six operators in the set displayed in Eq. (A.4) into three pairs of anticommuting operators, e.g., Since the spectra of all six operators are {±1} (with twofold degeneracy), we further have O i O i+1 2 ≤ 1, and the ACT theorem hence tells us that To see where these bounds can be of use, let us next examine fidelities with 2-qubit Bell states.

A.I.2. Average Bell State Fidelities
We now want to consider ways of quantifying how close a given 2-qubit state is to a maximally entangled Bell state. To this end, note that the density operators for the four Bell states can be written in a generalized Bloch decomposition as For any 2-qubit density operator ρ, we can then compute the fidelity with any of the Bell states. For this purpose we use the Uhlmann fidelity F, given by which reduces to F(ρ, |ψ ψ|) = ψ| ρ |ψ = Tr ρ |ψ ψ| (A.10) if one of the arguments is a pure state. For example, for the Bell state |ψ − one can use Eq. (A.8a) and the fact that all Pauli operators are traceless to find Since the only difference to the fidelities with any of the other Bell states are the relative signs between the different expectation values, we can immediately note that the fidelity of ρ with any of the four Bell states is bounded according to where we have dropped the subscript for the state ρ on the expectation values for brevity.

A.I.3. Nearest-Neighbour Average Bell Fidelity
When the system consists of more than 2 qubits, we can evaluate the fidelity with 2-qubit Bell states for any two of the constituent qubits. For simplicity, let us first consider the nearest neighbours for now and examine the case of 3 qubits. The average fidelity with arbitrary nearest-neighbour Bell states is then 1 2 F(ρ, ρ Bell,12 ) + F(ρ, ρ Bell,23 ) ≤F NN Bell , (A. 13) where we define the quantityF NN Bell as the upper bound but we refer toF NN Bell as the average nearest neighbour Bell fidelity from now on for simplicity. Next, we make use of the relation between the 1-norm || a|| 1 = n i=1 |a i | and the 2-norm || a|| 2 = n i=1 |a i | 2 1/2 in an n-dimensional vector space, i.e., the fact that where we have taken 1 = (1, 1, . . . , 1) T to be a vector whose components (w.r.t. whichever basis is chosen for a) are all equal to 1, and we have used the Cauchy-Schwarz inequality |( a, b)| ≤ || a|| 2 || b|| 2 . Combining this with the ACT theorem from Eq. (A.2) we find, e.g., Applying the same procedure to the other pairs of expectation values of anticommuting operators in Eq. (A.14), we arrive at the bound The average nearest-neighbour Bell state fidelity can of course be generalized to N qubits, i.e., the upper bound on the average nearest-neighbour Bell fidelity is However, when N is even, one triple of expectation values (w.l.o.g. for i = N − 1) remains unpaired and can only be bounded by Thus we arrive at the following upper bound on the nearest-neighbour average Bell fidelity for arbitrary Nqubit states, i.e., For, instance, for 3 qubits we have b 3 = 3 and the symmetric average Bell fidelity reads Here, we have arranged the expectation values such that it becomes immediately obvious that the triples of operators corresponding to expectation values listed directly below or above each other mutually anticommute. We can then apply the bound of Eq. (A.15) and the ACT of Eq. (A.2), e.g., as illustrated for the terms We thus arrive at the bound In fact, the same bound applies for arbitrary numbers of qubits, since all 3b N expectation values can be collected in groups of 3 mutually anticommuting operators. To see this, we use an inductive proof. Assume that we have found b N groups of three anticommuting operators for N ≥ 3 qubits and we wish to add another qubit. This means that we have to additionally consider the operators All columns contain three mutually anticommuting operators for N ≥ 3. If the original 3b N operators can be arranged in mutually anticommuting triples, then also the new set of 3b N +1 operators can be grouped in this way, which concludes the inductive step. We have already demonstrated that this statement is true for N = 3 and have hence shown that for any N ≥ 3 we have the bound Having established these general bounds that apply for arbitrary quantum states, we next examine how these bounds can be improved upon when the states in question are biseparable. This will allow us to formulate criteria for the detection of genuine multipartite entanglement.

A.II. GME Witnesses
In this section, we establish upper bounds for the nearest-neighbour and symmetric average Bell fidelities for biseparable states. These new upper bounds are below the respective bounds of Eqs. (A.22) and (A.28) and hence leave room for GME states in between. That is, any states for which the combinations of expectation values discussed above provide values beyond these biseparability bounds are GME. As we shall see, the biseparability bounds for nearest-neighbour Bell fidelities are not directly useful for detecting GME, but these bounds serve as a simple example for discussing the method of construction which will be helpful for identifying GME witnesses based on symmetric average Bell fidelities.

A.II.1. Outline of the Technique
In the following we consider bipartitions A|B of the set κ = {1, 2, . . . , N } of all N qubits, that is, we split κ into two sets, A = {a 1 , a 2 , . . . , a k |a i ∈ κ, a i = a j ∀i = j}, (A.29a) such that A ∪ B = κ and A ∩ B = ∅. For N qubits, one has 2 N −1 − 1 different bipartitions. Before we continue, let us briefly recall the definitions of biseparability and genuine multipartite entanglement. In general, a pure, N -partite state |ψ ∈ H 1, 2, . . . , N = H 1 ⊗ H 2 ⊗ . . . H N is called k-separable if it can be written as a tensor product with respect to some partition of H 1, 2, . . . , N into k ≤ N subsystems. As a special case of this definition, a pure state is called biseparable, if it can be written as a tensor product w.r.t. some bipartition, i.e., if there exists a bipartition A|B such that |ψ = |φ A |χ B . Conversely, a pure state |ψ ∈ H 1, 2, . . . , N that is not biseparable is called genuinely N -partite entangled. A mixed state with density operator ρ is considered to be genuinely multipartite entangled if it cannot be written as a convex combination of biseparable states, that is, if it cannot be written as where i p i = 1 with 0 ≤ p i ≤ 1 and |ψ (i) bisep are biseparable pure states. Note that the |ψ (i) bisep for different i can be separable w.r.t. different bipartitions. Now, consider a bipartition A|B and an operator O i O j such that i ∈ A and j ∈ B. If the system is in a pure state |ψ that is separable w.r.t. to this bipartition, i.e., if |ψ AB = |φ A |χ B , then we have When we have a triple of operators X i X j , Y i Y j , and Z i Z j for such a separable state across A|B, we have where we have used the Cauchy-Schwarz inequality in the second-to-last step and the subnormalization of the Bloch vector in the last step. The inequality (A.32) can be used to bound the Bell fidelities for pure biseparable states for different bipartitions.
As an example, consider again the nearest-neighbour average Bell fidelity for 3 qubits from Eq. (A.14). For the bipartition 1|23, we can apply (A. 32 Any pure 3-qubit state that is separable w.r.t. one or more of these bipartitions (any pure, biseparable state of 3 qubits) must hence satisfyF NN Bell ≤ 3 4 . Moreover, since any mixed state is considered to be biseparable when it can be written as a convex combination of biseparable pure states (not necessarily w.r.t. to the same bipartition), all mixed, biseparable states must also respect this bound. Conversely, the first 3 qubits of any state ρ for which are genuinely 3-partite entangled.

A.II.2. Nearest-Neighbour Average Bell Fidelity as GME Witness
In principle, the nearest-neighbour Bell fidelity could hence provide a detection criterion for GME that can be generalized to N qubits. However, at this point a remark on the detection power of this quantity is in order, since even some paradigmatic cases of genuinely tripartite entangled states for 3 qubits cannot be detected with this bound. That is, for the 3-qubit GHZ and W states |ψ (3) GHZ and |ψ (3) W (or the local unitarily equivalent 2excitation Dicke state |ψ (3) D,2 ), given by one finds nearest-neighbour average Bell fidelities of whereas the corresponding bound for detecting GME is One can hence try to improve the method or find an alternative. One way to improve the bound is by way of taking into account the purity of the biseparable states. That is, if we consider again the worst-case bipartition 1|23 for 3 qubits under the assumption that the state is separable w.r.t. this cut, i.e., that |ψ 123 = |φ 1 |χ 23 , we havē where we have used the Bloch vector b of the second qubit and the correlation tensor t = (t mn ) of qubits 2 and 3. In other words, the state |χ 23 can be written in a generalized Bloch decomposition as where σ = (σ n ) is the vector of Pauli operators (σ 1 = X, σ 2 = Y, σ 3 = Z). Now, since |χ 23 is a pure state, we have Tr(ρ 2 χ ) = 1, which translates to Since the bipartition 12|3 is equivalent and for 2|13 we have the lower valueF 2|13 NN Bell ≤ 1 2 , the "improved" biseparability bound for the nearest-neighbour average Bell fidelity for three qubits is 1 4 1 + √ 3 ≈ 0.683013. This value is still above the fidelity 2 3 obtained for pure GME states of 3 qubits in Eq. (A.44). In addition, we have also conducted a numerical search which did not reveal any pure 3-qubit states with nearest-neighbour average Bell fidelities beyond 2 3 . At the same time, one can find bisep-arable states that give values forF (3) NN Bell that are very close to 2 3 . For instance, for the state ρ bisep = |0 0| 1 ⊗ρ 23 , where |0 1 is an eigenstate of Z, and the (nearly pure) stateρ has Bloch vectors b = c = (0, 0, 0.447) T and a diagonal correlation matrix t = diag{0.894, −0.894, 1}, we findF (3) NN Bell = 0.654375. We therefore conclude that, in their present form, GME witnesses based on the nearestneighbour average Bell fidelity are practically irrelevant for 3 qubits and there is no reason to expect an improvement for more than 3 qubits.
We therefore now continue with an analysis of a different quantity, the symmetric average Bell fidelity.

A.II.3. Symmetric Average Bell Fidelity as GME Witness
In this section, we discuss the usefulness of the symmetric average Bell fidelity as a witness for GME. To this end, we again need to identify the bipartitions providing the worst (largest) upper bound forF Bell under the assumption of separability w.r.t. to the respective bipartition. Since the combination of expectation values that we consider now is symmetric under the exchange of any 2 qubits, this task is rather straightforward.
First, we consider the case of 3 qubits separately, where all three possible bipartitions (i.e., 1|23, 2|13, and 12|3) are equivalent. If the system state is pure and separable w.r.t. to any of these bipartitions, two of the triples of expectation values in Eq. (A.25) are "cut" by the bipartition and can be bounded by 1, while the remaining triple consists of three commuting observables, whose expectation values are jointly bounded by 3. For any labelling of the 3 qubits we hence havē F 1|23 Bell ≤ 1 12 3 + 1 + 1 + 3 = 2 3 . (A.43) As we discussed in Section A.II.2, this bound has to be compared with values achievable with pure GME states. For the 3-qubit GHZ-, W-, and 2-excitation Dicke states, we find symmetric average Bell fidelities which happen to coincide with the corresponding nearestneighbour average Bell fidelities of Eq. (A.44). We must hence try to improve the bound using a similar trick as before in Appendix A.II.2. Again assuming a biseparable pure state for the bipartition 1|23, we can writē Using numerical optimisation, we can also provide a pure biseparable state that comes very close to this bound. That is, for the state ρ bisep = |0 0| 1 ⊗ρ 23 , where |0 1 is an eigenstate of Z, and the pure stateρ has Bloch vectors b = c = (0, 0, 1 √ 2 ) T and a diagonal correlation Bell = 0.569036. As before, the pure state biseparability bound extends to mixed states via convexity. Thus, any 3-qubit state for which the combination of (moduli of) expectation values on the right-hand side of Eq. (A.25) exceeds 1 12 3 + √ 15 must be genuinely tripartite entangled.
Second, let us turn to the case of 4 qubits, where we are interested in bounding the quantitȳ F (4) For any pure state that is separable w.r.t. a bipartition into 1 versus 3 qubits, we find three triples of expectation values that are "cut" (each bounded by 1), while three triples pertaining to the same subsystem can be combined into mutually anticommuting triples, each bounded by √ 3, obtaininḡ F 1|234 Bell ≤ 1 24 6 + 3 + 3 As we have argued before, a biseparability bound for 3 qubits that is larger or equal to 2 3 is not very useful since even pure GME states (e.g., the 4-qubit Dicke state with two excitations) achieve only this value. We hence again turn to using the purity of the subsystems for a biseparable state |ψ 1243 = |φ 12 |χ 34 . In this case, the symmetric average Bell fidelity can be bounded bȳ Here, we encounter a different optimisation problem than before, since we no longer seek to maximize the sum of absolute values, but some quantities (e.g., | a 1 | and | a 3 |) are coupled. However, due to the symmetric form (w.r.t. the exchange of qubits 12 with 34) of the expression, we may writē We hence seek to maximize f (a, t) = 3t + 2a 2 under the constraints 2a 2 + 3t 2 = 3 and a 2 ≤ 1, which is achieved for a = 1 and t = 1 For more than 4 qubits, we can derive more general expressions using the method based on the ACT, while the exponentially increasing subsystem dimension makes bounds based on the subsystem purity unfeasible. Consider a system of N ≥ 4 qubits that is in a separable pure state w.r.t. to a bipartition into a single qubit versus the remaining N − 1 qubits. One may then identify N − 1 triples of expectation values that factorize and can be bounded by one, while the remaining b N −1 = b N −(N −1) triples form anticommuting sets of three which are each bounded by √ 3. We thus havē F 1|23...N Bell where we have made use of the fact that N −1 Note that, as required, this bound reduces to the result obtained in Eq. (A.48) for N = 4. Intuitively, it is now clear that other bipartitions will provide smaller upper bounds, since more operators are affected by the factorization. The exception being the case N = 4, where we have already seen that the separation into two sets of two provides a larger upper bound since each side then features unpaired expectation values.
To confirm this, let us briefly consider bipartitions into 2 and N − 2 qubits for N ≥ 5. In such a case, 2(N − 2) expectation values factorize for the respective pure, separable states, and one triple of operators pertaining to the two isolated qubits cannot be paired with anticommuting partners, whereas b N −2 = b N −(N −2) = b N −2(N −1)+1 triples of operators can be matched up in this way. Thus, we havē This expression provides a smaller upper bound when which is the case for N ≥ 5, as expected. We have also confirmed that this intuition holds for bipartitions into k versus N − k qubits for 3 ≤ k ≤ N − 3. We can hence formulate the biseparability bound based on the symmetric average Bell fidelity for arbitrary numbers of qubits in the following way. For any biseparable state of N qubits, the symmetric average Bell fidelity satisfies Conversely, any state that violates the inequality Eq. (A.58) is genuinely N -partite entangled. Before we discuss the practical usefulness of these witnesses, let us briefly analyze possible improvements in Appendix A.II.4.

A.II.4. Optimizing GME Witnesses Based on Bipartite Fidelities
To keep the notation simple during the derivations, we have thus far used only expectation values of pairs of the same Pauli operators, i.e., of the form | O i O j |. In practice, this corresponds to measuring the real part of certain off-diagonal elements of the density operator. To see this, consider a 2-qubit state ρ and note that Tr ρ(X 1 X 2 + Y 1 Y 2 ) = 4 Re( 01| ρ |10 ). (A.59) Of course, the off-diagonal element 01| ρ |10 need not be real for a given 2-qubit state. Here, one may note that the derivations of all bounds that we have considered so far are invariant under local unitary transformations. That is, we can replace the triple of operators for any unitary U i . This is the case because such a rotation maps a triple of anticommuting operators to another triple of anticommuting operators, and the length of the Bloch vectors also is left invariant. For instance, one could perform a rotation in the equatorial plane of the Bloch sphere, and map In the example of Eq. (A.59) this means we can pick θ 1 = 0 and θ 2 = − π 2 to obtain In particular, there exist rotation angles θ 1 and θ 2 such that Tr ρ(X 1X2 +Ỹ 1Ỹ2 ) = 4| 01| ρ |10 |. (A.63) In general, one hence has the freedom of N independent transformations U i ∈ U (2) to optimize the GME witnesses presented so far. In an experimental setting, this optimisation can be done a priori if the quantum state ρ that one expects to produce (approximately) in the experiment is known. However, if the underlying state is unknown, one may also measure all combinations of 2-qubit Pauli operators for all pairs of qubits within the set of N qubits (amounting to 9b N = 9 2 N (N − 1) 2-qubit measurements) and perform the optimisation on the experimental data. For instance, the results forF In addition to a posteriori optimisation, one may perform some of these 2-qubit measurements on different pairs simultaneously if the individual outcomes for each qubit are recorded. For instance, in a register of N ≥ 3 qubits one may obtain the expectation values of X 1 X 2 and X 2 Y 3 from measuring the first and second qubit in the eigenbasis of X and the third in the eigenbasis of Y and recording all three outcomes in each run. For mea-surements of this kind, such as have been performed in our experiment, data used to estimate different 2-qubit expectation values may be correlated, which has to be taken into account in the calculation of the estimate for the variance of the GME witness. This is explained in more detail in Ref. [ where the numerical values for 3 and 4 qubits are 3−1 ≈ 0.0915064, respectively. The gap between the bounds is hence largest for N = 3, and shrinks with increasing number of qubits. It is hence expected that there is some finite N for which no GME states exist that are detected by our witnesses, and at this point, we cannot say for which N this occurs.
For 3 qubits, we have already found examples of genuinely tripartite entangled states that can be detected, i.e., the 3-qubit W state |ψ (3) W and the 2-excitation Dicke state |ψ (3) D,2 which provide symmetric average Bell fidelities of 2 3 . The experimental results discussed in the main text (see Fig. 3) further show thatF (3) Bell is also a useful witness for mixed states produced in realistic situations. Beyond 3 qubits, the witnessesF (N ≥4) Bell (optimized only over rotations in the X-Y plane, see Appendix A.II. 4) have not been able to detect GME in our experimental setting. However, we know that 4-qubit states exist, e.g., the 4-qubit 2-excitation Dicke state |ψ (4) D,2 which could be detected in this way, sinceF (4) Bell (|ψ (4) D,2 ) = 2 3 . Unfortunately, the 2-excitation Dicke state for 5 qubits only provides a value ofF (5) Bell (|ψ (5) D,2 ) = 0.6 whereas F (N )bisep Bell = 1 20 (7 + 3 √ 3) ≈ 0.61. Tentative searches for other 5-qubit states for whichF (5) Bell exceeds the biseparability bound have been unsuccessful thus far. The question of whether GME states exist that can be detected with our method for N ≥ 5 hence remains open. We now discuss the genuine negativity (GMN) of Ref. [27] that we use to quantify GME in the experiment. We present its definition as a convex-roof construction and the alternative way of writing it in terms of a semidefinite program, which turns it into a numerically computable measure of entanglement for an arbitrary mixed state.
We start by introducing the notation and presenting preliminary definitions. Note that a bipartition of {1, . . . , N } can be specified by a subset A ⊂ {1, . . . , N } and its complementĀ = {1, . . . , N }\A. With this, we define the partial transposition on A for an operator X A ⊗ XĀ, where X A and XĀ act on the Hilbert spaces associated to A andĀ, respectively, as ( The definition then extends to any operator on the N -particle Hilbert space by linearity. Next, the negativity of a bipartite quantum state ρ with respect to the bipartition A|Ā is given by the sum of the negative eigenvalues of the partially transposed density matrix, i.e., where λ i (X) denotes the ith eigenvalue of an operator X. With this, the GMN of an N -particle state ρ is given by where the inner minimisation is over the possible bipartitions A|Ā of {1, . . . , N } and the outer minimisation is over decompositions ρ = i p i ρ i , where {p i } i is a probability distribution and ρ i are density matrices. For pure states, this reduces to N g (|ψ ψ|) = min A|Ā N A|Ā (|ψ ψ|).
For mixed states the GMN can still efficiently be computed using numerical tools from the field of semidefinite programming. This follows from the fact that the GMN is given by the optimal value of the following optimisation where Q, K A , and R A are operators acting on the Hilbert space. The GMN is zero for all (bi)separable states and, therefore, a nonzero value provides a way to certify GME. More precisely, the GMN is nonzero for any state that cannot be written as a PPT mixture. Recall that a multipartite state ρ is called a PPT mixture if it admits a mixed state decomposition ρ = A p A ρ A where {p A } A is a probability distribution and ρ A has a positive partial transposition (we say, ρ A is PPT) with respect to the bipartition A|Ā. That is, formally, ρ T A A ≥ 0. As noted earlier, a state is called biseparable if it can be written as a convex combination of states that are separable with respect to one bipartition A|Ā. Since any separable state is PPT, every biseparable state can be written as a PPT mixture. Consequently, a state with nonzero GMN is GME. Moreover, the GMN quantifies the entanglement in the sense that a state ρ is more entangled than σ if N g (ρ) ≥ N g (σ). The underlying mathematical property is that the GMN is nonincreasing under so-called full LOCC operations. In particular, from this property it follows that no GME state can be generated from a non-GME state with local operations only.
A further beneficial property of writing the GMN as in Eq. (B.2) is that this yields an entanglement witness, that is, an observable that provides ideally a sharp lower bound to the GMN as in Eq. (5) in the main text and may be accessed experimentally. For our purposes note that whether a witness can be measured depends on the measurements that are available in our experiment. We therefore discuss those measurements next. The procedure of obtaining the entanglement witnesses that are accessible for us is then described subsequently in Appendix B.III in more detail.

B.II. Accessible Measurements on k Neighbouring Sites
Here, we discuss the operators that can be measured locally on k consecutive sites with the data available in our experiment. Our starting point is hence the set of all possible observables whose measurement outcomes can be obtained from the 27 measurement settings that we mentioned in Sec. III. As we noted there, these measurement settings can be used to obtain estimates of the expectation values of all possible products of the identity and the three Pauli operators on all groups of three neighbouring qubits in the chain of 20 qubits. On such triples of neighbouring qubits this thus allows us to estimate the expectation value of any operator, in particular, any possible entanglement witness. However, the situation is different for more than three neighbouring sites. Regarding this case, recall that in the 27 settings we consider, the 20 qubits are measured simultaneously such that the accessible information turns out to be more than just knowing the three-body reductions of neighbouring qubits. Also as a consequence of how we choose these 27 settings, each of them has the property that the local bases in which sites i and i + 3 are measured are identical. To illustrate what this implies for the accessible operators, let us consider the example of four sites. In terms of projective qubit measurements, a basis of operators that we can measure with the 27 settings is given by B 4 = {P (4) s, α } s, α where P (4) s, α denotes the projector onto the state | s, α = |s 1 , s 2 , s 3 , s 4 α1,α2,α3,α1 , where s ≡ (s 1 , s 2 , s 3 , s 4 ) ∈ {↑, ↓} ×4 and α ≡ (α 1 , α 2 , α 3 ) ∈ {1, 2, 3} ×3 . In words, these operators comprise all projective qubit measurements where the first and the fourth qubit are measured in the same direction. Notably, with the Pauli operators X, Y and Z, this set of operators spans the same subspace as the operators {σ α1 ⊗ σ α2 ⊗ σ α3 ⊗ σ α1 } α , where α ∈ {0, 1, 2, 3} ×3 and σ 0 = 1, σ 1 = X, σ 2 = Y , σ 3 = Z. For the construction of the witness we will, however, use the projectors as our basis set. Accordingly, generalising this to arbitrary number of sites, we denote the projectors that form a basis of the operators that can be measured on k consecutive sites with the 27 settings by B k . As a further remark, let us mention that the number of elements in B k grows exponentially with k as |B k | = 27 × 2 k .

B.III. Accessible Quantitative Entanglement Witnesses
Next, we turn to the construction of an entanglement witness that is fully accessible from the available information provided by the 27 measurement settings. The main step is to solve the optimisation given in Eq. (B.6) below, which provides a quantitative witness. The witness can then be evaluated using the experimental data.
We distinguish between a simple entanglement witness, which is an operator Q that has a positive expectation value Tr(Qρ sep ) ≥ 0 for any separable state ρ sep and for which there exists at least one entangled state ρ with Tr(Qρ) ≤ 0, and a quantitative entanglement witness, i.e., an operator that fulfills the property of being an entanglement witness and additionally provides a lower bound to the GMN via (minus) its expectation value S = −Tr(Qρ) of the form S ≤ N g (ρ) (B.4) for any ρ, as noted in the main text. Now, considering k neighbouring sites in the chain of 20 qubits, the available information from the 27 settings is determined by the projectors in the set B k as described above. Then, an entanglement witnessŴ is accessible from this information if it can be written as a linear combination of operators from B k , i.e., if i; s, α ∈ R, since, in this case, it is fully determined by the set B k .
Here, for a given state ρ i; s, α in order to find a quantitative witness that provides the best lower bound to the GMN of ρ. As the computation of the GMN itself, see Eq. (B.2), this optimisation can be expressed as a semidefinite program. That is, with the definition In the results we present, we choose = 5 × 10 −3 and use three iteration steps. Except for the tolerance , this step does not deteriorate the bound. We then further add two simple constraints to Eqs. (B.6) and (B.7). First, we choose the coefficients as c (k),n i; s, α ∈ [−1, 1]. Second, we add the semidefinite condition K A ≤ 1 to Eq. (B.6). Note that this condition has previously been employed in an earlier version of the GMN, see Refs. [26,27]. However, there the quantitative witnesses are exact functions of the operators K A and R A , whereas here they are only related via an inequality. Solving Eqs. (B.6) and (B.7) with the constraints introduced in this section results in the bounds S (k) i presented in the main text. We observe that the additional modifications significantly help to improve the lower bounds.