Generation of genuine all-way entanglement in defect-nuclear spin systems through dynamical decoupling sequences

Multipartite entangled states are an essential resource for sensing, quantum error correction, and cryptography. Color centers in solids are one of the leading platforms for quantum networking due to the availability of a nuclear spin memory that can be entangled with the optically active electronic spin through dynamical decoupling sequences. Creating electron-nuclear entangled states in these systems is a difficult task as the always-on hyperfine interactions prohibit complete isolation of the target dynamics from the unwanted spin bath. While this emergent cross-talk can be alleviated by prolonging the entanglement generation, the gate durations quickly exceed coherence times. Here we show how to prepare high-quality GHZ$_M$-like states with minimal cross-talk. We introduce the $M$-tangling power of an evolution operator, which allows us to verify genuine all-way correlations. Using experimentally measured hyperfine parameters of an NV center spin in diamond coupled to carbon-13 lattice spins, we show how to use sequential or single-shot entangling operations to prepare GHZ$_M$-like states of up to $M=10$ qubits within time constraints that saturate bounds on $M$-way correlations. We study the entanglement of mixed electron-nuclear states and develop a non-unitary $M$-tangling power which additionally captures correlations arising from all unwanted nuclear spins. We further derive a non-unitary $M$-tangling power which incorporates the impact of electronic dephasing errors on the $M$-way correlations. Finally, we inspect the performance of our protocols in the presence of experimentally reported pulse errors, finding that XY decoupling sequences can lead to high-fidelity GHZ state preparation.


Introduction
Generating and distributing entanglement is one of the most fundamental yet non-trivial requirements for building large-scale quantum networks.Entanglement is the key ingredient of quantum teleportation, robust quantum communication protocols, or alternative quantum computation models such as the one-way quantum computer [1,2,3], or fusion-based quantum computation [4].In the context of cryptography, entangled states ensure secure communications via quantum key distribution or quantum secret sharing protocols [5,6,7,8,9,10,11].Error detection [12,13] and correction [14,15,16,17,18,19] schemes are based on the encoding of information onto entangled states, such that errors can be detected or corrected by utilizing correlations present in entangled states.At the same time, distributed entanglement increases the sensitivity and precision of measurements [20,21].
Defect platforms offer an electronic qubit featuring a spin-photon interface that enables the generation of distant entanglement [22,23,24], as well as nuclear spins that serve as the longlived memories needed for information storage and buffering.Several protocols based on remote entanglement or entanglement within the electron-nuclear spin system have already been demonstrated in these platforms, including quantum teleportation [25], quantum error correction [17,18,19], enhanced sensing [26,27,28], and entanglement distillation [29].The electronnuclear spin entanglement is usually realized through dynamical decoupling (DD) sequences; large entangled states can be created by applying consecutive sequences with the appropriate interpulse spacings.By tuning the interpulse spacings, one can select different nuclear spins to participate in entangling gates while decoupling the electron from the remaining always-on coupled nuclear spin bath [30].
A major issue is that entangled states often decohere faster than product states [31,32], so entanglement must be preserved for long enough times to complete quantum information tasks.Additionally, in defect platforms, the unwanted spin bath sets a lower bound on the duration of entangling gates, which should be long enough to suppress unwanted interactions that lower the quality of target entangled states.Despite this challenge, generation of electron-nuclear entangled states has been realized experimentally in NV [33], in SiV centers in diamond [34,35], and in SiC defects [36].In particular, Bradley et al. demonstrated experimentally in an NV defect in diamond electron-nuclear GHZ states involving up to 7 qubits [33] by combining DD sequences with direct RF driving of the nuclear spins, accompanied by refocusing pulses to extend the entanglement lifetime.This direct driving was necessary to improve the selective coupling of individual nuclear spins.Although successful, this method introduces experimental overhead and heating of the sample due to the direct RF driving of the nuclei, potentially creating scalability issues.Another problem is that as the number of parties contributing to the entangled state increases, a larger gate count and longer sequences are needed, exacerbating dephasing.Failure to provide optimal isolation from unwanted nuclei further deteriorates the electron-nuclear entangled states.Therefore, a more efficient approach to generating electron-nuclear spin entanglement within time constraints and verifying its existence is necessary for large-scale applications.
In this paper, we address these challenges by introducing a framework for preparing high-quality GHZ states of up to 10 qubits within time constraints.Using experimental parameters from a 27 nuclear spin register well-characterized by Taminiau et al. [37], we show how to improve sequential entanglement generation methods by minimizing the cross-talk induced by the unwanted nuclear spin bath.We find that it is possible to prepare GHZ-like states with maximal all-way correlations in excess of 95% and gate errors lower than 0.05% due to residual entanglement with unwanted nuclei.We show how to prepare GHZ-like states with single-shot operations, reducing the gate times at least twofold compared to the sequential scheme while offering decoupling capabilities comparable to the latter.We present a closed-form expression for the M -tangling power of the evolution operator and use this to develop a method for verifying genuine multipartite entanglement.Remarkably, this metric depends only on two-qubit Makhlin invariants and is closely related to the one-tangles we introduced in Ref. [38].This simplification allows us to systematically determine the DD sequences that maximize all-way correlations as desired for generating multipartite entangled states.Further, we analyze the entanglement of mixed electron-nuclear states and derive a non-unitary M -tangling power that captures residual entanglement links arising from unwanted nuclei.We incorporate electronic dephasing errors into the M -tangling power and derive a simple closedform expression.Finally, we study the impact of pulse control errors on the M -way correlations.
The paper is organized as follows.In Sec. 2, we discuss methods for generating electron-nuclear spin entangled states using DD sequences.In Sec. 3, we introduce the M -tangling power of an evolution operator.In Sec. 4, we show how to prepare GHZ M -like states through sequential or single-shot entanglement protocols.In Sec. 5, we quantify the entanglement of electron-nuclear mixed states by tracing out unwanted spins.In Sec. 6, we study the non-unitary M -tangling power of the electron-nuclear target subspace, which encodes correlations arising from the entire electron-nuclear system.Finally, in Sec.7 we study the effect of electronic dephasing errors as well as pulse control errors on the M -tangling power. 2 Establishing electron-nuclear entanglement

Multi
One way to generate electron-nuclear entanglement is through DD sequences, which are trains of π-pulses applied on the electron spin that are interleaved by free evolution periods.These sequences are constructed by concatenating a basic π-pulse unit of time t a certain number of iterations, N .Well-known examples are the Carr-Purcell-Meiboom-Gill (CPMG) [39,40,41,42] and Uhrig (UDD) [43,44] sequences, which have been used experimentally for example in Refs.[45,30,17,33,46], or proposed for defects in Ref. [47].In the defect-nuclear spin system, DD sequences serve a two-fold purpose; i) they average out the interactions between the electron and unwanted nuclei and, ii) under so-called resonance conditions, they selectively entangle a target nuclear spin with the electron [30].These resonances correspond to specific values of t, and they are generally distinct for each nuclear spin, as they depend on the hyperfine (HF) parameters of each nucleus.By composing consecutive sequences with different parameters (i.e., unit times t and iterations N ), we can thus create entanglement between the electron and multiple nuclei.We refer to this standard experimental approach of entangling one nuclear spin at a time with the electron as "sequential".In Ref. [38], we showed that alternatively single-shot operations can be used to entangle the electron with a subset of nuclei from the register, significantly reducing gate times compared to the sequential approach.Figure 1 summarizes the two approaches to generating electron-nuclear entanglement.
Figure 1(a) depicts the multi-spin scheme implemented by a single-shot entangling operation, and Fig. 1(b) shows the sequential scheme which requires M − 1 entangling gates to prepare GHZ M -like states.The M − 1 entangling gates are realized by composing sequences of different interpulse spacings and iterations of the unit, such that a different nuclear spin is selected from the register based on the resonance condition.
GHZ M states are created by initializing the electron in the |+⟩ state, polarizing the nuclear spins in the |0⟩ state [48,17], and then performing electron-nuclear entangling gates via DD sequences.Probabilistic (deterministic) initialization of the nuclear spins in |0⟩ can be achieved through measurement-based (SWAP-based) initialization [49].After state initialization, one can perform a DD sequence to generate entanglement.Setting the sequence unit time t to be a nuclear spin resonance time forces that nuclear spin to rotate along an axis conditioned on the electron spin state.If the axes are anti-parallel and along the ±x direction, and the unit is iterated an appropriate number of times, the nuclear spin can accumulate a rotation angle of ϕ = π/2, realizing a CR x (±π/2) electron-nuclear entangling gate [30].Repeating this process by changing the unit time t and iterations N of the subsequent sequence (and assuming the k-th target spin is ro-tated only during the k-th evolution; trivial nonentangling rotations could happen on the k-th nucleus during the remaining evolutions) creates a state equivalent to GHZ M = (|0⟩ ⊗M +|1⟩ ⊗M )/ √ 2 up to local operations: (1) where m k = ±1, m (k) j can be +1 or −1 for the kth nucleus, and we defined |±y⟩ = (|0⟩±i|1⟩)/ √ 2. For example, if all the gates are CR x (π/2), then the resulting state is (|0⟩|−y⟩ ⊗k + |1⟩|y⟩ ⊗k )/ √ 2. This is the idea behind the sequential entanglement protocol.The multi-spin entanglement protocol produces a more complicated state since the electron-nuclear evolution now becomes a CR xz gate (see Ref. [38]), but the resulting state is still locally equivalent to GHZ M .In this latter setup, multiple nuclei are simultaneously entangled with the electron via a single-shot operation.
Quantifying the quality of GHZ M -like states using a fidelity overlap is computationally costly since we would need to optimize over local gates.This optimization becomes more difficult for larger entangled states.In the following section, we introduce an entanglement metric that allows us to quantify genuine all-way correlations in an arbitrarily large electron-nuclear spin system, using only the information of the evolution operator generated by π-sequences.With this analysis we are able to translate the M -tangling capability of such an evolution operator irrespective of the total system size, into simple spin-spin correlation metrics between the central spin and each of the M −1 nuclear spins.Our analysis throughout the paper is completely general and is not restricted by the choice of DD sequence or by the choice of the defect system.It is further generically valid for any type of nuclear spins present in the register (provided they have spin I = 1/2).The code used to simulate all our following results can be found in [50].

M -tangling power
Genuine multipartite entanglement of GHZ Mlike states can be verified with metrics such as entanglement witnesses [51,52,53], concentratable entanglement [54], or the so-called Mtangles [55,56,57].We choose to work with the latter, which are usually defined in terms of a state vector (or, more generally, a density operator).The M -tangles, τ M (|ψ⟩) ∈ [0, 1], distinguish the GHZ entanglement class from other entanglement classes.For example, the three-tangle saturates to 1 for the GHZ 3 state, whereas it vanishes for the |W⟩ = 1/ √ 3(|100⟩ + |010⟩ + |001⟩) state.The M -tangles are invariant under permutations of the qubits and SLOCC, and are entanglement monotones1 [58].Since they are defined based on a state vector, their calculation requires that we make an assumption about the initial state of the register, which is then evolved under the π-pulse sequence.We circumvent this issue by focusing instead on the capability of a gate to saturate allway correlations and thus prepare states locally equivalent to GHZ M .
We introduce the M -tangling power of a unitary, which we define to be the average of the Mtangle ϵ p,M (U ) := ⟨τ M (U ρ 0 U † )⟩, over all initial product states, ρ 0 .The formulas of M -tangles expressed in terms of the state vector can be found in Appendix A.1 and in Refs.[55,57,56].
In our system, the total evolution operator (of the electron and M − 1 nuclear spins) produced by DD sequences has the form: where σ jj ≡ |j⟩⟨j| are projectors onto two of the levels in the electron spin multiplet, and is the rotation that acts on the l-th nuclear spin and, in general, depends on the electron spin state.Due to the fact that the evolution is controlled only on the electron, we find that we can compactly express the M -tangle of a given state as [see Appendix A.2]: where ρ ⊗2 is the density matrix of two copies of the M -qubit system living in the Hilbert space H M 2 .P is the product of projectors onto the antisymmetric (or symmetric) space of the i-th subspace and its copy, i.e.
where P (±) j,j+M = 1/2(1 ± SWAP j,j+M ), and with SWAP j,j+M = α,β∈{0,1} |α⟩ j |β⟩ j+M ⟨β| j ⟨α| j+M .Note that we have fixed system "1" to be the electron's subspace.(If the control qubit corresponds to a different subspace, then for odd M , the product of antisymmetric projectors needs to exclude the control system from the copies, and instead we apply P (+) on the sectors of the control qubit).For even M ≥ 4, Eq. ( 3), holds for density matrices ρ obtained by any arbitrary U , whereas for odd M , Eq. ( 3) holds strictly for states obtained by CR-evolutions.(We verified this numerically, comparing with the exact expressions of τ M that hold for states obtained by arbitrary evolutions.)By averaging Eq. ( 3) over all product initial states, we prove that the M -tangling power of an operator U has the simple expression [Appendix A.2]: with , where d = 2 (dimension of qubit subspace).Equation ( 5) is exact for any arbitrary U (M ), for M even and M ≥ 4, whereas for odd M it holds only for CRtype evolution operators of the form of Eq. ( 2).The simple structure of CR-type evolution operators generated by π-sequences allows us to derive analytically a closed-form expression for the M -tangling power for an arbitrarily large nuclear spin register coupled to the electron.In Appendix A.3, we prove that the M -tangling power of the operator U of Eq. ( 2), is given by: where G 1 is the so-called Makhlin invariant [59], which we have shown that for a π-pulse sequence is given by [38]: 1 ) is the rotation angle that the kth nuclear spin undergoes when the electron is in the |0⟩ (|1⟩) state, and n 1 characterizes the evolution of the kth nuclear spin under the DD sequence and is related to the bipartition entangling power of the nuclear spin with the remaining register, i.e., the nuclear one-tangle given by ϵ nuclear p = 2/9(1 − G 1 ) [38].
The nuclear one-tangle is a faithful metric of selectivity; its minimization ensures that the spin is decoupled from the register, while its maximization implies that the particular spin attains maximal correlations with the register.
The remarkable simplicity of Eq. ( 6) allows us to check whether a controlled unitary of the form of Eq. ( 2) is capable of preparing GHZ Mlike states.Saturating the bounds of all-way correlations between the electronic defect and the nuclear register translates into minimizing G 1 (or, equivalently maximizing ϵ nuclear p ) of the nuclei that will be part of the GHZ state, and implies that the maximum M -tangling power is At this point, we should mention a caveat of the M -tangle metric.In Ref. [56] it was mentioned that the four-tangle cannot discriminate the entanglement of a GHZ 4 -like state from that of bi-separable maximally entangled states, e.g., |Φ + ⟩ ⊗2 , where |Φ + ⟩ is one of the four Bell states; in both cases the four-tangle is 1.
Analogously, the even Mtangle cannot discriminate the entanglement of a GHZ M state from the entanglement of M/2 copies of two-qubit maximally entangled states |Bell⟩ ⊗M/2 (see also Ref. [60]).Nevertheless, in our system we have the extra condition that G 1 of all M − 1 nuclei should be minimized to saturate all-way correlations.Hence, if we fail to minimize at least one out of the M − 1 {G 1 } quantities, we know that the evolution operator will never be able to prepare genuine multipartite entanglement between all M parties.Another way to understand this is to consider the example of 4 qubits, and the entangling gate CR ).If we act with CR x (π/2) ⊗3 on arbtirary initial states, we can at most prepare a GHZ 4 -like state, but we will never be able to prepare two individual Bell pairs, since correlations are constrained to be distributed among all parties connected by the electron (since we are dealing with a central spin type system).Similar statements hold for the multispin method that utilizes CR xz gates.Therefore, we can safely use the M -tangle metrics to detect genuine multipartite electron-nuclear (and nuclear-nuclear) spin entanglement.
Let us further comment that the M -tangling power is derived by averaging over all product (and pure) states, so it only tells us how well the dynamics of U can saturate the all-way correlations of pure product states.No conclusion can be made, however, if one has a statistical mixture of such states (i.e., a mixed state).
In the following section, we show how to optimize the sequential or multi-spin schemes and guide the selection of nuclei from the register to generate genuine all-way correlations within time constraints.

Makhlin-invariant of j-th spin
Table 1: Meaning of the symbols that we use in the flowchart of Fig. 2 to describe the optimization steps for the sequential and multi-spin protocols.

Sequential scheme for GHZ 3 states
We first begin with the task of generating GHZ Mlike states using the sequential scheme, which requires M − 1 consecutive gates.Hereafter, we refer to the collection of M − 1 entangling gates as the composite evolution.
Under this scheme, to entangle the nuclei with the electron, we need to ensure that during each evolution we rotate conditionally only one particular nuclear spin, meaning we maximize its one-tangle.At the same time, all other nuclear one-tangles should be minimal in each evolution, so that we suppress cross-talk arising from the nuclear spin bath.The procedure of identifying optimal nuclear spin candidates and sequence parameters is shown in Fig. 2 and explained in more detail in Appendix G. Depending on the size of the GHZ state, we set different tolerances for unwanted/target one-tangles, and gate time restrictions within T * 2 of the nuclei (see Table 2 of Appendix G).At the end of this procedure, we are left with different options for selectively entangling a single nuclear spin with the electron.
After identifying the optimal sequence parameters and nuclear spin candidates, we compose the M − 1 entangling gates and create a "case" of a composite evolution operator, from which we can extract the Makhlin invariants G 1 associated with the evolution of each nuclear spin.This process allows us to calculate the M -tangling power of the composite gate (using Eq. ( 6)) and verify if it can prepare genuine multipartite entanglement.Further, as we showed in Ref. [38], we can analytically calculate the induced gate error due to residual entanglement links with the unwanted nuclei.This unresolved entanglement could make the target gate deviate from the ideal evolution, which we take to be the composite M − 1 entangling operations in the absence of unwanted spins; minimization of this gate error thus ensures a better quality of the GHZ-like states that we produce.
We start with the simplest case of generating GHZ 3 -like states using the sequential scheme.Throughout the rest of the paper, we focus on the CPMG sequence (t/4−π −t/2−π −t/4) N , where t is the unit time, π represents a π-pulse, and N is the number of iterations of the basic unit.For concreteness, we consider an NV center in diamond and define the qubit states of the electron to be |0⟩ ≡ |m s = 0⟩ and |1⟩ ≡ |m s = −1⟩.We further use the hyperfine parameters of the 13 C nuclear spins from the 27 nuclear spin register characterized by the Taminiau group [37], and following their conventions we label the nuclear spins as Cj with j ∈ [1,27].Despite these choices, the analysis that follows is general and valid for any π-sequence or electronic defect in diamond or SiC, and for nuclei with spin I = 1/2 (e.g., 13 C in diamond, 13 C and 29 Si in silicon carbide).
In Fig. 3(a), we show the nuclear one-tangles (scaled by the maximum value of 2/9) for 39 different composite evolutions (labeled as "cases") that can prepare GHZ 3 -like states.Each realization corresponds to a unique nuclear spin combination that is entangled with the electron (the

Reject case
Multi-spin scheme Sequential scheme Figure 2: Flowcharts that summarize the steps we follow to find optimal sequences that prepare GHZ M -like states, utilizing the sequential or multi-spin scheme.The various symbols are defined in Table 1.text above the bars indicates which 13 C nuclear spins are involved).The first (second) bar of each case corresponds to the one-tangle of the first (second) target nuclear spin at the end of the composite evolution of two entangling gates.Each entangling gate is implemened using an optimal unit time t of the sequence and an optimal number N of iterations of the unit.Since we optimize t close to the resonance time of the respective nuclear spin, the dot product of the nuclear axes is n 0 • n 1 ≈ −1.Thus, the maximization of the first (second) one-tangle happens when the first (second) nuclear spin rotates by ϕ = π/2 [see Eq. ( 7)], during the first (second) evolution.(Note that for each CPMG evolution it holds that ϕ 0 = ϕ 1 ≡ ϕ, see also Ref. [38]).
The composition of the two entangling operations gives rise to total rotation angles ϕ 0 ̸ = ϕ 1 for each nuclear spin.The nuclear rotation angles ϕ := ϕ 0 are shown in Fig. 3(b), and the dot products of the rotation axes are depicted in Fig. 3(c).Due to the composite evolution, the total nuclear rotation angles deviate from π/2, and the rotation axes are no longer anti-parallel for all cases.However, as we noted in Ref. [38], G 1 can be minimized as long as n 0 • n 1 ≤ 0. Due to the unwanted one-tangle tolerances we impose, we ensure that only one out of the M − 1 gates will conditionally rotate the target spin we intend to entangle with the electron.The remaining evolutions will approximately leave the entanglement of that specific spin with the electron intact.For example, the first (second) target spin will evolve approximately trivially (i.e., in an unconditional manner) during the second (first) evolution, hence preserving the high entanglement between all parties at the end of all the gates.This behavior is confirmed in Fig. 3(d) where we show the M -tangling power (scaled by [d/(d + 1)] M , which is higher than 0.99 across all different cases.In Fig. 3(e), we further show the total gate time, which we have restricted to less than 2 ms so that we are within the coherence times of the nuclei (T * 2 ∈ [3, 17] ms; see Ref. [61]).Further, in Fig. 3(f), we show the gate error due to residual entanglement of the target subspace with unwanted nuclei.We note that for various realizations, the gate fidelity can exceed 95% (e.g., cases "2", "7", "19", "23","30", "37"), implying small levels of cross-talk with the remaining 25 unwanted nuclear spins.

Multi-spin scheme for GHZ 3 states
Next, we proceed with the multi-spin scheme, which can prepare GHZ 3 -like states with a singleshot entangling operation.To identify optimal nuclear spin candidates and sequence parameters, we follow a slightly different procedure, shown in Fig. 2 and explained more in Appendix G.In this case, we require that two nuclear one-tangles are maximized for the same sequence parameters and under the gate time restriction of 2 ms.At the same time, we require small values of unwanted one-tangles to minimize the cross-talk with the remaining nuclear spin bath.We summarize the results of this method in Fig. 4. In Fig. 4(a), we show the one-tangles of two nuclei that are maximized by a single-shot gate for 15 different DD sequences.To accept these cases as "optimal", we require that the target one-tangles are greater than a threshold, here set to 0.9. Figure 4(b) shows the nuclear rotation angles, and Fig. 4(c) shows the dot products of the nuclear rotation axes for each case.Since now the sequence unit time will not in principle coincide with the resonance of both target spins (since their HF parameters are distinct), the nuclear rotation axes will not be anti-parallel for both nuclei.However, when n 0 • n 1 ≤ 0 for a particular nuclear spin, then if it is also rotated by the appropriate angle, the one-tangle can be maximal (see for example cases "2", "7", "8", "11" and Ref. [38] for a more detailed explanation).Thus, the feature of nonantiparallel axes can be compensated by rotating the nuclear spins by angles that deviate from π/2.On the other hand, if n 0 • n 1 > 0, G 1 cannot go to 0, and the one-tangle can never be 1 (see case "4", for nuclear spin C5).In Fig. 4(d) we show the M -tangling power (scaled by [d/(d + 1)] M ) of the multi-spin gate.As expected, the Mtangling power saturates to 1 when the nuclear one-tangles are also maximal.In Fig. 4(e), we show the gate time of the multi-spin operation, and in Figure 4(f), the gate error.We note that for various cases, the multi-spin gate duration is at least two times shorter (see for example, cases "4", "11", "22" ) compared to the sequential one.Thus, the multi-spin gates can be equally reliable as the sequential gates, but with the additional advantage of being significantly faster compared to the latter.

Generating GHZ M states
We can use similar ideas to extend the size of the GHZ M -like states.In Fig. 5(a), we show the Mtangling power of composite evolutions that can prepare GHZ M -like states up to M = 10 qubits.The nuclear one-tangles for all different realizations are provided in Appendix H. Figure 5(b) shows the gate error of the composite evolution due to the presence of the remaining 27 − (M − 1) nuclear spins, whereas Fig. 5(c) shows the total gate time of the composite sequences.As the size of the GHZ-like states increases, it becomes more difficult to eliminate the cross-talk between the target nuclei since we need to ensure that only one out of the M − 1 evolutions rotates conditionally one of the M − 1 target nuclei.However, we still find acceptable cases that can prepare up to GHZ 10 -like states with M -tangling power over 0.95.The infidelity tends to increase as we perform more entangling gates because it becomes more likely that at least one of the M − 1 entangling operators will induce cross-talk between the target subspace and the remaining nuclear spin bath.However, it is still possible to create GHZ 6 -like states within ∼ 2 ms, and even GHZ 9or GHZ 10 -like states within 4 ms, a significant improvement compared to gate times reported in Ref. [33], which could be as high as 2.371 ms for controlling only a single nuclear spin.
We also study the possibility of creating GHZ M -like states for M > 3 via the multi-spin scheme.In principle, due to the more constrained optimization that requires the simultaneous maximization of M −1 nuclear one-tangles, we expect that we will find fewer acceptable cases that respect the bounds we set for target/unwanted onetangles.In Fig. 6(a), we show the nuclear onetangles for various cases of single-shot entangling operations, which can create up to GHZ 9 -like states.Figure 6(b) shows the M -tangling power of the single-shot gate, whereas Fig. 6(c) shows the gate error.We note that the infidelity can in principle be less than 0.1, but the M -tangling power is not high enough in all cases, due to imperfect entanglement between target nuclei and the electron.
Figure 6(d) shows the gate time of the multispin operation.Once again, the gates are significantly faster compared to the sequential scheme.While for the sequential protocol increasing the number of parties of the GHZ M -like states implies durations of the total gate that only increase, the multi-spin operations follow a different trend.To understand this, note that the sequential scheme requires higher-order resonances (i.e., longer unit times) to suppress unintended couplings with unwanted nuclei.Higher-order resonances are also needed to ensure that only one target nuclear spin is rotated conditionally by each of the M − 1 sequential gates.Thus, increasing the GHZ size makes the requirement to resort to higher-order sequences, in principle, more stringent.On the other hand, increasing the size of the GHZ M -like state in the multi-spin protocol means that we need to gradually reduce how well we decouple the electron from the spin bath and allow interactions of the electron with increasingly more nuclear spins.In contrast to the sequential scheme, we thus need to shorten the unit time, which gives us acceptable cases of realizing GHZ M -like states.
Of course, the duration of the single-shot gate also depends on the number of times, N , we repeat the unit and on the constraints we impose on the target all-way correlations and unwanted nuclear one-tangles.In principle, a shorter basic unit can lead to faster operations.For example, it is remarkable that multi-partite entangled states can be prepared within 1 ms or even faster with only one entangling gate (see for example, the case of the GHZ 8 -like state).Experimentally, depending on the nuclear HF parameters, this method could be highly beneficial for preparing entangled states involving nuclear spin clusters fast, whereas the entanglement could be boosted using distillation protocols [29].
Overall, both the sequential and multi-spin protocols can generate high-quality GHZ M -like states with reasonable gate times.Our formalism based on the M -tangling power allows us to identify optimal scenarios of preparing entangled states with minimal cross-talk.Interestingly, a hybrid entanglement generation scheme involving both single-shot and sequential gates could offer a more realistic path to scalability by drastically reducing the number of entangling operations and gate times.

Entanglement of mixed states
So far we have made no assumption about the initial state of the system and focused instead on the capability of the gates to produce entan-gled states.We have used the gate error due to residual entanglement links and the unwanted nuclear one-tangles as a metric of mixedness in the GHZ M -like states prepared by DD sequences.We now work with a density matrix for the total system and inspect the entanglement of a target subspace after we trace out the unwanted nuclear spins.This partial trace operation will result, in general, in a mixed state for the target subspace, for which we cannot use directly the M -tangles that are defined for pure states.
In the case of mixed states, the M -tangles are calculated via so-called convex roof constructions [55,62]: where the minimum is taken over all ensembles with probabilities p i and pure states |ψ i ⟩ that reconstruct the mixed state density operator.The ensemble that gives the minimum value of the tangle is known as optimal.Such convex roof extensions are necessary to ensure an entanglement monotone, but since the minimum of a sum of convex functions is not always convex, we need to take the convex hull of Eq. ( 8) (see also Refs.[63,62]).
Finding the optimal ensemble is a non-trivial task as it entails searching over all possible decompositions of the mixed state.This task becomes even harder when the rank of the reduced density matrix is large.However, for our problem, we find that starting from an arbitrary pure state of the total system and tracing out any number of nuclei, the maximum rank of the reduced density matrix is 2.This is due to the form of the total evolution operator, which can produce at most GHZ M -like states (if extra single-qubit gates are not allowed) given the appropriate initial state.Further, the maximum rank proves that the creation of |W ⟩-like states is impossible if one allows only entangling gates generated by DD sequences.In Appendix I we derive the eigenvectors and eigenvalues of the density matrix, of which two are nonzero, provided the electron starts from a superposition of |0⟩ and |1⟩.
The analytical expression of the eigenvalues and eigenvectors of the reduced density matrix of the target subspace allows us to find the optimal ensemble that minimizes the M -tangle of the mixed state.This can be done using the methods developed in Ref. [62].The first step is to di-agonalize the reduced density matrix of the system.Based on the eigendecomposition we can then construct the trial state: where |v ± ⟩ are the eigenvectors of the rank-2 mixed state and λ ± the two nonzero eigenvalues.
Since any ensemble can be obtained by acting with unitaries on the diagonalized reduced density matrix, minimizing the entanglement of the trial state by varying the angle χ implies that we obtain the entanglement of the optimal ensemble that reconstructs the mixed state.The angle χ controls the relative phase of the eigenvectors.
Our system consists of 28 qubits, so simulating the total density matrix is computationally hard.The fact that we obtain the reduced density matrix analytically [see Appendix I] allows us to find the impact of the unwanted spins on the target subspace.Working with a density matrix necessitates that we specify the system's initial state.In the literature, it is often assumed that the nuclear spin bath starts from the maximally mixed state.This assumption, however, would result in maximal entanglement in the target subspace (provided we perform close to perfect gates and once the electron-nuclear target register is purified) free from cross-talk, since the maximally mixed state of the unwanted spins would remain invariant under any unitary evolution.This is far from true since the target subspace suffers from cross-talk from unwanted nuclei introduced by the entangling gates.On the other hand, considering a pure state for the entire electron-nuclear register is again unrealistic, but we choose this convention to find how classical correlations due to the partial trace operation manifest in the target subspace.
For the following analysis, we will use the 39 cases we found for creating GHZ 3 -like states (generalization to larger dimension of the target subspace is straightforward) from Sec. 4 via the sequential scheme.We further use the analytical expression of the three-tangle, first introduced in Ref. [55].To prepare an entangled state with genuine all-way correlations, we need to find the appropriate initial state for the target subspace, upon which we act with the entangling operations.In most cases, this state needs to be |+⟩|0⟩ ⊗M −1 .(Recall that the CR x (π/2) gates are locally equivalent to CNOT.)If this initial state does not give rise to a three-tangle of at least 0.95, we then optimize over the initial threequbit state (we assume a sampling of 0.05π for θ ∈ [0, π] and 0.1π for γ ∈ [0, 2π] for each qubit state |ξ⟩ = cos(θ/2)|0⟩ + e iγ sin(θ/2)|1⟩).It is thus possible that other initial states might give slightly larger three-tangles than the ones we will show later on.We further assume that the initial state of the bath is |0⟩ ⊗27−(M −1) .
We begin with the following setup.In Sec. 4, we found 39 cases of preparing GHZ 3 -like states using the sequential scheme, via the composition of two consecutive DD sequences.The composite evolution rotates both target and unwanted nuclei.Some unwanted spins might evolve slightly conditionally on the electron.Thus, we need to update the initial state of all spins under the composite evolution.Using the analytical expression of the reduced density matrix, we obtain its eigenvalues and eigenvectors and arrange the 39 cases in increasing value of the largest eigenvalue, p = λ + .In Fig. 7(a), we show the three-tangle of the pure state prepared using the composite evolutions of Fig. 3 by ignoring the presence of unwanted nuclei.Using the three-tangle, we verify that all-way correlations are maximal across all cases, as expected based on the fact that we approximately saturated the M -tangling power (and we also use the appropriate three-qubit initial state).We further show the three-tangle of the two eigenvectors, |v ± ⟩, as a function of the largest eigenvalue p, obtained by diagonalizing the mixed state (after tracing out unwanted nuclei).The three-tangle of |v + ⟩ coincides with the three-tangle of |v − ⟩, which means that the unwanted bath creates a mixture of two terms, which both have the same M -way entanglement and are GHZ 3 -like states.Another feature we observe is that the all-way correlations of the eigenvectors |v ± ⟩ coincide with those of the pure GHZ 3 -like state of the target electron-nuclear register.Hence, tracing out the unwanted spins makes the reduced state classically correlated, but the amount of entanglement in each term of the mixed state is unaffected.
In Fig. 7(b), we then use the two eigenvectors to build the trial state and find the minimum of the three-tangle as we vary χ for each value of p. Since we are tracing out 25 unwanted nuclei, we see that the entanglement of the mixed state reduces substantially.One simple way to understand this feature is to think of a trial (pure) twoqubit state which is a superposition of two Bell states.When we have the equal superposition, the two-qubit concurrence goes to 0, whereas for different weights between the two terms, the concurrence is 0 < C(|ψ⟩) ≤ 1.A similar thing happens here, with the difference that the superposition terms are two GHZ 3 -like states.In Appendix J, we show that the three-tangle of the superposition of two orthogonal GHZ 3 states is much more sensitive compared to the concurrence of the superposition of two orthogonal Bell states for the two-qubit case.Even for large values of p (i.e.p = 0.9) the minimum three-tangle is around ∼ 0.64 for the GHZ 3 mixture compared to the minimum concurrence which is ∼ 0.8 for the Bell mixture.In Fig. 7(b), we further show the convex hull of the minimum value of the three-tangle.In the scenario we are studying, the minimum value of the three-tangle is not convex since the pure state of the target subspace does not have perfect correlations due to slightly imperfect gates (i.e., gates that cause over-or under-rotation of the nuclei).If the gates create perfect all-way entanglement then both eigenvectors are perfect GHZ 3 -like states, and the minimum of the threetangle is convex in terms of p, and is given by We should further comment that our results depend highly on the initial state of the system.Since the evolution is conditional on the electron, if the initial state of the unwanted spins is not the appropriate one to generate the maximum possible entanglement, cross-talk will be suppressed.Also, as the system's initial state is generically mixed, the rank of the reduced density matrix after tracing out unwanted nuclei will be larger than 2.However, the simple model we assumed in this section captures qualitatively how the Mway entanglement changes as the target subspace becomes more mixed, as a result of residual entanglement with unwanted nuclei.
In the next section, we show another way of understanding the unwanted residual entanglement in terms of the non-unitary entangling power of the entire quantum channel, where now the partial trace operation is translated into the operator-sum representation [64].

Non-unitary entangling power
An alternative way to include the impact of unwanted nuclei on the target subspace is to use the Kraus operators associated with the partial trace channel.In Ref. [38], we derived the Kraus operators for an arbitrary number of nuclei coupled to a single electron spin.Here, we use this result to derive the M -tangling power of the evolution that the target system undergoes due to the presence of the unwanted nuclei.Because this is a generalization for non-unitary dynamics, we will refer to this metric for brevity as the non-unitary Mtangling power.
The density matrix of the total system evolves under the unitary U given by Eq. ( 2).If we trace out the unwanted nuclei, the target subspace evolves under the quantum channel E(ρ) = , where E k are the Kraus operators corresponding to unwanted nuclei, L is the number of total nuclei in the bath, and K is the number of target nuclear spins.In this case, the M -tangling power (where M = K + 1) of the quantum channel reads [65]: (10) with the same definitions for Ω p0 and P we introduced in Sec. 3. Since the total evolution is controlled on the electron, it is possible to derive a closed-form expression for the M -tangling power of the channel, which reads (see Appendix A.4): (11) with W (j) given by: x,1 n x,0 n We find that the M -tangling power of the channel is upper-bounded by the unitary M -tangling power of the target space, i.e.
The equality is satisfied when the unwanted subspace evolves trivially (i.e., irrespective of the electron's state), in which case it holds that G (j) 1 = 1 and W (j) = 0, ∀j.Therefore, whenever ϵ p,M (E) < ϵ p,M (U ), we know that the unwanted spin bath possesses nonzero correlations with the target subspace.This is an alternative way to see that the entanglement of the mixed state is lower compared to the pure case, a feature we already observed in Sec. 5.The remarkable simplicity of Eq. ( 11) allows us to obtain full information of the entanglement distribution within the electron-nuclear register, irrespective of the total number of qubits in the system, or number of nuclear spins we are tracing out.Using ϵ p,M (E), we can find how well the gate we design saturates the M -way correlations of the target subspace, given the fact that there might be residual unwanted correlations linking the target subspace with unwanted spins; those unwanted correlations are encoded in the Kraus operators.An extra advantage of this approach is that we never need to assume an initial state of the electron-nuclear register, as we did for example in Sec. 5. To compare the non-unitary M -tangling power against ϵ p,M (U ), we consider as an example the optimal cases we identified for preparing GHZ 4 -like states from Sec. 4 for either the sequential or the multi-spin schemes.The results for the sequential scheme are shown in Fig. 8(a).The dark blue bars show the unitary M -tangling power across 29 realizations, and the pink bars show the non-unitary Mtangling power.As expected, ϵ p,M (E) is lower in all cases due to unresolved crosstalk, and this deviation is enhanced when the gate error is larger [see Fig. 5(b)].A similar feature is observed in Fig. 8(b) for the 7 cases of preparing GHZ 4like states using the multi-spin scheme.Specifically, for case "5", for which we found the lowest gate error emerging from residual entanglement across the 7 cases [see Fig. 6(b)], the non-unitary M -tangling power is closer to the unitary one.Hence, this verifies our previous observation that the mixed electron-nuclear states we create in the target subspace remain GHZ M -like, but their entanglement is reduced whenever correlations between the target subspace and the unwanted spin bath are present.
In Fig. 8, we also display with green bars an approximate expression for ϵ p,M (E), where we set W (j) = 0, ∀j.We note that the approximate expression agrees very well with the exact expression of ϵ p,M (E).
Finally, let us discuss the computational complexity related to the various entanglement metrics we introduced.The calculation of the G 1 Makhlin invariants and hence of target/unwanted one-tangles scales linearly with the number of nuclear spins, since each G 1 quantity describes the evolution of a single nuclear spin.Thus, the M -tangling power, ϵ p,M (U ), can be calculated without computational difficulty.Similarly, ϵ p,M (E) scales linearly with the size of the total size of the register (including target and unwanted nuclei), and can also be calculated without difficulty.The infidelity of the target gate due to the presence of unwanted spins requires a summation over 2 L−K Kraus operators.Therefore, the starting point to identify optimal cases for preparing GHZ M -like states is to calculate the target/unwanted onetangles in order to maximize (minimize) target (unwanted) correlations.To further inspect the optimal cases in terms of gate error, one can then proceed with this metric to obtain full information about the dynamics of the system.

Additional error analysis
In this section, we study the effect of pulse control errors and dephasing errors that can occur on the electronic qubit.To find the impact of those errors, we consider only their contribution to the M -tangling power of the target subspace.

Target subspace M -way entanglement under electronic dephasing
We begin by considering dephasing errors for the electronic qubit.The dephasing error can be represented via the Kraus operators [66]: In Appendix B, we prove that the "dephased" M -tangling power of the target subspace (consisting of the electron and target nuclei) has a simple closed-form expression.In particular, for a CPMG decoupling sequence of unit time, t, and of N iterations, it takes the form: ), (15) where R is given by: In the case where we compose CPMG sequences of different unit times t j and iterations N j (as we do, for example, in the sequential scheme), the function R takes the form: We verify that in the limit of λ 0 = 1 and λ 1 = 0, we recover the ideal ϵ p,M (U ) for the target subspace.We further find that the M -tangling power in the presence of electronic dephasing errors is bounded from below by 50%.We begin by studying the impact of the electronic dephasing for the 39 cases of preparing GHZ 3 -like states via the sequential scheme.The results are shown in Fig. 9(a), where we consider λ 0 = 0.98, λ 1 = 1 − λ 0 , and dephasing angles ranging from θ −1 = 400 µs to θ −1 = 50 µs.We see that for θ −1 = 400 µs, we preserve high Mtangling power over 90% for most realizations.For a dephasing angle θ −1 = 100 µs, we find several cases with M -tangling power around 80-85%.For a dephasing angle θ −1 = 50 µs, the M -way correlations can drop to less than 60% across many cases.
In Fig. 9(b), we repeat the calculations for the 27 cases of preparing GHZ 3 -like states using the multi-spin protocol.The M -tangling power of the multi-spin protocol remains higher in the presence of electronic dephasing compared to the sequential scheme.We even observe that for the smallest dephasing angle of θ −1 = 50 µs, we have several cases which exceed 80% M -way correlations.This behavior is expected since the multispin scheme is, in principle, faster than the sequential protocol, and hence it is less sensitive to electronic dephasing errors.
Figure 9(c) shows the M -tangling power of case "1" of the sequential scheme, as a function of λ 0 and θ −1 , whereas Fig. 9(d) shows the M -tangling power of case "1" of the multi-spin scheme.In the sequential scheme, we control the nuclei C1 and C18 (with a gate time ∼ 1843.2 µs), whereas, in the multi-spin protocol, we control the nuclei C1 and C2 (with a gate time ∼ 1366.5 µs).Although the results of Fig. 9(c) and Fig. 9(d) should not be compared directly due to the different total timings, we note that the multi-spin scheme is more robust to dephasing errors.In Appendix C, we provide another example where we compare the  performance of the multi-spin scheme with the sequential, but now for case "16".Although case "7" of the multi-spin scheme has a total sequence time longer than the sequence time of the sequential scheme, we find again that the multi-spin scheme is more robust to dephasing.We attribute this feature to the more complicated dephasing channel of the sequential scheme.The expression of ϵ p,M (E deph ) for the sequential scheme is described by different unit times t j and iterations N j since we are composing CPMG sequences to control each nucleus.We believe that this feature can potentially lead to M -way correlations that are more sensitive to electronic dephasing errors.

Target subspace M -way entanglement under pulse errors
We now consider the effect of pulse errors during the control of the electronic spin.We assume that the pulse error results in an over-or under-rotation along the x-axis of the electron.We model such errors by modifying the perfect R x (π) acting on the electron into R x (π + ϵ) = e −iπ/2(1+ϵ)σx .Such rotation angle errors yield an evolution operator which is no longer blockdiagonal.Consequently, we cannot use the Mtangling power we found in Eq. ( 6).Additionally, we cannot study the M -tangling power for an odd-sized system since, as we mentioned in Sec. 3, the expression of Eq. ( 5) is only applicable to CR-type evolution operators.This difficulty in describing the M -way correlations on the evolution operator level for odd M arises from the fact that averaging over all initial states is not straightforward due to the expression of the odd M -tangle we start with.Nevertheless, by modeling the rotation angle error in this way, we still preserve unitary dynamics, and we can find the impact of such errors on the M -tangling power of even-sized systems numerically by making use of Eq. ( 5).First, we consider systematic errors and study their impact on the cases of preparing GHZ 4like states via the sequential or multi-spin protocols.In Fig. 10(a) and Fig. 10(b), we show the M -tangling power for the 29 cases of preparing GHZ 4 -states, assuming a systematic overrotation of ϵ = 2% and ϵ = 8% respectively.The dark blue bars depict the error-free case.The light blue bars correspond to the XY2 sequence, whereas the orange bars to the CPMG sequence, including the over-rotation errors.In the case of 2% error, we see that CPMG can provide high values of M -way correlations only for a few cases (e.g., case "5", case "11", case "17", case "24").However, if we consider the XY2 sequence whose one unit is t/4 − (π) X − t/2 − (π) Y − t/4, we find that the impact of the error on the quality of the GHZ states is negligible.For an over-rotation error of 8%, the M -way correlations accumulated via the CPMG sequence are lower than 50%, whereas the XY2 sequence still provides high entanglement across almost all cases (all cases besides case "9" and case "20").
We observe a similar behavior of systematic over-rotation errors for the multi-spin scheme.We depict the results for the multi-spin protocol in Fig. 10(c) and Fig. 10(d), where we assume a systematic over-rotation error of 2% and 8%, respectively.For 2% error, we find that ϵ p,M (U ) obtained via the XY2 scheme is extremely robust, although the pulse iterations per case are in general high: N = 128 for case "1", N = 76 for case "2", N = 124 for case "3", N = 136 for case "4", N = 78 for case "5", N = 223 for case "6", and N = 243 for case "7".For an 8% error shown in Fig. 10(d), case "6" and case "7" display the most prominent reduction in M -way correlations across all cases, which is consistent with the large number of decoupling unit iterations.
The robustness of XY-sequences to pulse errors has already been reported extensively in the literature [42,67,68,69].In Ref. [69], a quantum process tomography protocol called boostrap was used to determine control pulse errors for the π and π/2 pulses.It was mentioned therein that the π/2 pulses show twice as much variation than πpulse errors indicating that the pulse edges have a larger impact on shorter pulses.The pulse errors were taken as constant during each run and the same for different runs.The reported error values for π-pulses were ϵ x = ϵ y = −0.02(rotation errors for π X and π Y pulses), q x = 0.005 (x-axis component of π Y rotation), q y = 0 (ŷ-axis component of π X rotation) and q z = ±0.05(ẑ-axis component of π X and π Y rotations).More importantly Ramsey, spin-echo, XY4, CPMG, and UDD simulations with those parameters were in strong agreement with experimental data.The estimated pulse error was 1%.In Ref. [68] XYsequences were used to cancel out the systematic pulse errors, suppress decoherence of the electron,  as well as suppress artificially injected magnetic noise, which was introduced in the same stripline that was used for the control pulses.The fidelity of a single π-pulse on the electron was again estimated to be 99% using calibration techniques from Ref. [41,67].Therefore, since the π-pulse microwave control can typically exceed 99% for defect platforms, we expect that the M -tangling power we evaluate for the 2% systematic rotation error will be close to the experimentally observed one.In Appendix E, we perform another simulation with the aforementioned error parameters and verify that the M -tangling power we obtain for the ideal CPMG sequence is close to the one obtained via the XY2 sequence in the presence of pulse imperfections.This reveals that combining our protocols with the XY2 decoupling sequence is sufficient to prepare high-fidelity GHZ M states in the presence of realistic experimental errors.
We continue the pulse error analysis, assuming random over-/under-rotation pulse errors, which we sample from a normal distribution with standard deviation σ = 0.01.In Fig. 11(a), we show the mean M -tangling power for the XY2 (light blue bars) and for the CPMG (orange bars), using the 29 cases of preparing GHZ 4 -like states with the sequential scheme.In each case, we run 500 independent trials where we sample anew from the normal distribution (per iteration of the decoupling unit), and we evaluate the mean ϵ p,M (U ) across the 500 trials.For comparison, we also display the ideal ϵ p,M (U ) per case.CPMG and XY2 perform on par when random pulse errors are present.The error bars mark the range of ϵ p,M (U ) for the 500 trials per case.We observe that for the sequential scheme, the M -tangling power in the presence of random pulse errors is higher than 0.7 across all cases.
In Fig. 11(b), we show the mean M -tangling power for the 7 cases of preparing GHZ 4 -like states using the multi-spin scheme.We sample errors from the normal distribution with standard deviation σ = 0.01 and repeat the calculation 500 times per case to collect statistics and evaluate the mean M -tangling power.Once again, we observe that the XY2 and CPMG sequences perform on par in the presence of random over-/under-rotation errors.The error bars per case correspond to the range of ϵ p,M (U ) for a total of 500 random trials.
Overall, we find that the generation of high-quality GHZ M states is possible with the current state-of-the-art experimental control and constraints.Our formalism provides a detailed analysis of the entanglement generation in the presence of errors and allows us to address the problem of designing optimal entangling gates that saturate all-way correlations.

Conclusions
Genuine multipartite entangled states are an essential component of quantum networks and quantum computing.Exploiting the full potential of nuclear spins in defect platforms for largescale applications requires precise and fast entanglement generation.We showed under what conditions decoupling sequences produce gates capable of maximizing all-way correlations and quantified the entanglement capability of the gates through the M -tangling power.Using this formalism, we guided the selection of sequence parameters and nuclear spin candidates to prepare high-quality GHZ M -like states by appropriate driving of the electron spin.We improved the sequential entanglement generation scheme, pushing the gate time to as low as 4 ms for preparing GHZ-like states of up to 10 spins.We also studied the possibility of direct entanglement generation, which performs on par with the sequential scheme in decoupling capabilities, with the extra advantage that the latter approach drastically reduces the gate count and speeds up the entanglement generation.Further, we studied the entanglement of mixed states, revealing that the M -way correlations of a target subspace are sensitive to residual entanglement with the unwanted nuclei.We introduced a non-unitary M -way entanglement metric, which additionally captures correlations between the target subspace and the unwanted nuclei and showed that it is upper-bounded by the unitary M -tangling power.We derived a nonunitary M -tangling power for the target subspace in the presence of electronic dephasing errors, revealing that the multi-spin protocol can have in principle superior performance compared to conventional schemes due to its shorter implementation time.Finally, we studied the impact of pulse errors on the target subspace M -tangling power and showed that our protocols combined with XY decoupling sequences can provide highfidelity preparation of GHZ states in the presence of realistic experimental errors.Our results pave the way for the systematic and efficient creation of multi-partite entanglement in spin defect systems for quantum information applications.

A Mathematical treatment of M -way entanglement A.1 Formulas of M -tangles expressed in terms of the state vector
For any even M ≥ 4 the M -tangle is given by [55,56]: where in the last line, we have linearized the equation of the M -tangle by introducing a second copy and projecting onto the antisymmetric space of subsystems i and i + M , defined as For M = 3, the three-tangle can be expressed in the following form [55]: where ) is the one-tangle of the electron (A is the system of the electron partitioned from the space of the two nuclei represented by systems B and C), and (similar definition holds for τ A|C ), with ρAB = σ ⊗2 y ρ * AB σ ⊗2 y .τ A|B (or τ A|C ) is also known as the entanglement of formation and is alternatively given by τ where λ j are the square roots of eigenvalues (in decreasing order) of ρ AB ρAB , with ρ AB = Tr C (ρ).
An alternative way to write the three-tangle (as in Ref. [71]) is: This is a complicated expression and averaging this quantity to find the M -tangling power of an arbitrary evolution operator for M = 3 is a difficult task.However, for CR-type evolution operators we find that it holds ⟨ψ|1 ⊗ σ y ⊗ σ y |ψ * ⟩ = 0 = ⟨ψ|σ z ⊗ σ y ⊗ σ y |ψ * ⟩.Considering ⟨ψ|1 ⊗ σ y ⊗ σ y |ψ * ⟩ we can show that this vanishes for CR-type evolutions: where we have used the fact that R † n j σ y (R n j ) * = σ y as well as the property ⟨ϕ|A|ϕ⟩ = ⟨ϕ|σ y C|ϕ⟩ = ⟨ϕ|σ y |ϕ * ⟩ = 0, where A is an anti-linear operator whose expectation value vanishes for all states |ϕ⟩ ∈ H and C denotes complex conjugation.Analogously, one can prove that ⟨ψ|σ z ⊗ σ y ⊗ σ y |ψ * ⟩ = 0. Thus, these results simplify the expression of the three-tangle for states generated by CR-type evolutions, meaning that we can express it as: We can linearize the above formula by noticing that we can write the three-tangle by extending the Hilbert space into a 6 qubit system: The above expression can be understood as a "vectorized" form of Eq. ( 22), since we note that |σ y ⟩ = vec(σ y ) = −i(|01⟩ − |10⟩), which means P (−) = 1/2(1−SWAP) between the sectors i and i+M can be represented as P (−) = 1 2 |σ y ⟩⟨σ y |.We also found that Eq. ( 23) can be generalized to odd M > 3, for arbitrary initial (pure) product states evolved under a CR-type evolution: The most general expression for the odd Mtangle for states that undergo arbitrary unitary evolution can be found in Ref. [57].We also verified by numerical inspection that Eq. ( 24) agrees with the definition of the odd tangle of Ref. [57] when the evolution is CR-type.
Averaging over the uniform distribution in H M 2 , we find that the M -tangling power reads: where i,i+M , with d = 2 and P (+) i,i+M = 1/2(1 2M ×2M + SWAP i,i+M ).To prove the expression for Ω p0 , it suffices to note that the uniform distribution of product states factorizes, and hence we can consider the total average as averages of the sectors i and i + M , ∀i.Thus, we have Ω p0 = M i ω i,i+M .Since Ω p0 is symmetric under the exchange of i and i + M , we then have ω i,i+M = cP (+) i,i+M , where c is a constant equal to c = (d + 1) −1 (see also Ref. [72]).Alternatively, the combined integral of the sectors i and i + M can be expressed in terms of an integral over the unitary group where the initial state is fixed and we vary the unitary acting on i and i + M systems: where we have used the fact that Tr[ρ 0 ] = 1 and Tr[ρ 0 SWAP i,i+M ] = 1.
Note that for even M , the M -tangling power holds for arbitrary gates U , whereas for odd M it holds only for controlled evolutions of the form

A.3 Proof of M -tangling power for CRevolution
Before going into the proof, it will be useful to find the action of products of symmetric/antisymmetric projectors (which project onto the sectors i and i + M ) on a trial product state.
Suppose we have a collection of states {|ϕ l ⟩, |ϕ l+M ⟩} M l=1 .Suppose further that we have the ordered ket If we act with the projectors i,i+M , on the above state we find: Note that for the above equality to hold it is important that both the kets and the symbols inside the kets are labeled by an index; otherwise if the kets are not labeled the above expression holds only up to re-ordering of the unlabeled kets.
To derive a closed-form expression for the Mtangling power of CR-type gates in defect systems, we start by considering odd M .For odd M , we start from the general expression: and we define the vector |m; The action of the projectors l,l+M on the above vector yields: Next, we define: and where we have suppressed the symbol of summation.In the above expressions we have defined , for l ̸ = 1.Using Eq. ( 31) and Eq. ( 30), we evaluate 2 M (U † ) ⊗2 P |m; m + M ⟩: (33) We then find the action of [2(d + 1)] M U ⊗2 Ω p0 on the bra ⟨m; m + M |: Thus, we find that where we have used the property {m i },{m i+M } i = i {m i },{m i+M } .We now work with the first expression of the Kronecker delta's which gives: ), (36) where we have used the fact that , as well as 4 − (Tr[Q The second term with the Kronecker deltas reads: For similar reasons the third term of Kronecker deltas also vanishes.Working with the last term of Kronecker deltas we find: ), (38) where we made use of (−1) M −1 = 1 since M − 1 is even.Therefore, our final expression for odd M is: ), (39) which concludes our proof.
For even M , note that P = M i=1 P (−) i,i+M .Again, we follow a similar procedure and start by calculating the action of 2 M (U † ) ⊗2 P from the left on the ket |m; m + M ⟩, which gives: Now, we combine Eq. ( 34) with Eq. ( 40), to obtain First, we focus on the first term of Kronecker deltas, which gives: ). ( It is further easy to show that the second and third terms of Kronecker deltas in Eq. ( 41) evaluate to 0, for similar reasons as in the odd case.The fourth term of Kronecker deltas produces: ). ( Therefore, the expression 2 M Tr[U ⊗2 Ω p0 (U † ) ⊗2 P ] for even M reads:

A.4 M -tangling power of non-unitary quantum evolution
In this section we consider a setup where we have L nuclear spins with K of them being the target ones, and L−K being the unwanted spins.Thus, the target subspace has M = K + 1 qubits, consisting of the electron and K target nuclear spins.In Ref. [38] we showed that the Kraus operators associated with the partial trace operation of the L − K unwanted spins have the closed form expression: with f (i) j given by: where for the i-th Kraus operator, the first product is taken over the unwanted spins comprising the environment that are in |0⟩, and the second product is taken over the unwanted spins that are in |1⟩.If all spins are in |0⟩, then the second product is 1, whereas if all spins are in |1⟩ the first product is 1.Further, we have The evolution of the target subspace is thus described by the Kraus operators via the quantum channel To include the impact of the unwanted subspace into the M -way entanglement that is generated by the quantum evolution, we define the M -tangling power of the non-unitary quantum channel: Clearly, in the limit of a unitary channel i.e., E r = E s = U (the summation is over the unique Kraus operator) we recover the M -tangling power of a unitary.
To derive a closed-form expression, we follow a similar procedure as in the unitary case.We begin by defining Starting with the odd M case, we find that where we have made use of the fact that only the first and last terms of products of Kronecker deltas survive.Therefore we find that ϵ p,M (E) for the odd case is: ). ( whereas for the even case following similar steps we obtain: ).
(52) In the above expressions we have made use of the completeness relation i.e.: and hence we can write: We note that the expressions for even/odd number of qubits in the target subspace are identical.
Further, the non-unitary M -tangling power is upper bounded by the unitary M -tangling power.
Let us now simplify the expression that appears due to the summation over Kraus operators.We first define the set S = {1, 2, . . ., L − K} of unwanted spin indices.Let P (S) be the power set of S, i.e.: This power set is useful in order to write down which unwanted spins are in |0⟩ or |1⟩ in the Kraus operators.Let p ∈ P(S) be an element of the power set and l ∈ p be an unwanted spin index within the subset p.We now start from: In the last line we have defined α l = ⟨0|(R n 1 |0⟩.To evaluate the above expression, let us look into some examples.If we have one unwanted nuclear spin, then P(S) = {{}, {1}}, and so it holds: If we have 2 unwanted nuclear spins, then we find: Similarly, if we have 3 unwanted nuclear spins, the summation over the power set evaluates to: operators (neglecting global phases): We label the Kraus operators of one unit with the bitstring m = (m 1 , m 2 , m 3 ), where each m j can take the value of 0 or 1.The only non-zero terms occur for (q, j, r) = (0, 0, 1) and (q, j, r) = (1, 1, 0), respectively.Thus, the Kraus operator reads: where we have defined: Also, we define ⊕ as addition modulo 2. It is easy to verify that the Kraus operators for one decoupling unit satisfy m 1 ,m 2 ,m 3 [D CPMG deph,m ] † D CPMG deph,m = 1 as follows: where we have used the fact that λ 0 + λ 1 = 1.
For multiple decoupling units, the Kraus opera- , where N is the number of iterations of a single unit.In this notation, each m (k) consists of indices (m 3 ), and each one of them takes the value 0 or 1.Each Kraus operator is distinguished by a different ordering of 0s and 1s in the total vector m.For N iterations, we have 2 3N such Kraus operators.For brevity of notation, we will label the composite Kraus operators that result after N CPMG iterations as D N,q , where q = 1, ..., 2 3N labels which Kraus operator we are referring to out of all Kraus operators.
Let us now look into the M -tangling power, where we focus on the subspace of the electron and the K target nuclear spins, and we ignore contributions from unwanted nuclei.This target subsystem now evolves under the non-unitary channel E deph (ρ) = r D N,r ρ[D N,r ] † .We first start by defining D N,r = j g(r) j σ jj ⊗ l R (l) n j .Note, that here R (l) n j is the rotation after N iterations, acting on the l-th nuclear spin.Thus, we can express the non-unitary M -tangling power of the channel as: 1 ) * g(r) 1 (g ), where we have made use of the results we found in the previous section.We have also defined 1 ) * g(r) 1 (g (s) 0 )].Note that in this case, the entangling power expresses only the evolution in the target subspace.In other words, we do not include the correlations that arise from the unwanted spins.In the limit of no dephasing, i.e., λ 0 = 1, λ 1 = 0 and θ = 0, the only Kraus operator that survives corresponds to the bitstring of m = [0, . . ., 0] T , and we recover the unitary channel.We can verify this from the above equation, since we will have r = s = 0 and g(r) 0 = g(r) The expression for g(r) j is defined as follows: Using the above definitions and considering only 1 CPMG iteration, we find the following expression for the sum: For N CPMG iterations we can show that R is given by: In the first line we defined f (r) = (g (r) 0 ) * g(r) 1 .Using the above result, we find that the expression of the M -tangling power of the target subspace in the presence of dephasing errors is: ). (72) It is easy to verify that in the limit of λ 0 = 1, (λ 0 = 0) λ 1 = 0 (λ 1 = 1) we get the unitary M -tangling power as expected.
If we are composing p CPMG sequences of different unit times t j and iterations N j then we will have: Based on the above expression we verify that in the limit where all p sequences have the same time t = t j , ∀j ∈ [1, p], we recover the previous expression for R. In particular, the above expression is the one we use for the sequential protocol when the electron undergoes dephasing.

C Comparison of sequential with multi-spin schemes in the presence of electronic dephasing
In this section, we compare the performance of the multi-spin with the sequential scheme in the presence of electronic dephasing errors.In Fig. 12(a) we show ϵ p,M (E deph ) of case "16" of preparing GHZ 3 -like states for the sequential scheme.In Fig. 12(b) we show again ϵ p,M (E deph ) of case "7" of preparing GHZ 3 -like states for the multi-spin scheme.The nuclei we control with the sequential scheme are C5 and C13, and with the multi-spin scheme, C4 and C5.The total sequence time for the sequential scheme is 1357.9µs, and for the multi-spin is 1916.6 µs.
Although the multi-spin scheme requires slightly longer time than the sequential one, it is again more robust to dephasing errors.This feature could be attributed to the different unit times t j and sequence iterations N j that enter the expression of ϵ p,M (E deph ) for the sequential scheme since we are composing CPMG sequences of different timings and iterations in the sequential protocol.

D Comparison of XY2 and CPMG performance to pulse errors
In this section, we study over-/under-rotation errors for the electronic control pulses.We consider the optimal cases we found in the main text for preparing GHZ 4 -like states via the sequential or multi-spin scheme.The systematic pulse errors are simulated as R x/y (π + ϵ) = e −iπ/2(1+ϵ)σ x/y for the CPMG and XY2 sequences.We calculate the deviation between the M -tangling power of the target subspace in the absence of errors, ϵ p,M (U ) ϵ=0 , with the one in the presence of errors, ϵ p,M (U ) ϵ̸ =0 .We define this deviation as ∆ϵ p,M (U ) = |ϵ p,M (U ) ϵ=0 − ϵ p,M (U ) ϵ̸ =0 |.In Fig. 13(a), we show ∆ϵ p,M (U ) for the 7 cases of preparing GHZ 4 -like cases via the multi-spin scheme, assuming the XY2 sequence.The different colored lines correspond to each one of the 7 cases, and the labels correspond to the number of total sequence iterations.We see that for small pulse errors (≤ 3 − 4%), the deviation in estimating this quantity is on the order 10 −3 − 10 −2 .Additionally, we expect that as the number of pulses increases, the deviation is, in principle, enhanced [e.g., see dark red line versus green line].This feature, however, does not always hold since the dynamics are case-dependent and tend to corrupt the GHZ preparation differently.In Fig. 13(b), we repeat the same calculation assuming the CPMG sequence for a range of systematic pulse errors from 0.1 − 0.8%.As expected, the ability of the CPMG sequence to prepare high-quality GHZ 4 -like states is severely impacted by the pulse errors.
In Fig. 13(c), we display in four panels the deviation in M -tangling power for the 29 cases of preparing GHZ 4 -like states using the sequential protocol.The labels correspond to the total number of iterations, 3 j=1 N j .In this scenario, we are assuming the XY2 sequence.We notice a similar performance as for the multi-spin scheme, and the results again reveal the robustness of XY2 under moderate pulse errors (≤ 3 − 4%).

E Rotation angle and rotation axis errors
In this section we consider the effect of rotation angle and rotation axis errors on the M -tangling power.We consider as an example the 7 cases of preparing GHZ 4 -like states via the multi-spin scheme.To showcase the robustness of the XY2 decoupling sequence, we consider both rotation axis and rotation angle errors, assuming the estimated parameters from Ref. [69], which we also mentioned in Sec.7.2.In particular, we model the erroneous π X and π Y pulses as: In Fig. 14 we show the ideal M -tangling power (dark blue bar) assuming a perfect CPMG decoupling sequence.The light blue bars correspond to the XY2 decoupling sequence which includes the errors ϵ x = ϵ y = −0.02,q z = 0.05, q x = 0.005.We observe that the deviation in the expected M -way correlations in the presence of both rotation axis and rotation angle errors is sufficiently small, and hence the XY2 decoupling sequence can be reliably used in experiments to prepare high-fidelity GHZ states.

F Uncertainty in HF parameters
In this section, we consider the effect of errors due to uncertainty in the experimentally measured HF parameters of the nuclear spins.We use the same optimal sequence parameters we found for the HF parameters we considered in the main text and additionally shift the HF parameters A and B of the target nuclear spins by 0.01 or 0.05 kHz.In Fig. 15, we consider a 0.01 kHz shift and plot the deviation ∆ϵ p,M (U ) := ϵ p,M (U )−ϵ uncertain p,M (U ) for the cases of preparing GHZ states up to 6 qubits for the sequential scheme in Fig. 15(a) and for the multi-spin scheme in Fig. 15(b).We observe that the presence of a 0.01 kHz uncertainty in the HF parameters only produces a negligible uncertainty on the order of 10 −3 relative to our ideal calculations for ϵ p,M (U ).In Fig. 16, we repeat the same calculations for the sequential scheme in Fig. 16(a) and for the multi-spin scheme in Fig. 16(b), assuming a shift of the HF parameters by 0.05 kHz.We note that the error in the estimation of ϵ p,M (U ) is on the order of 10 −2 − 10 −3 .Thus, based on the experimentally achievable accuracy in the measured HF parameters [37], we can estimate the M -tangling power to a relatively narrow confidence interval.Optimization within a small range of the optimal parameters t * and N * that we find for definite values of HF parameters can guide the experimental calibration of the timing of the pulses and optimal sequence iterations.

G Optimization process for GHZ M -like generation
Here we highlight the optimization procedure for generating GHZ M -like states using the sequential or multi-spin protocols.In Fig. 2 of the main text, we show a diagram of the optimization procedure for the sequential protocol.First, we need to provide as inputs: • the size of the GHZ-like state we want to prepare (|GHZ|), • the maximum resonance number (k max ) to search for each nuclear spin, • the maximum time for the total sequence and of individual gates (T max ), • the tolerances of unwanted/target nuclear one-tangles, • the gate error tolerance.
Then, for each nuclear spin, we search over all resonances by varying the unit time around each resonance by δt = ±0.2µs, and with a time step of 5 × 10 −3 µs and select as optimal the unit time t where n 0 • n 1 is as close as possible to the value −1.Having fixed the unit time, we then use the analytical expressions for the minima of G 1 which give us the number of iterations that maximize the nuclear one-tangle [see Ref. [38] for the expressions of the minima].Since G 1 is periodic, there are multiple numbers of iterations that can minimize G 1 , and so we usually truncate to about 15 maxima of the one-tangles, which we also postselect such that N • t ≤ T max .Note that here T max could be the gate time restriction we impose for the gate time of a single decoupling sequence, rather than the total gate time restriction for the composition of all sequential gates.We then inspect all elements in the sets of {t, N } and keep those that give a sequence with gate time smaller than T max .After this step, we calculate the remaining nuclear one-tangles and check if they are smaller than the unwanted one-tangle tolerance.If this is not satisfied, we reject that particular (t, N ) case, whereas if this is satisfied, we store the unit time and number of iterations.We repeat this process for each nuclear spin, such that we have all the possible unit times and iterations that can give maximal entanglement of the target spins with the electron while keeping the unwanted one-tangles minimal.At this stage, we have multiple unit times and iterations that satisfy these requirements.
We then combine the unit times and iterations corresponding to sets associated with the selection of any |GHZ| − 1 nuclei out of the entire nuclear register, such that the total decoupling sequence does not exceed the maximum time, T max .For each such nuclear spin combination, we evolve all nuclei individually under the composite evolution and obtain their nuclear one-tangles.At this step, each nuclear spin combination is associated with multiple unit times and iterations that we could consider, so we need to choose which t and N we keep for each distinct nuclear spin combination.We choose those t and N which give rise to a maximal M -tangle at the end of the composite |GHZ|-1 entangling gates.Different choices could be made here, e.g., selecting the t and N that give the shortest gate time or penalizing both gate time and the deviation from the maximal possible M -tangle using an appropriate cost function.At this stage, we have narrowed down nuclear spin candidates and associated each nuclear spin combination with a particular sequence composed of times t j and iterations N j that gives rise to maximal M -way entanglement.For each of these cases, we inspect whether the unwanted nuclear one-tangles of the composite evolution are below the unwanted one-tangle tolerance we imposed.If this is satisfied, we accept this case and proceed with the calculation of the gate error.Finally, we accept this case if the gate error is lower than    Table 3: Tolerances (tol.) and relevant parameters for the optimization of the generation of GHZ M -like states using the multi-spin protocol.
the gate error tolerance.In Table 2, we provide the tolerances and relevant parameters we set for each GHZ size.
Regarding the multi-spin protocol, we are not composing gates since we are interested in a single-shot operation that generates direct entanglement of the electron with multiple nuclei.We explain here the procedure we follow according to the flowchart of Fig. 2. We start with one nuclear spin selected from the entire register and find its resonance time for some particular resonance number k.We vary the unit time t within ±0.25 µs and define an upper bound for the maximum number of iterations that respects the gate time restriction for the particular unit time.Then, for all unit times and N ∈ [1, N max ] we obtain all the nuclear one-tangles using the knowledge of rotation angles and the dot product of the nuclear axes for each nucleus given one iteration.Then, for all possible times and iterations, we check how many one-tangles are above the target one-tangle tolerance and how many are below the unwanted one-tangle tolerance.If we have |GHZ|−1 one-tangles above the target one-tangle tolerance, and all other nuclear one-tangles are below the unwanted one-tangle tolerance, then we accept this case.We repeat this process by choosing a different resonance number k, up to some k max .Then, we select a different nuclear spin from the register as our starting point and repeat all the aforementioned steps.
After completing the above stage, we have multiple times and iterations for which we maximize |GHZ| − 1 target one-tangles simultaneously.We then rearrange all these cases corresponding to different unit times and iterations in terms of unique spin combinations.For each spin com-bination, we select the time and iterations of the single-shot operation, which give the maximal M -tangle.Similar to the sequential scheme, the choice of narrowing down particular t and N could be based on minimal gate time or a cost function that both minimizes gate time and maximizes the M -tangle for the particular spin combination.Finally, for each spin combination we evolve all unwanted nuclear spins individually, such that we pass this information to the calculation of the gate error.If the gate error is below the tolerance we imposed, we then accept this case.

H One-tangles of Sequential Scheme for generation of GHZ M -like states
Here we present the nuclear one-tangles corresponding to the sequential entanglement scheme of Fig. 5.The nuclear one-tangles for all different realizations of entangling M − 1 nuclei with the electron to prepare GHZ M -like states are shown in Fig. 17.The labels above each bar of each case refer to some 13 C nuclear spin of the register labeled as Cj, with j ∈ [1,27].

I Eigendecomposition of mixed state
To keep the discussion general, we consider L nuclear spins, with K of them being the target ones and hence, L − K being the unwanted nuclei.To derive the optimal ensemble for the calculation of entanglement of the mixed state, we make the assumption that the initial state of the system is any arbitrary product pure state: (77) Thus, the full density matrix reads: ρ = j,k∈{0,1} σ jj ρ el σ kk ⊗ l R (l) n j ρ (l) nuc (R † n k ) (l) .(78) Next, suppose that we wish to trace out the last L − K spins from the density operator (assuming that we have ordered the basis such that these appear last in the Kronecker product): Tr[R (l+K) (80) Thus, we can easily find the eigenvalues of the electron's reduced density matrix which read: (81) whereas the eigenvectors are given by: with c = α * βf 10 (provided c ̸ = 0).Now, to find the eigenvectors of the total reduced density matrix, we need to reapply the controlled gates of the target subspace on the matrix whose columns are the two eigenvectors: (83) Thus, we now have the diagonalized density operator: which means that we can use the trial state: to find the entanglement of the mixed reduced density matrix of the electron and the target nuclei.Note that given arbitrary initial states, the highest rank that the reduced density matrix can have (irrespective of the initial pure state of the system, the number of total nuclei, or the number of nuclei we trace out) is 2, due to the form of the controlled evolution operator.
In the case when the electron starts from a nonsuperposition state e.g., when β = 0, we find that the two eigenvalues are (86) In this scenario, since λ − = 0, the reduced density matrix corresponds to a pure product state.

J Concurrence and three-tangle
Here we compare the entanglement of rank-2 mixed states involving two qubits with the entanglement of rank-2 mixed states involving three qubits.For two-qubit pure states the concurrence is defined as ), where ρ A(B) is the reduced density matrix of system A (B).The entanglement of a mixed two-qubit state can be computed with similar methods as in the main text by minimizing the concurrence of a trial state.We assume that we have a trial state of the form: where |Φ ± ⟩ = 1/ √ 2(|00⟩ ± |11⟩) are the two Bell states.We further consider a rank-2 mixed threequbit state which is the mixture of two |GHZ ± ⟩ = 1/ √ 2(|0⟩ ⊗3 ±|1⟩ ⊗3 ) states.Thus, similarly to the main text, we define the trial state: (88) For each value of p ∈ [0, 1] we vary the relative phase of the two terms, χ, for both |ψ trial,2 ⟩ and |ψ trial,3 ⟩, and obtain the minimal concurrence, or three-tangle, respectively.We plot the results as a function of p in Fig. 18.We observe that the minimum concurrence, as well as the minimum three-tangle are both convex functions.In fact, the minimum concurrence is given by C = 2|p − 1/2| and the minimum three-tangle by (1 − 2p) 2 .We note that for all values of p (except for p = 1/2 where we have maximal mixed states) the minimum three-tangle is lower than the minimum concurrence, and hence the threeway entanglement is more sensitive to the relative ratio of the superposition terms in the mixed state than the concurrence of two-qubit mixed states.

J.1 Parameters for 27 nuclear spins
The HF parameters for the 27 nuclear spin register we consider in the main text are presented in Table 4 and can also be found in Refs.[37,61].

Figure 1 :
Figure 1: Schematics of two protocols for generating GHZ M -like states (shown for M = 4 and using the CPMG sequence).(a) The multi-spin scheme is capable of generating direct entanglement between the electron and a subset of nuclei from the nuclear spin register in a single shot.(b) The sequential scheme requires M − 1 consecutive entangling gates to prepare a GHZ M -like state.

Figure 3 :
Figure 3: Generation of GHZ 3 -like states via the sequential protocol.Each case # corresponds to a different composite DD sequence that sequentially entangles a pair of nuclei with the electron to generate GHZ 3 -like states.(a) Nuclear one-tangles (scaled by 2/9) after the two entangling gates.The first (second) bar of each case corresponds to a particular nuclear spin, labeled as Cj.(b) Nuclear rotation angles (ϕ := ϕ 0 ) and (c) dot products of nuclear rotation axes at the end of the composite evolution of the two entangling gates.(d) M -tangling power for each case of composite evolution.(e) Total gate time of the sequence and (f) gate error of the composite evolution due to residual entanglement with the remaining nuclear spins that are not part of the GHZ 3 -like state.

Figure 4 :
Figure 4: Generation of GHZ 3 -like states using the multi-spin protocol.Each case # corresponds to a different DD sequence that simultaneously entangles a pair of nuclei with the electron to generate GHZ 3 -like states.(a) Nuclear one-tangles (scaled by 2/9) after the single-shot operation.(b) Nuclear rotation angles and (c) dot products of nuclear rotation axes at the end of the single-shot entangling operation.(d) M -tangling power for each case.(e) Total gate time of the sequence and (f) gate error of the evolution due to residual entanglement with the remaining nuclear spins that are not part of the GHZ 3 -like state.

Figure 5 :Figure 6 :
Figure 5: Preparing states locally equivalent to GHZ M via the sequential scheme.Each case # corresponds to a unique composite DD sequence that selects M − 1 nuclei from the register.(a) M -tangling power of the composite M − 1 entangling gates for various cases.(b) Gate error of the composite evolution due to the presence of unwanted nuclei, and (c) gate time of the total sequence for each case and each different size of GHZ-like states.The panels from top to bottom correspond to GHZ 4 , GHZ 5 , GHZ 6 , GHZ 7 , GHZ 8 , GHZ 9 and GHZ 10 states.The nuclear onetangles for each case and value of M can be found in Appendix H.

Figure 7 :
Figure 7: (a) Three-tangle of the pure state of the target subspace ignoring unwanted nuclei (purple), and threetangle of either eigenvector (blue) of the mixed state for each of the 39 cases of Fig. 3.The 39 cases are sorted in terms of largest eigenvalue, λ + = p, shown along the x-axis.(b) Three-tangle of the mixed electron-nuclear state (red), obtained by minimizing the entanglement of the trial state of Eq. (9) for each value of p.The dashed line shows the convex hull.

Figure 8 :
Figure 8: (a) Comparison of the unitary, ϵ p,M (U ), with the non-unitary, ϵ p,M (E), M -tangling power for the 29 cases of Fig. 5 that prepare GHZ 4 -like states in the target subspace, via the sequential protocol.(b) Comparison of ϵ p,M (U ) with ϵ p,M (E) for the 7 cases of Fig. 6 for preparing GHZ 4 -like states via the multi-spin protocol.In both panels, the dark blue bars show the unitary M -tangling power, the pink bars the non-unitary M -tangling power, whereas the green bars show an approximation of the latter.

Figure 9 :
Figure 9: M -tangling power of the target subspace when the electronic qubit experiences dephasing.(a) M −tangling power for the 39 cases of preparing GHZ 3 -like states with the sequential protocol.(b) M -tangling power of the 27 cases of preparing GHZ 3 -like states with the multi-spin scheme.The darker blue bars in (a) and (b) correspond to the no-dephasing scenario, whereas the remaining colors include electronic dephasing with dephasing angles θ −1 = (400, 100, 50) µs.For the cases where we consider dephasing, we set λ 0 = 0.98.(c) M -tangling power for case # 1 of the sequential scheme for different dephasing angles and λ 0 parameters.(d) Same as in (c) but for case # 1 of the multi-spin scheme.

Figure 10 :
Figure 10: Generation of GHZ 4 -like states in the presence of systematic over-rotation pulse errors.(a), (b) Mtangling power using the 29 cases we found for preparing GHZ 4 -like states via the sequential scheme, assuming a systematic error of 2% in (a), and of 8% in (b).(c), (d) M -tangling power for the 7 cases of preparing GHZ 4 -like states using the multi-spin scheme assuming a 2% error in (c), and an 8% error in (d).In all panels, the bars with the highest value of ϵ p,M (U ) correspond to the error-free case.The orange bars show the CPMG including pulse errors, whereas the light blue bars show the XY2 including pulse errors.All bars are scaled by the maximal value of (2/3) 4 .

Figure 11 :
Figure 11: Generation of GHZ 4 -like states in the presence of random over-/under-rotation pulse errors.(a) Mean M -tangling power using the 29 cases we found for preparing GHZ 4 -like states via the sequential scheme, assuming rotation angle errors sampled from a normal distribution with standard deviation σ = 0.01.(b) Mean M -tangling power for the 7 cases of preparing GHZ 4 -like states via the multi-spin scheme assuming the same standard deviation.In both (a) and (b), the bars with the highest value of ϵ p,M (U ) correspond to the error-free case.The orange bars correspond to the CPMG case including pulse errors, whereas the light blue bars show the XY2 case including pulse errors.All bars are scaled by the maximal value of (2/3) 4 .Each error bar for CPMG and XY2 in (a) captures the range of ϵ p,M (U ) where we run 500 trials per case, with different random samplings per sequence iteration.In (b), we perform 1000 random trials per case to display the error bars.

Figure 12 :Figure 13 :
Figure 12: M -tangling power of the target subspace when the electron undergoes dephasing.(a) ϵ p,M (E deph ) for case # 16 for preparing GHZ 3 -like states via the sequential scheme as a function of the dephasing angle θ and λ 0 .(b) Same as in (a) for case # 7 of the multi-spin scheme.

Figure 14 :
Figure 14: Robustness of XY2 decoupling sequence in the presence of pulse and rotation axis errors.The dark blue bars show the ϵ p,M (U ) we find via the ideal CPMG sequence for the 7 cases of preparing GHZ 4 -like states via the multi-spin scheme.The light blue bars show the ϵ p,M (U ) obtained via the XY2 decoupling sequence assuming both rotation axis and rotation angle errors.

Figure 17 :
Figure 17: One-tangles of nuclear spins after composing M − 1 sequential entangling gates to create the cases of Fig. 5.Each case corresponds to a distinct DD sequence that selects different nuclei from the nuclear spin register.Each bar shows the one-tangle of a particular nuclear spin, and the text above the bars corresponds to the nuclear spin labels.From top to bottom, we show the nuclear one-tangles for the generation of GHZ 4 -like up to GHZ 10 -like electron-nuclear entangled states.

Figure 18 :
Figure 18: Minimum concurrence of a mixture of Bell states (blue), and minimum three-tangle of a mixture of |GHZ ± ⟩ states (red).The black dashed line is the analytical expression of the minimum three-tangle.

Table 2 :
Tolerances (tol.) and relevant parameters for the optimization of the generation of GHZ M -like states using the sequential protocol.

Table 4 :
Hyperfine parameters of the13C atoms we considered in the paper.