A many-body singlet prepared by a central spin qubit

Controllable quantum many-body systems are platforms for fundamental investigations into the nature of entanglement and promise to deliver computational speed-up for a broad class of algorithms and simulations. In particular, engineering entanglement within a dense spin ensemble can turn it into a robust quantum memory or a computational platform. Recent experimental progress in dense central spin systems motivates the design of algorithms that use a central-spin qubit as a convenient proxy for the ensemble. Here we propose a protocol that uses a central spin to initialize two dense spin ensembles into a pure anti-polarized state and from there creates a many-body entangled state -- a singlet -- from the combined ensemble. We quantify the protocol performance for multiple material platforms and show that it can be implemented even in the presence of realistic levels of decoherence. Our protocol introduces an algorithmic approach to preparation of a known many-body state and to entanglement engineering in a dense spin ensemble, which can be extended towards a broad class of collective quantum states.


INTRODUCTION
Controlling quantum properties in many-particle systems, whether for technological advantage or foundational studies, can be reduced to controlling the relative participation and phase of the system's eigenstates. In most cases, this leads to entanglement amongst the system's particles [1]. The initialization of a quantum system to a pure state is a necessary starting point from which to engineer entanglement [2]. Initialization through traditional cooling techniques brings a quantum system in contact with a bath whose temperature is below that of the system's characteristic energy scale, bringing the target system to its ground state [3]. An equivalent picture exists for driven systems, where directionality within an energy hierarchy of dressed states can bring the system to an effective ground state of the dressed system [4]. The latter approach is more versatile as it allows the design of tailored ground states. Crucially, it also requires access to the degree of freedom one intends to cool. Versions of this approach appear in driven-dissipative state preparations of photonic systems [5,6], diamond color centers [7][8][9], epitaxial quantum dots [10,11], and trapped atoms in an optical resonator [12].
Central-spin systems, typically consisting of a single electronic spin coupled to an ensemble of nuclear spins [13], have been a particularly potent testing ground for the active approaches to state preparation [14,15]. Firstly, this is by necessity: nuclear-spin energy scales make ground state preparation well beyond the reach of modern refrigeration techniques. Secondly, the one-to-all coupling of the central spin qubit to an ensemble of spins yields a convenient proxy for the ensemble spins [16]. In the few-spin regime of diamond color centers or rare earth ions, dynamic nuclear polarization can initialize a set of proximal spins to a high-purity polarized state [17]. This is then the starting point for qubit storage in an ensemble excitation [18], two-body singlet engineering [19], or spin-by-spin quantum computation [20].
In the limit of dense ensembles, the constituent spins are indistinguishable when interacting with the central spin. Entanglement is then readily generated and measured by controlling the central spin's dynamics [16,21,22]. Collective states of the ensemble become natural targets of a driven purification technique [15]. In addition, an effective all-to-all coupling mediated by the central spin [23] leads to a highly-correlated behavior of the ensemble which, in principle, can be harnessed for state engineering. Despite the absence of individual spin control, the dense systems of interest [14,24,25] operate in a truly many-body regime of 10 4 to 10 6 interacting ensemble spins -far exceeding the number of physical qubits present in near-term quantum simulators [26,27].
A many-body singlet is a superposition of ensemble spins with zero angular momentum. The singlet state is a hyperfine vacuum and renders the bath invisible to the central spin, protecting the latter from interacting with its environment [7,28,29] -a decoherence-free subspace for the central spin qubit. Owing to destructive interference between pairs of ensemble spins, the ensemble is also protected from noise whose length scale is larger than the system size. This makes it a particularly useful and readily detectable first target state to prepare in a arXiv:2301.10258v1 [quant-ph] 24 Jan 2023 central-spin system. Further, many-body singlet preparation offers insights into the evolution of entanglement under slower intra-bath interactions, as well as anomalous spin diffusion at longer time scales [30].
To this date, state engineering efforts in dense central spin systems remain confined to tuning the mean field degrees of freedom, such as ensemble polarization and its fluctuations [15,21,31,32], and classical correlations amongst ensemble spins [22]. A degree of purification has been achieved via polarization of 80% through optical techniques [33] and nearly 50% via central-spin control [22] -such approaches proving to be generally challenging due to a dependence of the central spin transition frequency on the polarization of the bath. Meanwhile stabilizing the ensemble polarization via the central spin [15] can approach the quantum limit of polarization stability. Despite these efforts, the resulting nuclear states remain highly mixed, featuring little inter-particle phase coherence. A degree of quantum correlation among spins was observed when probing a partially polarized ensemble via the central spin [22], which suggested the possibility of purification via reduced total angular momentum states, so-called dark states [34]. Theoretical proposals to engineer the state of dense central spin systems have focused on dissipative phase transitions [4], or quantum memory [10,34] and spin squeezing schemes [35] relying on fully polarized initial states. A protocol for direct control over the ensemble's inter-particle phase, allowing for purification and entanglement engineering at low total ensemble polarization, is still missing.
In this work, we use a simple and realistic form of symmetry breaking to gain control over an ensemble's inter-particle phase. Leveraging this control, we propose a three-stage state engineering protocol that utilizes the effective interaction between two distinct spin species, naturally present in real physical systems [14,24,25], and proceeds via control of the central spin exclusively. The first stage of the protocol locks the total polarization of the system to zero [15]. The second stage initializes the two spin species to an anti-polarized state with near-unit purity. The third stage involves a sequence of unitary gates that drives the system into a many-body singlet via phase engineering. Having in mind a near-term experimental realization, we quantify the protocol robustness as a function of model parameters and identify candidate physical platforms in which it could be successfully implemented.

Symmetry and the total angular momentum representation
Central spin systems feature high-dimensional Hilbert spaces. In the dense limit, the central spin's interaction with the ensemble does not distinguish individual spins. This leads to collective symmetries and corresponding constants of motion, which we focus on here.
In this simplest scenario of a perfectly homogeneous spin bath (see Fig. 1a), a general system Hamiltonian is unchanged under the re-ordering of ensemble-spin indices. Such invariance is the highest symmetry that a spin bath can exhibit, and it results in the emergence of a constant of motion that dictates the rate of collective dynamics, as in Dicke superradiance [36]. Within the bath of spin-1/2 particles, this constant is identical to the magnitude of the total angular momentum, I, directly related to the eigenvalue of: where i k is the single spin operator of the k-th spin in the ensemble, and N is the total number of ensemble spins. Ensembles of higher-spin particles, such as spin-3/2, can be described by equivalent, albeit more numerous, constants of motion. Under this symmetry, coupling to the central spin cannot change the magnitude of the total angular momentum, only the polarization, and coherent control over the ensemble is thereby limited. In particular, efforts to prepare a many-body singlet (I = 0) are futile and the ensemble dynamics are governed by thermally (with β ∼ 0) occupied states of I, dominated by the highest degeneracy states near I ∼ √ N . We propose making use of the simplest form of reduced symmetry to gain control: breaking the system into two distinguishable but equally abundant ensembles (see Fig. 1b), which we call species. The groups of spins of the same species are characterized by their individual total angular momenta: I 1 and I 2 . Their magnitudes, I 1 and I 2 , become the new constants of motion. In the discussions to follow, we assign them the value N/2, which is the most likely value for a fully mixed initial state. Our results hold for all (I 1 , I 2 ) values in a ∼ N/2 vicinity [37], thus capturing all experimentally relevant dynamics. We also note that this easily extends to more than two species.
This situation of two distinguishable spin species makes it possible to alter the magnitude of the total angular momentum, I = |I 1 +I 2 |, and therefore to significantly reduce it (see Fig. 1c). When controlled via a directional pumping process (i.e. cooling), the spin ensembles can be initialized into a pure collective anti-polarized state [38]: with M = min(I 1 , I 2 ). Such a state contains no coherence between the species, and represents a classical limit to the total angular momentum reduction, featuring a noise of I 2 cl = (I 1 + I 2 )(|I 1 − I 2 | + 1) ∼ √ N ; a factor √ N lower than that of a thermal state. Creating entanglement via a controlled phase between the two species can lower the magnitude of the total an- FIG. 1. State engineering in the central spin system a, Perfectly homogeneous central spin system, symmetry-protected from changing its total angular momentum. b, Symmetry-broken central spin system, consisting of spin species I1 and I2, no longer protected from changing its total angular momentum. c, Reduction of the total angular momentum magnitude via the state purification protocol, resulting in a narrowed steady-state distribution (blue-shaded curve), peaked around I = N −1/4 . d, Effective Jaynes-Cummings ladder of states |I z 1 = −M + n, I z 2 = M − n , with anti-correlated I z 1 and I z 2 , where M = min(I1, I2). The thick and faint yellow arrows illustrate the dominant, and 2∆ω-detuned three-body interactions, respectively. The curvy orange arrows represent the central spin resets. The pale blue arrow displays a net direction of the phase space flow within this effective sideband-cooling process. e, Stages of anti-polarized state preparation. The red and blue arrows within the generalized Bloch spheres denote the total angular momenta of the two species, I1 and I2. Stage one (left panel): Locking of the total polarization, I z = I z 1 + I z 2 , to zero. Stage two: activation of the three-body interaction (middle panel), equivalent to the back-action sensing, followed by central spin reset (right panel). gular momentum further down to I = |I 1 − I 2 |. In particular, for I 1 = I 2 , the quantum limit is reached after the preparation of a many-body singlet: expressed uniquely in the |I z 1 = −M + n, I z 2 = M − nbasis using Clebsh-Gordan coefficients. The many-body singlet is characterized also by a full noise suppression: System Hamiltonian: control via the central spin We take the general Hamiltonian for a dense central spin system in an external magnetic field [13], and split the spin ensemble into two spin species i = 1, 2: We consider the regime of high magnetic fields, as defined by a dominant Zeeman interaction ∝ ω c of the central spin, S. The second term captures the internal energy structure ω i of the spin ensembles. The last term represents a Heisenberg-type hyperfine coupling (∝ a i < ω i ) between the central spin and the bath spins. We consider a symmetry breaking between the two spin ensembles ω 1 = ω 2 , which in real systems can take myriad forms. Without loss of generality [39], we assume a 1 = a 2 ≡ a and ω 1 > ω 2 . In this high field regime, the central spin quantization axis is pinned to the z-direction, and the hyperfine interaction reduces to [40,41]: The leading and only first-order hyperfine term is the collinear interaction (used in total polarization locking). The second term renormalizes the nuclear Zeeman interaction of the ensemble by a negligible amount of −a 2 /4ω c ω i . The third and fourth terms contain the effective two-body and the three-body interactions, respectively. The former will remain off-resonance during our protocol, while the latter will be critical in cooling and correlating the two spin ensembles.
Exclusively to stabilize the ensemble polarization [21], we also consider that, apart from the Zeeman interaction (∝ ω i ), the ensemble spins are subject to a small eigenstate-mixing interaction (∝ ν i ω i ) [31], which results in an effective non-collinear term: Finally, we consider a control of the entire system exclusively via the central spin, represented by a resonant drive of the central spin with strength Ω ω c , and work in a rotating frame of reference, under Rotating Wave Approximation [37].

Preparation of a pure anti-polarized state
Prior to the protocol execution, the spin ensemble is found in a fully-mixed state characterized by a density matrix ρ ∝ 1. The first two stages of our protocol initialize the ensemble to a pure collective state -an antipolarized state of the two spin species (Eq. 2) -which enables the generation of a many-body singlet of the whole ensemble in the third stage (last section of this article).
Stage 1: polarization locking -The initial state exhibits maximal uncertainty of the total polarization of the bath ∆ 2 I z ∼ √ N , where I z = I z 1 + I z 2 . The first stage of the protocol reduces this uncertainty to zero using a previously developed technique [15]. In brief terms, the collinear hyperfine interaction (the first term in Eq. 5) allows the central spin to sense the polarization-deviation from the I z = 0 lock-point which induces a Larmor precession around the z-axis. Subsequently, this deviation is corrected by a resonantly activated (Ω = (ω 1 +ω 2 )/2 -c.f. Ref. [42]) non-collinear interaction (using Eq. 6) which translates the acquired phase into a change in the total ensemble polarization. Finally, the state of the central spin is reset and this stage is repeated until the ensemble reaches the limit of I z = 0 and ∆ 2 I z = 0. From this point on, the non-collinear interaction remains offresonant.
Stage 2: full purification -At the end of the first stage, the ensemble's locked zero-polarization state is a mixture of states |I z 1 = −M + n, I z 2 = M − n where n = 0, 1, .., 2M , for which the remaining uncertainty lies in the polarization of individual species, I z 1 and I z 2 . The second stage of the protocol removes this uncertainty to produce a pure collective state. Doing so relies on the driven resonant activation Ω = ∆ω ≡ ω 1 − ω 2 of the threebody interaction (fourth term of Eq. 5). For the central spin initialized in state |↓ x , this interaction activates the transition: This can be visualized in an effective Jaynes-Cummings ladder of states ( Fig. 1d) parameterized by the principal quantum number n corresponding to the |I z 1 = −M + n, I z 2 = M − n state. Equation 7 then represents a single quantum n ↔ n − 1 transition down the ladder. Combined with a central spin reset |↑ x → |↓ x applied every half-period of the three-body interaction, repeating the transition and reset forms a directional pumping process -equivalent to sideband cooling in harmonic systems [43] -towards the ladder's ground state n = 0. This ground state consists of the fully antipolarized ensemble state |I z 1 = −M, I z 2 = M , which is a pure collective state of the ensemble.
Maintaining the directionality of the pumping process relies on a detuning of 2∆ω between the selected transition |I z 1 , I z 2 |↓ x ↔ |I z 1 − 1, I z 2 + 1 |↑ x and the unwanted transition |I z 1 , I z 2 |↓ x ↔ |I z 1 + 1, I z 2 − 1 |↑ x (thick and faint yellow arrows in Fig. 1d, respectively), which drive the system towards the opposite ends of the ladder. This is sustained as long as the three-body interaction strength, N a 2 /2ω c , remains smaller than ∆ω. We characterize the directionality with their ratio, which ultimately sets the purification limit: We summarize this two-stage cooling process in Fig. 1e, where the initial locked state has two spin ensembles with arbitrary but opposite orientations. Semi-classically, the second stage (three-body interaction) is equivalent to the central spin sensing the Larmor precession beatnote between the two spin ensembles (recall ω 1 = ω 2 ) in the hyperfine back-action, which acts on the central spin akin to a classical driving field. The timing of the central spin reset favors a particular phase of this beatnote, and thus prepares the two ensembles in a specific orientation. , as a function of the protocol time. ρ0 denotes the state of the bath at the beginning of the second stage of the protocol. The x-axis is shared across the panels.

Ideal system dynamics
We verify the convergence of the second stage of the protocol towards the ground state by calculating the full quantum evolution of the system numerically. We treat the ideal case, for which the three-body interaction is fully coherent and the central spin reset is instantaneous. We repeat this second stage as many times as necessary to reach steady state. For convergence towards the steady state under reasonable computational resources, we choose to work with N = 32 bath spins and focus on the I 1 = I 2 = 32/2 manifold [37]. We take the first stage of the protocol to be fully capable of confining the ensemble's dynamics to the I z = 0 subspace [15], spanned by {|I z 1 = −M + n, I z 2 = M − n , n = 0, 1, 2, .., 2M }, to which we accordingly restrict our simulation. At the beginning of the simulation, we set the ensemble's den-sity operator, ρ, to an equal statistical mixture of the subspace basis states. Figure 2a shows the central spin's evolution under the activated three-body interaction and confirms the coherent oscillation in the {|↓ x , |↑ x }-basis, with a half-period of τ 0 ∼ 2πωc a 2 (I1×I2) -inversely proportional to the fastest three-body interaction rate [37]. The state of the central spin is instantaneously reset every τ 0 , and during the consecutive iterations the magnitude of the acquired |↑ x -population drops down to zero as the ensemble approaches its ground state. Semi-classically, this corresponds to a drop in the magnitude of the noise I 2 sensed by the central spin (Fig. 2a). We verify this independently as shown in Fig. 2b, where I 2 saturates to its classical ∼ √ N limit as steady state is reached. Figure 2c shows the polarizations of each of the two species, I z 1 and I z 2 as a function of iteration time. We see that they saturate to maximal and opposite values I z 1 = + min(I 1 , I 2 ) and I z 2 = − min(I 1 , I 2 ). Importantly, at the longest iteration times their uncertainties approach zero, as expected for a pure anti-polarized state.
We quantify the quality of state preparation following our protocol using the bath state impurity: where ρ is the density operator of the bath. For a pure state this measure is known to reach zero. In our simulation, we selected a directionality parameter of κ = 5 (Eq. 8) as an example case. Figure 2d shows that for this value of κ the impurity reaches a steady-state value of ≈ 10 −4 , which represents a negligible initialization error. After settling, the system is trapped in an oscillatory limit-cycle, resulting from the competition of This error can be made arbitrarily small by increasing κ, which we address in the following section.

Dependence of protocol performance on system parameters
In the absence of dephasing, there remains a fundamental trade-off between impurity and the convergence time to steady-state T c . This is because reducing the three-body interaction -increasing κ -ensures a smaller contribution from off-resonant processes, thus reducing the initialization error, but slows the dynamics to steadystate. For κ → ∞ one could expect bringing the impurity arbitrarily close to zero, but this comes with a prohibitively long convergence time, T c .
To visualize the achievable steady-state impurities, and related convergence times, we run a complete simulation of our pulsed protocol for a range of values of κ and for N = 8 and N = 128. Figure 3a confirms the clear trend of purity improvement with an increase in κ, where substantial degrees of purification are achieved past κ 1. Figure 3b shows the inverse N -dependence of the half-period of the collective three-body exchange, τ 0 ∼ 4πω c /N a 2 , and the number of stage iterations required to reach convergence proportional to N . The convergence time is a simple product of these two quantities and thus combines to indicating that the convergence time is dependent purely on the three-body interaction strength, and not the system size.
Verifying the above trends in the large-N limit becomes computationally prohibitive. As an efficient way to extend our results into this regime, we turn to a steadystate solver of a quantum master equation in which the coherent three-body exchange proceeds continuously and simultaneously with a central spin reset whose rate is Γ op = 2π/τ 0 . Figure 3c displays the steady state impurity, , for N up to 10, 000 and for the same range of values of κ as in Fig. 3a. This confirms effective purification for κ 1 for large N . The solid black curve in Fig. 3c is the prediction from a simple rate equation[37], i.e. N → ∞, which overlaps with the steady-state model in the large-N limit. Both models feature a ∼ κ −2 rolloff, as obtained analytically from a first-order expansion of the impurity in the rate equation model (see [37]). As shown in Fig. 3d, eigenvalue analysis of the rate equation allows us to calculate the convergence time, T c , agreeing with the size-independent T c ∼ 16πω c /a 2 behavior observed in the pulsed model (Fig. 3b).
Comparing the pulsed (Fig. 3a) and continuous (Fig. 3c) protocol performance, we note that maintaining the temporal separation between the central spin reset and activation of the three-body interaction with the ensemble allows reaching lower steady-state impurities. This property, illustrated by comparing experiments done in the continuous [21] and pulsed [15] regimes, can be straight-forwardly explained: continuously measuring the central spin reduces its ability to sense the ensemble noise, and in turn limits the achievable purity.

Resilience to system imperfections
Real physical systems deviate from the idealized behavior considered so far, as they are affected by central spin and ensemble inhomogeneous dephasing, as can occur if the ensemble spins have a spread of Larmor precession frequencies √ ∆ 2 ω i > 0 [15,16,37]. Each stage of our purification protocol relies on the exchange dynamics between the central spin and the ensemble. Working in this strong coupling regime renders the protocol equally sensitive to both the central spin and ensemble dephasing mechanisms [37], proceeding at rates Γ c and Γ b , respectively. We thus use the total dephasing rate Γ d = Γ c +Γ b as the relevant parameter and explore the robustness of our protocol against such system imperfection. Using the rate equation approach [37], we calculate the steady-state impurity, , as a function of the dephasing rate normalized by the three-body interaction rate, Γ d τ 0 /(2π), as shown in Fig. 4a (black curve). A significant alteration of the protocol performance is observed when the three-body interaction time τ 0 becomes longer than the typical dephasing time, i.e. Γ d τ 0 /(2π) 1. Nonetheless, for large values of the directionality parameter κ (here, κ = 10), the engineered steady state impurity remains lower than unity by orders of magnitude. In this low-impurity regime (and for N → ∞), this is quantitatively captured by a simple expression [37]: The rate of the activated three-body exchange is also affected by dephasing, which results in an increased protocol convergence time, T c (blue line in Fig. 4a). As expected from a time-energy uncertainty principle, dephasing also comes with a broadening of the three-body resonance Ω = ∆ω, as seen in Fig. 4b.

Benchmarking candidate material platforms
Throughout this work, we have identified the two platform-agnostic control parameters determining the protocol performance: the directionality parameter, κ, and the normalized total dephasing rate, Γ d τ 0 /(2π). We now evaluate these parameters for candidate physical systems and quantify the corresponding bounds on steadystate bath impurities, , as well as the protocol convergence times, T c , as shown in Fig. 5.
We restrict our analysis to the physical platforms that naturally realize dense central spin systems (see the Hamiltonian of Eq. 4), and in which the rudimentary protocol ingredients, like the control and reset of the central spin state, have been demonstrated previously. Our nonexhaustive selection of the candidate systems includes Gate-Defined Quantum Dots (QDs), Lattice matched GaAs-AlGaAs QDs, Stranski-Krastanow InGaAs QDs, and Rare Earth Ions (REI). The system-specific parameters used in the calculations are tabulated in the supplementary materials [37].
For QD systems the two spin species that break ensemble symmetry are simply two nuclear-spin species with different gyromagnetic ratios: gallium and arsenic. Calculations for QD systems have taken into account the spin-orbit magnetic-field B dependence of the centralspin lifetime, ∝ B 3 [44] and ∝ B 5 [45] for the Gate Defined QDs and the epitaxial QDs (both InGaAs and GaAs-AlGaAs), respectively. Beyond magnetic fields of 10 T, this spin-orbit effect becomes the most performance- limiting factor and caps the impurity-convergence time trade-off. Spectacular purification to < 10 −4 can be achieved for GaAs and Gate Defined QDs, due to the high values of the directionality parameter, κ, and the low total dephasing rates, Γ d τ 0 /(2π), within the typical range of externally applied magnetic fields, B (see Fig. 5a). Gate-Defined QDs feature ∼100-fold larger directionality parameters, κ, than the GaAs QDs at the same externally applied magnetic field, for the same (unnormalized) total dephasing rates, Γ d [14,24]. This generally leads to higher degrees of purification, but significantly longer convergence times for Gate Defined QDs (see Fig. 5b).
InGaAs QDs feature a larger effective nuclear dephasing rate, Γ d , due to nuclear spectral inhomogeneity arising from strain-induced quadrupolar broadening [25]. In their case, three-body interactions could still lead to a weak purification, as consistent with the observation of enhanced spin-wave modes at low ensemble polarization [22].
We now turn our attention to REI systems -specifically to 171 Yb 3+ : YVO 4 , used recently to demonstrate quantum-state transfer between a 171 Yb electronic central spin and the second shell of the nearest 51 V nuclei [18]. This system belongs to the regime of small but dense central-spin systems and offers little tunability with an external magnetic field as all nuclear spins are of the same species. However, the nuclear-spin shells surrounding the central spin can be distinguished via quadrupolar shifts, allowing a set of values for the effective ∆ω [37]. The required central-spin mediated three-body interaction arises as a second-order magnetic dipole-dipole coupling between two of the first, second, and higher shells of 51 V nuclei interfaced with the 171 Yb. The measured nuclear coherence time is three times shorter than the interaction time τ 0 , but this still allows for generating high-purity anti-polarized states (see Fig. 5).

Preparation of a many-body singlet
The singlet state |I = 0 is a superposition of all |I z 1 = −M + n, I z 2 = M − n eigenstates with a ∝ (−1) n phase on each state, as shown in Eq. 3. Having established a pure anti-polarized state of the two subensembles (i.e. n = 0) following the first two stages of the protocol, the third and final stage will prepare the singlet state by weaving an alternating phase into the Jaynes-Cummings ladder, as shown in Fig. 6a. We stress that the many-body singlet state is not an eigenstate of our system Hamiltonian; however, it refocuses every ∆t = 2π/∆ω following the protocol termination [37].
We consider only the ideal execution of this final stage involving unitary gates, free of any dephasing. We work in a frame co-rotating with ΩS x +∆ω(I z 1 −I z 2 )/2. For the sake of clarity, we outline the ideal protocol steps assuming that the effective Jaynes-Cummings ladder has 2M + 1 = 2 K rungs, corresponding to a total number of ensemble spins N ∼ (2 K − 1) 2 /2; the generalization to an arbitrary N is straightforward. In the first instance, we take a simplified scenario in which the three-body-interaction is independent of the |I z 1 = −M + n, I z 2 = M − n -state; we will later take the dependence on n into account.
We apply the following sequence of unitary operations, starting from the state (i) a central spin (π/2) z -gate, giving state as illustrated in the second level diagram from the left in Fig. 6a; (ii) a three-body-interaction π-gate, U (τ ) = U π , where τ is the interaction time, and which up to a global phase leaves the system in the state as shown in the third level diagram from the left in Fig. 6a. Together, this pair of gates (π/2) z and U π injects entanglement into the spin ensemble via the central spin. We combine them into a composite gate S 1 , as shown in Fig. 6a. Application of this gate doubles the overlap with the singlet state, | ψ 1 | (|↓ x |I = 0 ) | 2 , from the initial state's, | ψ 0 | (|↓ x |I = 0 ) | 2 . To increase this overlap further, we apply an extended composite gate, S 2 , which contains two steps: (i) entanglement injection using S 1 (the fourth and fifth level diagrams from the left in Fig. 6a) and (ii) phase redistribution onto the |I z 1 = −M + n, I z 2 = M − n -states with n = 0, 1, 2, 3 using a central spin (π) z -gate and a three-body-interaction π-gate (the two right-most level diagrams in Fig. 6a). We generalize this sequence to a composite gate S K , for which the phase redistribution is applied 2 K−1 − 1 times (gray box in Fig. 6b). Each application of a phase redistribution sub-sequence brings the next highest rung on the ladder into the ensemble superposition state. As a result, the overlap of system state |ψ j with the many-body singlet state | ψ j | (|↓ x |I = 0 ) | 2 doubles at each step of a sequence of composite gates S j , for j = 1, 2, .., K.
The modular structure of this algorithm is advantageous in terms of minimizing the impact of the central spin dephasing on the ensemble state preparation. It is sufficient for the central spin to stay coherent over a given S j gate duration, after which it can be safely reinitialized. We note that the phase redistribution is the most operation-costly sub-sequence, as it involves 2 j − 2 gates within each S j gate. Nevertheless, the protocol complexity remains linear with the number of states along the ladder 2 K as the total number of gates in a complete sequence is 2(2 K − 1).
We demonstrate the performance of this algorithm in Fig. 6c, using the trace distance from the engineered state to the singlet state as a function of the total number of gates (gray squares). For K = 4 (N = 112), we can visualize the steady progression from |ψ 0 towards the singlet state |↓ x |I = 0 as the algorithm steps through the sequence of composite gates S 1 to S 4 , culminating in an exact preparation where the final trace distance reaches zero. The classical limit I 2 cl ∼ √ N is overcome as soon as the sequence begins to inject entanglement into the ensemble, and reaches the quantum limit of I 2 qu = 0 at its termination. Interestingly, this happens at the expense of raising the uncertainty in I z 1 and I z 2 towards their thermal values ∼ √ N , akin to squeezing. We now turn to the more realistic description of the spin ensemble for which the threebody-interaction π-gate time depends on the state  I2). The illustrated sequence of six gates prepares a manybody singlet for K = 2. The inset displays a relative phase color-coding applied throughout the panel. The global phase is factored out, leaving the lowest energy state with a zero reference phase. b, Concatenation of the composite Sj gates turns the anti-polarized state into a many-body singlet. The gray box contains a quantum circuit corresponding to the SK gate, for an arbitrary integer K. c, Trace distance from a singlet state as a function of the number of singlet-preparing gates used within the simplified (gray squares) and the real (circles) systems with K = 4 (i.e. N = 112), following the exact and variational-searched protocols, respectively. The red-to-blue gradient indicates that the structure of the optimal protocol (discussed in the Ref. [37]) varies with the number of used gates.
ior complicates the algorithm implementation, the core structure used to inject entanglement and redistribute phase remains in place. With the right choice of the three-body interaction times τ and the central spin zgate phases, φ, it is possible to reach the singlet state. To do so we treat the sequence parameters {τ i , φ i } variationally for each gate to minimize the (Hilbert-Schmidt) trace distance to the singlet state from Formally, this involves an optimization over a ∼ 2 Kdimensional space of parameters, for which we employ a gradient-descent algorithm [37,46]. The resulting trace distances for the sequences of increasing length are presented in Fig. 6c as the colored circles (red-to-blue gradient is consistent with the color-coding used in ref. [37]). Strikingly, our algorithm arrives within a trace distance of a few % from the singlet state. A process optimization in larger systems, aided by machine learning algorithms, could prove to be a viable route for creating arbitrary state superpositions [47].

CONCLUSIONS AND OUTLOOKS
In this work, we have proposed a protocol that can initialize a spin ensemble into a pure anti-polarized state and steer it towards a many-body entangled singlet state -exclusively by using a single central-spin qubit. To do so, we have made use of the three-body interaction naturally present in dense spin ensembles and shown that it can be harnessed by breaking the ensemble into two spin species. We have suggested several platforms where this algorithm would be realizable, and where significant purity can be achieved even in the presence of dephasing. We note that in these systems, breaking the spin ensemble can take multiple forms including, but not limited to, nuclear-spin species with different gyromagnetic ratios [25] or high-spin species which split into two effective qubit ensembles under the influence of electric-field gradients (e.g. from strain) [48].
From the perspective of an electron spin qubit hosted in a material with non-zero nuclear spins, a singlet state of its surrounding spin ensemble would dramatically boost its coherence -both homogeneous and inhomogeneous noise sources [25] would be quenched altogether. From the perspective of leveraging this spin ensemble as a quantum memory resource [34], initialization to a pure collective state is sufficient to run an algorithm with unit fidelity, and the availability of two anti-polarized species could even be extended to a two-mode register. The state-engineering recipes that we have established could be extended to more elaborate computational [49] and error-correcting [20] algorithms. Fundamentally, tracking a many-body state in the presence of tuneable interactions can reveal the entanglement dynamics in and out of the central-spin system, opening an experimental window onto quantum information scrambling and area laws for entanglement entropy. example, measuring the transition rate asymmetry 2 of the Ω = ω 1 activated processes: (2) and the transition rate asymmetry of the Ω = ω 2 activated processes: constrains both I 1 and I 2 , which can be used to conditionally discard the post-protocol state, or to adjust the gate parameters in the third stage of the protocol.

II. FULL SIMULATION OF QUANTUM DYNAMICS
To simulate the quantum dynamics of the composite system in a reduced Hilbert space (see the manuscript) we use a master equation and steady state solvers from the QuTiP 3,4 library in Python.
The exchange between the central spin and the ensemble is modelled by the following Hamiltonian, written in a frame rotating with a δ-detuned laser drive, following a Rotating Wave Approximation: The individual terms are defined in the manuscript. Within the first stage of the protocol setting Ω(t) = 1 2 (ω 1 + ω 2 ) activates the non-collinear (∝ S z I x i ) interaction. Within the second and third stages of the protocol, Ω(t) = ∆ω activates the three-body interaction (∝ S z (I + 1 I − 2 + I − 1 I + 2 )). At all times δ = 0. The instantaneous central spin resets within the second stage of our protocol were effectuated by the following operation on the composite density operator: where Tr c stands for the partial tracing with respect to the central spin degrees of freedom. The steady state solver was applied to a continuous approximation of the protocol in which the central spin was reset continuously at an optical pumping rate Γ op = 2π/τ 0 , where τ 0 corresponds to the time between subsequent instantaneous central spin resets in the exact protocol. The collapse operator used in modelling that process was Γ op |↓ x ↑ x |. We verify the protocol robustness in a I 1 = I 2 manifold by performing an identical simulation to that from the manuscript section 'Ideal system dynamics'. Both κ and τ 0 are fixed to the same values, and the only difference lies in the choice of I 1 = 32/2 − 1 and I 2 = 32/2 + 1. The results are illustrated in the Fig.2, and organized in a oneto-one correspondence with the Fig. 2 of the manuscript. As displayed in the Fig. 2a, the S x population (that is, the only direct observable during the protocol) saturates to −1/2, like in the Fig. 2a of the manuscript. The polarizations I z 1 and I z 2 are driven towards the maximal possible opposite values, that is ± min(I 1 , I 2 ) (see Fig.  2c), as anticipated. Dynamics feature equal amount of purification after settling to the limit cycle (see Fig. 2d), and the only significant difference comes in the degree of I 2 reduction. Still, the protocol reaches the anticipated classical limit I 2 cl = (I 1 + I 2 )(|I 1 − I 2 | + 1).

B. Equal impact of the ensemble and the central spin dephasing on the protocol's performance
We studied the effects of dephasing of the central spin and the ensemble on the protocol's performance. The former process was modelled using a √ Γ c S x collapse operator, and the latter made use of two distinct collapse operators: Γ b /2I z 1 and Γ b /2I z 2 . Fig. 3 illustrates the steady state impurity calculated in the continuous protocol approximation for a range of dephasing rates Γ c and Γ b . It is evident that both processes have quantitatively identical impact on the state preparation, which motivates introducing a total dephasing rate, Γ d = Γ b + Γ c , as a figure of merit in quantifying the resilience to system's imperfections.

III. RATE EQUATION MODEL
A. Evolution of populations and the scattering rates A good quantitative understanding of the protocol dynamics can be gained from a simple rate equation model. Within this model, the dynamics is again restricted to the I z = 0 ladder of states, as a result of perfect total polarization locking. For convenience, we label the n th ladder state, |I z 1 = −(M − n), I z 2 = +(M − n) , with a principal quantum number: |n . The model assumes that at coarse-grained timescales, the coherences between different states across the ladder vanish, or in other words: The evolution of the population of the |n -state, p n , is then captured by the following equation: where r ± n stands for a rate of the population flow from |n -state to |n ± 1 -state. The ∝ p n term in the equation describes the 'out-flow' processes proceeding at a total rate r − n + r + n , as illustrated in the Fig. 4a. The first of the two processes involves a coherent exchange |↓ x |n → |↑ x |n − 1 driven by a three-body interaction, and detuned by ∆ − = Ω − ∆ω, followed by the central spin reset after time τ 0 -here approximated as proceeding concomitantly with the exchange, at the optical pumping rate Γ op = 2π τ0 . The second process consists of the same two steps, except that the dynamics proceeds between |n and |n + 1 states, for which the three-body interaction detuning is ∆ + = Ω + ∆ω. The rates of the outflow processes, r ± n , are then proportional to the populations of the excited states |↑ x |n ± 1 , where the constants of proportionality are given by the central spin reset rate, Γ op . This yields: We approximate the expectation value from Eq. 8 by a steady-state excited state population of a two level system constituted by |↓ x , n and |↑ x , n ± 1 states; it follows that: where Γ = 1 2 Γ op +Γ d incorporates the sum of phenomenological ensemble and central spin dephasing rates, Γ d . The effective drive strength, α ± n , corresponds to the three body interaction rate, given by: Its dependence on I and n is a result of collective enhancement, dependent on the total angular momenta of of sub-ensembles (I 1 = I 2 = I ∼ N/2). The second term in the equation 7 captures the effect of the 'in-flow' processes. The |n -state can be populated by ∆ + -detuned process originating from |n − 1 -state, or a ∆ − -detuned process originating from |n + 1 -state, as shown in the Fig. 4b. The relevant scattering rates are found using Eq. 9.

B. Steady state solution
The system reaches steady state whenṗ n = 0 for all n. This generates the following recursive relation between the |n -state populations: This fully determines steady state populations, up to a normalization factor, which can be constrained using n p n = 1. The computational complexity of solving for steady state populations is O(I); this is to be compared with O(I 6 ) complexity of solving for steady state of quantum master equation with dissipation.
C. Scaling of impurity, , with directionality parameter, κ, and dephasing Γ d τ0/(2π) The preparation performance can be characterized by the impurity of the reduced density operator, ρ b = Tr c ρ. The rate equation model approximates this impurity by: In case of the perfect state preparation we have p 0 = 1, and = 0. The preparation performance starts to drop when the |1 -state acquires a finite population, leading to an initial impurity increase (calculated at at Ω = ∆ω): where we introduced α ≡ α + 0 = α − 1 . Since α 2 ∝ N , and Γ op Γ ∝ N 2 (as Γ op = N a 2 /(2ω c ), for the optimal value of τ 0 ), in the thermodynamic limit, this expression simplifies further: where κ: as defined in the manuscript. The requirement for the cooling to proceed optimally is therefore κ 1 and Γ d /Γ op 1.

D. Convergence Time, Tc
The rate equation formalism allows to estimate the convergence time of the protocol, T c , from spectral analysis of the matrix Λ, that generates the rate equation as: where p is a vector of |n -state populations. We notice that due to the contractive nature of dynamics, the real parts of the Λ-matrix eigenvalues are all non-positive, and can be arranged according to: We then identify T c = 2π/| Re(λ 1 )|, as | Re(λ 1 )| is the slowest convergence rate in the model.  11 . The three-body interaction in this system is a central spin mediated second-order dipole-dipole, and it couples the nuclei in the second shell (subscript i) with those in the outer (subscript j) shells with an effective strength ∼ aiaj/ωc. The ∆ω corresponds to the difference in quadrupolar frequencies of the nuclei in the second shell (i.e. a frozen core) and the outer shells.

V. OPTIMIZING PREPARATION OF MANY-BODY SINGLET
A. Central spin gate We introduce the following notation for a single qubit gate: where σ x , σ y , and σ z are 2 × 2 Pauli matrices.

B. Three-body interaction gate
Upon matching the Ω = ∆ω resonance condition, the evolution of the system in the frame rotating with 1 2 ∆ω(I z 1 − I 2 2 ) + ΩS x is dictated by the following Hamiltonian: where I ± = I ± 1 I ∓ 2 are the non-linear ladder operators for the effective Jaynes-Cummings ladder of |n -states, whose action is captured by: and I − = (I + ) † . Evolution of the system in that frame over time τ , generated by H exc , realises the following gate: cos a 2 e n τ 4ω c − 1 |n n| i sin a 2 e n τ 4ω c |n + 1 n| .
In particular, it is readily seen that the exact π-gate time, τ π , is dependent on n since the enhancement factors e n vary across the ladder of states. In contrast, for the simplified system with e n = const, the exchange interaction π-gate would be: FIG. 5. Optimizing gate times and phases for singlet preparation. a, Hilbert-Schmidt distance (i.e. our cost function, C) during simplified (squares) and optimized (circles) protocols of varied length. The colorcoding is consistent with Fig.6c of the manuscript, which only concerns the protocols' endpoints. b, Evolution of the trace distance in considered protocols. c, Evolution of the expectation value of total angular momentum squared in considered protocols. d, Optimal Z-gate phases in protocols of varied length. The beige circles correspond to initial condition for the gradient descent optimization. e, Optimal exchange gate times in protocols of varied length. The beige circles correspond to initial condition for the gradient descent optimization. 'cpl' in the label stands for coupling strength, a 2 /(4ωc).

C. Optimizing the gates via gradient descent
In order to tailor the sequence parameters in a real system, we minimize the appropriately constructed cost function, which takes the following quantum state: as an input (a method originally devised in Ref. 12 ). The arrow over the product sign represents direction of stacking the consecutive terms with j = 0, 1, 2, .. in the product. The minimization is done with respect to the variational parameters {τ j , φ j } j=0,1,2,.. , which correspond directly to the three-body interaction activation times (τ j ), and phases of central spin gates (φ j ). To reach the minimum of the cost function we apply the RMSprop gradient descent algorithm in the space of {τ j , φ j } j=0,1,2,.. .

Choice of the cost function
We choose a Hilbert-Schmidt distance between the singlet state, χ = |I = 0 I = 0|, and the reduced density operator of the ensemble, ρ b = Tr c (|ψ ψ|), as our cost function: The main advantage of working with this cost function is an ease of calculating its gradient analytically, which speeds up the optimization procedure. To see it, we first act with the differential operator, on the |ψ -state from the Eq. 23, and find: for variational parameters v j ∈ {τ j , φ j } and operators X j ∈ {H exc , σ z /2}. We then notice that extending this result to the cost function from Eq. 24 amounts to a simple application of a chain rule.

Hyperparameters and convergence of an RMSprop algorithm
The update rule for a given variational parameter, v, in our implementation of the RMSprop algorithm is: where ∂ v C is the partial derivative of the cost function, C, with respect to the variational parameter, v. The hyperparameters used throughout the optimization routines were ζ = 0.015, ξ = 10 −8 , and β = 0.85. For all the sequences of lengths shorter than the maximal considered, the convergence was reached in less than 1000 optimization epochs. For the longest considered sequence, optimization was terminated after 7000 epochs, however, the cost function could likely be decreased further. The difficulty of this task comes most likely from a presence of a barren plateau in a cost function landscape -i.e. region around the minimum, where the gradients of C vanish exponentially fast 13 .
The evolutions of the cost function, C, the trace distance, and the expectation value of total angular momentum squared during the optimal protocols of varied length, are plotted in the Fig. 5a, Fig. 5b, Fig. 5c, respectively (circles). Their simplified system's counterparts are plotted alongside, for reference (squares).

Initial condition for gradient descent
To speed up the convergence of the gradient descent optimization, we start the procedure with a physicallymotivated initial guess for variational parameters.
The initial z-gate phases are chosen to be identical to those in the optimal protocol for the idealized system (see the main text, and the beige circles in the Fig. 5d). The initial activation times for the three body interaction are chosen as: where e j stands for the enhancement factor from the Eq. 20 -see the beige circles in the Fig. 5e. The intuition behind this choice, is that it realizes the ideal amplitude transfer for each consecutive |n -state brought into the superposition.

D. Singlet auto-refocusing in a non-rotating frame
Singlet state prepared in the frame rotating with 1 2 ∆ω(I z 1 − I z 2 ) + ΩS x , will coincide with a singlet state in a non-rotating frame at times t = 2πk/∆ω for k = 0, 1, 2, 3..; indeed in a non-rotating frame the evolution of the prepared state (up to a global phase) is given by: