Entangling nuclear spins in distant quantum dots via an electron bus

We propose a protocol for the deterministic generation of entanglement between two ensembles of nuclear spins surrounding two distant quantum dots. The protocol relies on the injection of electrons with definite polarization in each quantum dot and the coherent transfer of electrons from one quantum dot to the other. Computing the exact dynamics for small systems, and using an effective master equation and approximate non-linear equations of motion for larger systems, we are able to confirm that our protocol indeed produces entanglement for both homogeneous and inhomogeneous systems. Last, we analyze the feasibility of our protocol in several current experimental platforms.

In this work we follow this research line, and aim to design a protocol that produces entanglement between two distant nuclear ensembles deterministically.This has been pursued in Refs.[25,26] for the two nuclear ensembles surrounding a double quantum dot.For an analogous operation at longer distances, as explored in Refs.[27,28] for entanglement generation between distant electrons, one may use a coherent electron bus, such as a longer array of quantum dots [29][30][31][32], a chiral quantum-Halleffect channel [27], or surface acoustic waves [33][34][35].The idea behind most of these proposals is that the moving qubits (a current of electrons) act as a joint bath, whose overall effect is to induce the decay of the system towards a target steady state.By engineering the bath we can modify the properties of this steady state, in particular, we can ensure that the two parties (the two electrons or two ensembles of nuclear spins) are entangled.
Our work comprises the first proposal for entanglement generation between distant ensembles of nuclear spins connected via an electron bus.Two different regimes are of potential interest to apply our scheme: On the one hand, the large nuclear ensembles that surround QDs fabricated with GaAs, InGaAs, and 29 Si-rich Si, with a number of spins ranging from 100 to 10 6 [36]; on the other, nuclear ensembles with just a few nuclei (< 10) present, for example, in purified Si [24,37].We also discuss how to detect the entanglement experimentally using different entanglement witnesses suitable for each case.

II. DESIGN OF THE PROTOCOL A. System
Our protocol is designed for a device such as the one shown schematically in Fig. 1: Two QDs, henceforth labelled as QD1 and QD2, are connected to electron reservoirs so that electrons with definite polarization can be injected into each QD, and they can be released from either QD at will.Furthermore, the two QDs are connected with each other via a coherent electron bus, which can transport electrons from QD1 to QD2.
If QDj is occupied by an electron, the nuclei in a small volume around it will interact with it via the hyperfine interaction.We consider the Fermi contact interaction between an s-type conduction band electron and the nuclear spins, which is given by the Hamiltonian [1] H j = g j A j • S j = g j A z j S z j + Here, S j denotes the electronic spin, while A j = i g ij I ij is a collective nuclear spin operator; I ij is the spin of the ith nucleus in QDj.The ladder operators are defined as usual S ± j = S x j ± iS y j (analogously for A ± j ).We consider spin-I nuclei, and a total of N j nuclei in each ensemble.The dimensionless coupling constants g ij may FIG. 1. Schematics of the system under consideration and the manipulations performed in each step of the protocol.
differ depending on the electronic wavefunction at the positions of the nuclei and the nuclear species, they are normalized such that i g 2 ij = 1/(2I) [38].The constants g j = A H / i g ij set the overall strength of the hyperfine interaction in each QD; A H is a constant that depends on the material and the density (concentration) of nuclei around the QD [1,39].
Within this model there are implicit several approximations: We assume that the electron injection/retrieval process is fast on the scale of the hyperfine interaction, so that H j is switched on/off abruptly whenever an electron is added to/removed from QDj, and that it is spinindependent [34].We also neglect other effects, such as non-contact terms of the hyperfine interaction and nuclear dipole-dipole interactions, which usually are much weaker compared with the contact interaction [40].We also neglect spin-orbit coupling, which is weak in GaAs, and all electron-spin decoherence processes.These simplifications are discussed and justified in Section V.

B. Protocol
The protocol consists of a sequence of interactions that generate a dissipative dynamics whose principal (Lindbladian) terms stabilize a spin-squeezed state of the two ensembles [41,42].The interactions are designed such that the main unwanted dephasing terms cancel out.Specifically, the protocol consists in the repetition of a cycle comprised of 2 + 2k steps (k ≥ 1).At the beginning of each step, a new electron with definite polarization is injected into QD1.This electron interacts with the nuclei in QD1 for a certain period of time τ 1 and then it is transported to QD2, where it interacts with the nuclei there for a period τ 2 , see Fig. 1.At the end of each step the electron in QD2 is released from the device, leaving it ready to start a new step.The dwell times τ 1 and TABLE I. Electronic polarizations and dwell times for the steps that make up a single cycle of the protocol.The third and fourth rows are steps repeated a total of k times each.

step(s)
|sin g1τ1 g2τ2 τ 2 , and the polarization of the injected electrons at each step are required to fulfill a certain pattern: in the first step (injected spin "up") g 1 τ 1 = µ and g 2 τ 2 = ν, with τ 1 > τ 2 (|µ| > |ν|), while in the second step (injected spin "down") g 1 τ 1 = ν and g 2 τ 2 = µ.In the remaining 2k steps each nuclear ensemble evolves independently under the influence of k consecutive electrons, all with the same polarization ("down" for QD1 and "up" for QD2), that spend each a short time τ j ∝ k −1 trapped in QDj.All this information is summarized in Table I.The rationale behind this particular choice of the dwell times and injected electron states is the following: For short dwell times, we can neglect higher-order terms in the hyperfine interaction and each electron flips at most once when interacting with the nuclear ensembles.In the first two steps of the protocol, each electron is transported coherently between QD1 and QD2, such that interference between the "flipped in QD1" and "flipped in QD2" paths occurs.This uncertainty regarding where the electrons have flipped will be mapped onto an entangled nuclear state after tracing out (discarding) the electrons.Simultaneously, these two steps produce unwanted dynamics due to the zz hyperfine terms.terms.On the nuclei, these latter terms act like an (in general inhomogeneous) effective magnetic field (the Knight field) that can be cancelled with additional local terms, which are produced in the remaining steps of each cycle.
This specific pattern for the dwell times may be hard to fulfil in actual experiments if the coupling constants g j are not well known.However, the protocol also works if the pattern is satisfied only approximately.We are also assuming a perfectly coherent electron transfer from QD1 to QD2, but this again is not strictly necessary.For the sake of clarity here we will focus on the "ideal" implementation of the protocol, as given in Table I, while the results including imperfect dwell times and electron transfer are shown in the Supplementary Material.
The induced nuclear dynamics is Markovian by construction, since discarding the electron after each step results in a "memoryless" environment for the nuclei.According to the process described, the nuclei's reduced density matrix at the end of the nth step ρ n , only depends on the nuclei's reduced density matrix at the end of the previous step, In this equation χ n−1 = ρ n−1 ⊗|s in s in | is the density ma-trix of the entire system at the beginning of the step, and |s in is the state of the injected electron.The Liouvillians that appear in it correspond to the coherent evolution given by the hyperfine interaction If the dwell times are small, τ j < H j −1 (τ j < 2I −1 A −1 H ), we can approximate the dynamics by truncating the Taylor expansion of the exponentials up to second order.Then, the change in the density matrix after a single cycle ∆ρ is given by (see Supp.Material) where

C. Entanglement criteria
To detect entanglement between the two nuclear ensembles in an experiment one can employ a spin-squeezing criterion [25,41] based on the following principle [43]: Given two observables X j and Y j of each ensemble j = 1, 2 satisfying the commutation relation [X i , Y j ] = δ ij Z j , and linear combinations u = αX 1 + βX 2 and v = αY 1 − βY 2 , then for any separable state the inequality holds.Choosing α = β = 1, X j = J x j and Y j = (−1) j+1 J y j , where J ν j is the ν-component (ν = x, y, z) of the total nuclear angular momentum J j ≡ i I ij , we conclude that if ∆ EPR < 1, where is the so-called EPR uncertainty, then there is entanglement between the two nuclear ensembles.Notice that this is a sufficient condition for the state to be entangled, but it is not necessary, i.e., there are entangled states that cannot be detected by this criterion.Experimentally, this criterion requires the measurement of the nuclear polarizations of each ensemble.This may be done using NMR techniques [44], although the small relative size of even the largest ensembles in typical QDs makes it challenging to detect the nuclear signal of only the relevant nuclei-those that couple to the quantum dot electrons.On the other hand, for small ensembles (N j < 10) one may use other techniques [45,46] that allow single nucleus addressing Instead, the most natural way to measure observables of the nuclear spin ensemble in a QD is via the interaction with the electron spin, i.e., by measuring the Overhauser field [3].This allows to magnify the small nuclear signal and ensures that only spins actually coupled to the electron contribute.This technique, however, does not give access to the total nuclear spin operators if the couplings are inhomogeneous.But the principle described above, Eq. ( 4), is more general and also applies to the collective operators A j ; the variances of non-local operators involving these then have to be compared to expectation values of such that ∆ < 1 implies that the two ensembles are entangled.To avoid confusion, we will refer to the expectation value J z j /(N j I) simply as the polarization of the jth ensemble, while we will refer to 2 A (z,2) j as the generalized polarization (the range of both quantities is [-1,1]).Note that for an homogeneous system with identical ensembles (N 1 = N 2 ) both ∆ EPR and ∆ are equal.
One way to measure A z j in a gate-defined single QD is to measure electron spin resonance in an external magnetic field B z , e.g., by determining the resonance frequency at which spin blockade of transport through the QD is lifted [47].By rotating the nuclear spins using a suitable NMR pulse the same method can be used to measure A x,y j (and their variances).However, there seems to be no direct access to A (z,2) j .One possibility is to exploit the Jaynes-Cummings-like coupling, H ⊥ j , of the electron spin to the nuclear ensemble.Neglecting H z j for the moment, the probability for an electron spin prepared in the state "down" to flip after a short time τ < 1/g j is, to leading order in τ , approximately given by τ 2 g 2 j A + j A − j /4, and similarly the one for the flip of an initially spin-"up" electron is τ 2 g 2 j A − j A + j /4.The difference of the two gives The presence of H z j in the hyperfine Hamiltonian leads to corrections of third and higher order, −iτ 3 g 3 [A z j , {A + j , A − j }] /16, which could be made even smaller if we add an external magnetic field (pointing in the z-axis) so as to compensate the Overhauser field A z j +gµ B B ext = 0.In that case, the correction would be given by −iτ 3 g 3 [δA z j , {A + j , A − j }] /16 (δA z j ≡ A z j − A z j ), which is small for the highly polarized states we are considering.
Besides the criteria discussed here, there are many other ones one could employ to detect bipartite entanglement experimentally, each of which comes with its own challenges in how to implement the measurements and to which states the criterion is most sensitive.For example, for two single-spin-1/2 ensembles one could detect entanglement measuring the singlet occupation [48].Another more general option is to employ criteria based on the covariance matrix [49,50].

III. HOMOGENEOUS CASE
To demonstrate that this protocol indeed produces entanglement, we begin analyzing the dynamics in the  11), for the subspace with maximal total angular momenta Jj = NjI (rj = 1).(b) Characterization of the steady state in the subspace with maximum total angular momenta as a function of the number of nuclei, for two identical (N1 = N2 ≡ N ), homogeneous ensembles of spin-1/2 nuclei.We compare the exact dynamics (circles) with the result given by the second-order master equation (squares), which does not depend on the individual values of µ and ν but only on their ratio.For all the curves plotted this ratio is the same, ν/µ = 0.8.The dashed horizontal line marks the value obtained by the Holstein-Primakoff approximation in the limit Jj → ∞, ∆ HP EPR .(c) Steady-state ∆EPR for an initial mixed state: blue circles and dots correspond to an initial equal-weight mixture of all states with J1, J2 ≥ Jmin, while orange squares and crosses correspond to an initial state with a Gaussian probability distribution, pj(J) ∝ exp[−(J − J) 2 /(2w 2 )], of the total angular momenta with mean J and variance w 2 = 3.The system is the same as in panel (b), with N = 40 spins.The ratio ν/µ = 0.8 is the same for all the points.The results have been obtained using the Holstein-Primakoff approximation and the 2nd-order master equation.
homogeneous case, i.e., when all nuclear spins belonging to the same ensemble couple with the same strength to the electronic spin (g ij = 1/ 2IN j , for all i = 1, . . ., N j ).In that case, the collective nuclear spin operators form a spin algebra A j ∝ J j , and the total angular momentum of each ensemble J j ∈ {N j /2, N j /2 − 1, . . ., (N j mod 2)/2} is conserved throughout evolution.Then, the density matrix can be split in blocks corresponding to different sectors (subspaces) with different values of (J 1 , J 2 ) that evolve independently, and the steady state has no correlations between these blocks.
The last two dissipators of Eq. (3) drive the spin ensembles towards an entangled steady state, while the first two dissipators locally depolarize the nuclear spin ensembles and are detrimental for entanglement generation.If evolution were dictated just by the second line of Eq. (3) (limit k → ∞), then the steady state in the subspaces with equal total angular momenta J 1 = J 2 = J would be a two-mode spin squeezed state [51], given by the pure (unnormalized) state written here in terms of Dicke states |J, M j [52], which are eigenstates of J 2 j and J z j with eigenvalues J(J + 1) and M respectively.In this state, for sufficiently different µ and ν the two nuclear ensembles are highly polarized in opposite directions, and the EPR uncertainty turns out to be independent of J, ∆ EPR = (µ − ν) 2 /|µ 2 − ν 2 |, hence ∆ EPR < 1 for µν > 0, so the steady state is entangled.For unequal total angular momenta J 1 = J 2 , the steady state corresponding to the dissipators in the second line of Eq. ( 3) is not pure, although still entangled [26].For non-maximal J 1 , J 2 the steady state is not unique, but the (spin permutation) quantum numbers associated with this non-uniqueness [52] are conserved by the homogeneous dynamics and do not affect the EPR uncertainty, so they can be ignored.
By increasing the number of steps per cycle (make k larger), we are able to bring the actual steady state of Eq. (3) closer to |Ψ TMS (see Supp.Material).There is, however, a physical limitation on how small the dwell times can be made while still retaining control of the process.For example, in Ref. [34], the injection/ejection pulses take ∼ 90 ps, and the time-jitter is quoted as 6.6 ps.Nonetheless, as we show in the following, the protocol also works for small values of k (even for k = 1).
As extracted from the analysis of the dynamics in the limit k → ∞, the two nuclear ensembles become polarized in opposite directions.This is also true when considering a finite value of k.In this case, we can obtain an analytical estimate of ∆ EPR , valid in the limit of large J j , bosonizing the excitations of each ensemble over the fully polarized states |J 1 , J 1 1 and |J 2 , −J 2 2 using the Holstein-Primakoff (HP) transformation [53], which maps spin operators J j to bosonic operators a j , a † j as Note that the role of the ladder operators is reversed in QD2.If a † j a j J j , we can approximate and similarly , where r j = J j /(N j I) is the ratio between the total spin and the maximum possible total spin for the jth ensemble.Rewriting Eq. ( 3) using this transformation we obtain a new master equation whose steady state can be found analytically (see Supp.Material).For k = 1, the EPR uncertainty in the steady state reads where κ ≡ (µν − ν 2 )/4 > 0. Considering the relation between the arithmetic and geometric means, we can readily see that ∆ HP EPR is smaller the more similar the total spins J 1 and J 2 are, and also the more similar the ratios r 1 and r 2 are.In Fig. 2(a) we plot this quantity for the subspace with maximal total angular momenta (J j = N j I), showing the values of N 1 /N 2 and ν/µ for which we expect to detect entanglement in the steady state.
In addition, this bosonization approach allows us to compute the dynamics, provided the conditions a † j a j J j are satisfied at all times.This is indeed the case for sufficiently small ν/µ and for the initial fully polarized state denote the states in which all nuclei in QDj are polarized "up" or "down" respectively).The resulting dynamics approaches the steady state polarization exponentially with rate κr j in each QD (see Supp.Material).If both ensembles are initially highly polarized (r 1 , r 2 1), we can take the associated time T hom stab ≡ 1/(2κ) as an estimate of the timescale in which the system reaches the steady state.
In Fig. 2(b) we compare the steady-state polarizations and ∆ EPR obtained by means of the exact Eq. ( 2) with the ones obtained using the second-order master equation, Eq. ( 3), for a system consisting of two identical, homogeneous ensembles (initially fully polarized).While the polarizations do not present visible differences, there is some discrepancy between both in terms of entanglement, which is reduced by making the interaction times smaller.As we can see, ∆ EPR decreases monotonically as we increase N (J j ), and the approximation ∆ HP EPR agrees well with the result of the master equation already for N 25.Thus, in order to obtain the maximum amount of entanglement it is preferable that initially the subspace with largest occupation be the one with maximum total angular momenta (J 1 , J 2 ) = (N 1 I, N 2 I).This requirement on the initial state could be enforced polarizing completely the two nuclear ensembles, a task for which many protocols have been devised [e.g.16, 54-57] and towards which significant experimental progress has been reported [5,[58][59][60].Some limitations to dynamical nuclear polarization have also been predicted [61].
We remark that the protocol also works for nonmaximal J j and even mixed initial states.For example, we may consider an initial state of the form ρ 0 = ρ QD1 ⊗ρ QD2 , with ρ QDj = J p j (J)ρ j (J).Here, ρ j (J) is some normalized density matrix that belongs to the subspace of states of the jth nuclear ensemble with total angular momentum J, and p j (J) is some discrete probability distribution.Since the steady state within each fixed angular momenta subspace is unique (up to variations of the permutation quantum numbers), the steady-state value of ∆ EPR only depends on the initial distributions p j (J).In Fig. 2(c) we show the resulting steady-state ∆ EPR for an initial Gaussian distribution and also for an equal-weight mixture of all states with J 1 , J 2 ≥ J min , demonstrating that maximal polarization of the nuclear ensembles is not strictly required for the protocol to work (see Supp.Material for further details).

IV. INHOMOGENEOUS CASE
As explained in the beginning of Section III, the homogeneous case is numerically simpler because of the conservation of total angular momentum.However, in the inhomogeneous case, the total angular momentum of each nuclear ensemble is not conserved and we have to use the full Hilbert space in order to compute the dynamics.This limits the direct application of Eqs. ( 2) or (3) to very small systems.Henceforth, we will consider spin-1/2 nuclei, I α ij = σ α ij /2 (α = x, y, z).In Fig. 3 we show the expected dynamics for a system with only a few spins in each QD, considering both homogeneous and inhomogeneous couplings.For inhomogeneous couplings the criterion based on ∆ fails to detect entanglement, yet the steady state is entangled as can be demonstrated by the PPT criterion [62][63][64].Elapsed cycles 3. Simulation of a small system with N1 = N2 = 2 nuclei in each ensemble computed using the exact dynamics, Eq. ( 2).The initial state is |ΨFP .The rest of parameters are µ = 0.1/ √ N1, ν = 0.8µ, k = 5.Comparison between the homogeneous case, gi1 = gj2 = 1/ √ 2 (blue), and a case with a small inhomogeneity, g11 0.714, g21 0.7, g12 0.725, g22 0.689 (orange).Even though ∆ > 1 for the inhomogeneous case, the steady state is entangled, as it is shown in the bottom panel, where we plot the lowest eigenvalue of partial transpose nuclear density matrix with respect to the nuclei in QD2, ρ T 2 .Note how the small difference between the couplings in the two cases considered only affects the dynamics at long times.
For larger ensembles, we concentrate instead on the quantities required to compute ∆, i.e., A (z,2) 1,2 and var(A x,y 1 + A x,y 2 ).They depend on the elements of the covariance matrix and lower-order terms like σ x ij , whose equations of motion (EoMs) can be derived from Eq. ( 3) as (see Supp.Material) As can be seen in these equations, the elements of the covariance matrix depend on zz-correlations, λ for (i, j) = (i , j ), and higher-order correlations.Due to the spin-operator commutation relations, and since the Lindblad operators depend on linear combinations of σ ± ij only, the only higher-order terms occurring in the EoM for γ ii jj and λ ii jj are of the form σ + i σ z j σ − k .In order to obtain a tractable system of equations, it is necessary to introduce some factorization assumption which effectively approximates these higher-order correlations in terms of lower-order ones.There are several approximation schemes one can use for this purpose [54].Among all the approximations we have investigated (see Supp.Material), the one which best approximates the dynamics is for (i, j) = (k, l) = (m, n) and (m, n) = (i, j).
In Fig. 4(a) we employ these EoMs to compute the dynamics of a large inhomogeneous system (larger than the one shown in Fig. 3, but still much smaller than realized in experiments with typical GaAs or Si QDs).The results show that for not too inhomogeneous couplings the steady state is entangled, and that this entanglement can be detected experimentally via ∆.Furthermore, for all cases considered ∆ < 1 after some initial cycles, albeit at later times ∆ takes on larger values and deviates away from this minimum, reaching ∆ > 1 for the more inhomogeneous cases.For moderately inhomogeneous systems there is a clear distinction between two timescales.On the one hand, there is the timescale in which the system reaches the minimum ∆, which roughly corresponds to the stabilization time of the equivalent homogeneous system, T ent ∼ (2κ) −1 .On the other hand, there is the timescale in which the system reaches the steady state, T stab , which is now much longer than in the homogeneous case.Also, estimating this timescale is harder in the inhomogeneous case, as the dynamics do not follow a simple exponential decay.The dynamics of the generalized polarization of each ensemble, A (z,2) j , is perhaps the simplest to analyze, as it depends just on the coupling constants within the given ensemble.Numerically, we can show that the stabilization time scales with the smallest difference of the coupling constants of the corresponding ensemble as −2 , see Fig. 4(b), where T stab has been computed as the time in which Besides computing the dynamics, we can solve directly for the steady state, that is, solve the system of equations ∆γ jj ii = 0 and ∆λ jj ii = 0. Note that for the initial state ρ 0 = |Ψ FP Ψ FP |, many of the covariance matrix elements and zz-correlations are equal: γ 11  ii = 1, γ 22 ii = 0, γ jj ii = 0, for (i, j) = (i , j ); λ jj ii = 1, and λ 12 ij = −1.We realize that for this initial condition, if two nuclei that belong to the same QD have the same coupling, say g ij = g i j , then γ jl ik = γ jl i k and λ jl ik = λ jl i k for all k ∈ QDl throughout the entire time evolution.Therefore, we can reduce the number of variables considering just one representative for each possible value of the couplings.This is what we will refer to as "shell model", the shells being the sets of nuclei within each QD with the same coupling constant.
Numerically, we find for this initial state that correlations between nuclear spins in different shells of the same QD vanish in the steady state.As a consequence, the steady-state polarization of each shell is independent of the others, and can be found analytically solving a simple quadratic equation (see Supp.Material).Furthermore, the steady-state polarization of each shell only depends on the values of µ, ν, and the number of nuclei in the shell, but not on the specific value of the coupling constants.On the other hand, the steady state value of the transverse spin variance depends continuously on the specific values of the coupling constants.
The condition that two nuclei have precisely the same coupling constant is somewhat exceptional in a real experiment, yet the steady state depends on this condition, as it determines the number of shells and the number of nuclei in each shell.In any case, small variations in the couplings only affect the long-term dynamics, making it possible to have an entangled state at intervening times for almost any system with sufficiently homogeneous couplings, as Figs. 3 and 4 suggest (see the Supp.Material for further calculations supporting this conclusion).

V. EXPERIMENTAL REALIZATIONS
The characteristic timescales of our protocol depend strongly on the materials used to make the device.Table II summarizes at a glance the typical energy-and timescales for different experimental setups.One of the most common material for QD fabrication is GaAs.All isotopes of both Ga and As have nuclear spin I = 3/2 [1].Since the dwell times are inversely proportional to A H , the large value A H ∼ 91.2 µeV (average value considering the different nuclear isotopes and their natural abundance) [60] in GaAs imply short dwell times in the order of a few to tens of picoseconds.
Another well-known material for QD fabrication is Si.Si has two isotopes, but only one of them is magnetic, 29 Si, with a natural abundance of around 5% and nuclear spin I = 1/2.Furthermore, the concentration of 29 Si in Si can be controlled by different methods of enrichment or purification.This allows great flexibility in the design of the device, with ensemble sizes ranging from a few nuclei to tens of thousands of them.The smaller hyperfine constant A H 4.3 µeV (for Si with a 100% concentration of 29 Si) [65] results in dwell times that can range from a few nanoseconds to a few microseconds depending on the concentration of spinful nuclei.
It is important to take into account other effects we have neglected up to now, which lead to nuclear spin decay and decoherence.The T 1 and T 2 times of nuclear spins are long, but finite, and protocols that are slow compared to these times cannot be expected to be well-described by our model (and thus, likely, are not able to produce entanglement).In this regard, the most relevant processes that would hinder the generation of entanglement are nuclear dipole-dipole interactions.In GaAs nearestneighbor dipole-dipole interactions are estimated to be ∼ 10 −5 µeV [1,55] corresponding to T dd ∼ 0.1 ms; in Si, where the distance between nuclei is larger and they have ), whereas the time in which we expect to produce entanglement is Tent ∼ (2κ) −1 cycles, which is of the order of ∼ 11IN cycles.In the last column we show the nuclear dephasing time due to nuclear dipole-dipole interactions.weaker magnetic moment, nearest-neighbor dipole-dipole interactions can be estimated as ∼ 4 × 10 −7 µeV (for 29 Si at concentrations larger than 12.5%), corresponding to T dd ∼ 1 ms.Furthermore, reducing the concentration of 29 Si would decrease also the strength of these unwanted interactions, leading to longer coherence times for the nuclei.In GaAs QDs, entanglement is generated after tens of microseconds, while one needs hundreds of microsends for Si devices.This result is almost independent of the concentration of 29 Si because T ent ∝ N/A H .According to these estimations, it seems that reaching the highly entangled transient state is indeed feasible in both GaAs and Si.

Material
Ferromagnetic leads or spin filters can be used for spin polarized electron injection [66][67][68].Concerning the coherent electron bus, great progress has been achieved towards this goal using surface acoustic waves in piezoelectric materials [33,34,69], and transport across long QD arrays [29][30][31][32].Although Si is not piezoelectric, the former bus could be implemented combining Si with other materials in a layered structure [70].One important imperfection that may arise in experimental implementations is the dephasing of the electron spin during transport from QD1 to QD2.If the transported spin would fully dephase, it would induce two independent jump operators µA + 1 and νA + 2 instead of the combination µA + 1 + νA + 2 (and similarly for νA − 1 + µA − 2 ).This could be modelled including a dephasing channel [71, p. 384] between the action of e τ1L1 and e τ2L2 in the first two steps of each cycle.The new master equation would be as Eq. ( 3), with the non-local dissipator terms D(µA + 1 + νA + 2 ) and D(νA − 1 +µA − 2 ) re-scaled by a factor (1−p) and additional terms: p D(µA + 1 ) + D(νA + 2 ) + D(νA − 1 ) + D(µA − 2 ) /4.Note that two of these new terms add up to the local polarizing terms that we considered already, whereas the other two constitute new local depolarizing terms.If the strength of these new terms is much smaller than that of the local polarizing terms, pµ 2  (µ − ν) 2 /k, then the polarizing terms dominate and we expect the nuclear dynamics to be very similar to the one in the limit of perfect coherent transport.One could estimate the phase-flip probability from the ensemble-averaged spin dephasing time of the electron spin during transport, which is averaged due to its motion (motional narrowing).The numbers reported in Ref. [69] show that surface-acousticwave transport in GaAs itself adds little dephasing, and an estimate of the dephasing time due to the motionally averaged Overhauser field of the nuclei in the channel suggests a negligible phase-flip probability p ∼ 1% for a channel of length 5 µm.Further details about the effect of electron-spin decoherence during transport are given in the Supplementary Material (sections A, C, and Fig. C1), where we show that even for a phase-flip probability p ∼ 5% it is possible to obtain an entangled steady state (for an homogeneous system) depending on the rest of system parameters and the choice of dwell times.
The dissipative nature of our protocol makes it inherently robust against perturbations such as imperfect control of the dwell times [27], and the preparation of the initial state.For the same reason, it will drive inhomogeneous nuclear ensembles towards an entangled state as long as the inhomogeneity is not too strong.In this sense, this type of protocol is expected to be more robust than coherent protocols such as the one presented in Ref. [72].

VI. CONCLUSIONS
In this work we propose a protocol for the deterministic generation of entanglement between two distant ensembles of nuclear spins, which are located around two QDs connected via a coherent electron bus.The protocol relies on the ability to inject electrons with definite polarization in the QDs and the coherent transfer of the electrons from one QD to the other.Controlling the time the electrons spend trapped in each QD, we engineer a bath that drives the system towards an entangled steady state.The dissipative nature of our protocol makes it inherently robust against imperfections in the dwell times and the electron transport.
Using a bosonization technique we have fully characterized the steady state in the homogeneous case within the large total angular momentum limit.Furthermore we have analyzed the dynamics in the inhomogeneous case using the EoMs for the covariance matrix.There, we have identified two different time scales characterizing the protocol: a faster one related to the collective hyperfine coupling on which the spins are driven into an entangled state, as evidenced by spin-squeezing of the collective nuclear spin operators, and a slower one related to the inhomogeneity, which leads to partial dephasing of the entanglement.
While the protocol requires the coherent transport of electron spins-a resource which also allows to distribute electron-spin entanglement-employing the nuclear spins makes an additional resource available for quantum information processing: in a hybrid approach, constantly generated nuclear-spin entanglement could be a resource accessed by the electron-spin quantum register at need.To take full advantage of the long nuclear coherence time, it would be appealing to store quantum information in the nuclear system, using the electron spins only to operate on the nuclei.Note that our protocol generalizes directly to the case of n spin ensembles.This is readily seen in the homogeneous and bosonic limit, where n linearly independent terms akin to the two non-local ones in the master equation Eq. (3) will stabilize an n-mode squeezed vacuum state, which comprise, for example, Gaussian cluster states that allow universal measurement-based quantum computation [73].Alternatively, pairs of entangled spin ensembles can be merged into multipartite entangled states by performing, for example, a measurement of a joint collective spin component of two of them; combined with full unitary control over the nuclear spins (as conceivable for few or single nuclear spins [24]) this also allows the preparation of computationally universal [74] multipartite entangled states.
while the ±-terms can be rewritten as Thus, after the first step of the protocol, the nuclei's reduced density matrix evolves as while after the second step, the nuclei's reduced density matrix evolves as The concatenation of these two steps yields: For the second part, we first work out the concatenation of m steps with |s in = |↓ , τ 1 g 1 = (µ 1 − ν 1 )/k and τ 2 g 2 = 0. Up to second order we find Similarly, the concatenation of k steps with |s in = |↑ , τ 1 g 1 = 0 and τ 2 g 2 = (µ 2 − ν 2 )/k, produces Thus, the concatenation of k steps of the first kind, followed by k steps of the second kind results in Last, putting together the equations for the first and second part of the protocol, all z-terms cancel and we get where we have defined ∆ρ n ≡ ρ n+2+2k − ρ n .Eq. ( 3) of the main text corresponds to the case where p = 0, µ 1 = µ 2 = µ and ν 1 = ν 2 = ν.

B. FIDELITIES FOR HOMOGENEOUS SYSTEMS
In Fig. B1, we show the fidelity of the steady state with respect to the two-mode spin squeezed state [Eq.( 7) in the main text] as a function of the number of steps in the cycle.The system consists of two identical homogeneous ensembles of spin-1/2 nuclei.The initial state belongs to the subspace with maximal total angular momenta J j = N j /2.As expected, the fidelity increases for increasing k.
2nd-order exact FIG.B1.Fidelity of the steady state with respect to the two-mode spin squeezed state for an homogeneous system of spin-1/2 nuclei as a function of the number of repetitions k performed in the last two steps of the protocol.The parameters of the system are N1 = N2 = N , µ = 0.5/ √ N and ν = 0.8µ.The steady state has been computed using the approximate 2nd-order master equation, Eq. ( 3) in the main text, and also using the exact time evolution superoperator, Eq. ( 2) in the main text.
For small systems with N j < 10 the EPR uncertainty may not detect entanglement if k is too small, see Fig. B2(a).Of course, making k larger brings ∆ EPR closer to the value in the ideal state, which does show entanglement irrespective of the number of nuclei considered.Nonetheless, according to the PPT criterion [64] the system is indeed entangled for all values of k and N , see Fig. B2(b).But, how can this entanglement be detected experimentally?An entanglement witness for the simplest case of just two qubits (N 1 = N 2 = 1, I = 1/2) can be constructed based on the following inequality, valid for separable states [48], Here, m ≡ J (m ≡ m ) is the expected value of the total spin, in this case J = I 11 + I 12 , and p S ≡ Ψ − | ρ |Ψ − is the singlet occupation.Since (1 − m 2 ) ≤ 1, we have that p S > 1/2 implies that the two spins are entangled.Thus, we define the following entanglement witness such that W < 0 for entangled states.In Fig. B2(c) we show that indeed this entanglement witness is able to detect the entanglement produced by our protocol.
FIG. B2.EPR uncertainty (a), minimum eigenvalue of the partial transpose (b), and entanglement witness based on the singlet occupation (c), for the steady state of a homogeneous system of spin-1/2 nuclei with parameters: N1 = N2 = N , µ = 0.1/ √ N , and ν = 0.8µ.For small k and N the EPR criterion fails to detect entanglement, although, the negativity of the partial transpose shows that the system is indeed entangled.The steady state has been computed using the approximate 2nd-order master equation (filled symbols), Eq. ( 3) in the main text, and also using the exact time evolution superoperator (black symbols), Eq. ( 2) in the main text.

C. HOLSTEIN-PRIMAKOFF TRANSFORMATION
In the limit of large total angular momenta and large polarizations, we can study the nuclear dynamics analytically replacing the collective total angular momentum operators by bosonic operators, and 2 ) in Eq. (A17).Doing so, defining μj ≡ µ j √ r j and νj ≡ ν j √ r j , we obtain the following master equation From this master equation, using the identities one can obtain the equation of motion of any operator expressed in terms of {a j , a † j } j=1,2 .Due to the bosonic commutation relations, the equation of motion of any product of these bosonic operators will only depend on products of the same or lower order.Hence, to find the steady-state expectation value of any set of monomials in these bosonic variables we will just have to solve a finite linear system of equations.We are interested in computing the EPR uncertainty, which depends on the absolute value polarizations | J z j | J j − a † j a j and variances var(J x/y 1 + J x/y 2 ) Thus, we have to compute the steady-state expectation values of all monomials up to second order.Their equations of motion are: where we have defined The adjoint operators evolve according to ∆X † = (∆X) † .By definition, the steady-state expectation value of an operator X, X ∞ , satisfies ∆ X ∞ = 0.For the operators we are considering, the only non-zero expectation values turn out to be Since many of the expectation values vanish, the formulas for the variances can be simplified as Substituting these back into the formulas for | J z j | and var(J x/y 1 + J x/y 2 ), we can obtain an analytical expression for ∆ EPR in the steady state.In Eq. ( 11) of the main text we show the analytical formula obtained for ∆ EPR in the case of an ideal protocol (p = 0, µ 1 = µ 2 = µ and ν 1 = ν 2 = ν).In Fig. C1 we show the values of ∆ EPR as a function of various parameters [cf.Fig. 2(a) of the main text].In the light of it we conclude that, at least for homogeneous ensembles, the protocol is able to produce entanglement for a wide range of values of the dwell times.These results also demonstrate that the protocol is robust against certain amount of decoherence during electron transport.
Moreover, from the corresponding EoMs, we can obtain the evolution of the expected number of bosons in each mode (and, therefore, the polarization of each ensemble).For an ideal protocol they are given by where we have defined κ ≡ (µν − ν 2 )/4 and η ≡ µ 2 /4 − κ.From these equations, we can deduce an analytical expression for the stabilization time In Fig. C2, we present results for the stabilization times of homogeneous systems.We show that this time depends on the initial state chosen, and for the fully polarized initial state |Ψ FP = |⇑ 1 |⇓ 2 it follows the prediction of the HP approximation.N and ν = 0.8µ, performed in an homogeneous system of spin-1/2 nuclei.Filled symbols correspond to the numerically obtained stabilization times for two different initial states, whereas the line corresponds to the formula T = 1/(2κ).

C.1. Mixed initial state
In Fig. 2(c) we consider a product initial state ρ 0 = ρ QD1 ⊗ ρ QD2 , where the initial state of each ensemble is mixed ρ QDj = J p j (J)ρ j (J).Here, ρ j (J) is a normalized density matrix that belongs to the jth-ensemble subspace of states with total angular momentum J. Normalization of ρ 0 requires J p j (J) = 1.As explained in the main text, for homogeneous systems the steady state will be given by ρ ∞ = J1,J2 p 1 (J 1 )p 2 (J 2 )ρ ∞ (J 1 , J 2 ), with ρ ∞ (J 1 , J 2 ) some steady-state within the (J 1 , J 2 )-subspace (the steady state is unique up to variations of the permutation quantum numbers).It can be computed numerically as the kernel of the Liouvillian corresponding to the evolution prescribed by Eq. ( 2) or the master equation Eq. ( 3), or analytically using the HP transformation.However, this latter approach would only yield a good approximation if the initial distributions p j (J) are peaked towards higher angular momenta.
In this manuscript we consider two different kinds of total angular momentum distributions.On the one hand, we consider an equal-weight mixture of all states with angular momenta J j ≥ J min j .In this case the probabilities p j (J) are given, up to a normalization constant, by where corresponds to the degeneracy of the Dicke states for non-maximal total angular momentum J. Since the number of states with total angular momentum J increases as J decreases, this probability distribution is peaked at J min j and decays towards higher values of J. On the other hand, we have considered Gaussian distributions with mean value J and variance w 2 , In Fig. C3, we compare both kinds of distributions.

D. NON-LINEAR EQUATIONS OF MOTION
In this section we show how to obtain the EoMs, and asses the validity of several different approximation schemes to close them.From the 2nd-order master equation, Eq. ( 3) in the main text, using the identity tr (BD(A)[ρ]) = [A † , B]A + A † [B, A] /2, we obtain the following general EoMs, valid for any operator X, Substituting X by σ + ij σ − i j and σ z ij σ z i j we obtain the equations of motion presented in the main text, Eqs. ( 12) and ( 13).However, these equations are quite involved and do not allow a clear vision of the evolution of the system.For this, we will rather investigate first the polarization dynamics in each ensemble separately, and then investigate the dynamics of the transverse spin variance.

D.1. Polarization dynamics
From the 2nd-order master equation, we find the following general equation of motion valid for any operator X 1 that acts non-trivially only on the nuclei belonging to the first ensemble.For an operator X 2 that acts non-trivially only on the nuclei belonging to the second ensemble, we obtain the same equation, replacing or, equivalently, changing κ → −κ and the subindex 1 → 2. We can realize that if the coupling constants {g i2 } and {g i1 } are the same, that is, if the two ensembles are identical, the steady-state polarization of the second ensemble must be the same as that of the first, but with opposite sign.This can be seen applying a transformation Thus, the EoM for Jz 2 would be the same as that for J z 1 (also with the same initial value if the initial state is |⇑ 1 |⇓ 2 ).Undoing the transformation, the steady-state polarization in QD2 is J To compute the polarization of the first ensemble, we particularize Eq. (D2) for X 1 = σ z i1 (i = 1, . . ., N 1 ), obtaining Here, and in the following equations, we have dropped the QD index to avoid cluttering.It should be understood that all operators and coupling constants refer to the first ensemble.As it turns out, the individual polarizations depend on all the matrix elements of the covariance matrix γ ij ≡ σ + i σ − j , which, in turn, depend on higher order correlations as Rewriting the terms in the second sum as those in the first, we obtain where λ ij = σ z i σ z j .Note that higher-order correlations can be simplified in case that some indices repeat, Therefore, The EoMs for the zz-correlations can be obtained in an analogous way, For the nuclei in the second ensemble we obtain equations similar to those presented so far, the only difference being the change of sign κ → −κ.
The different approximation schemes that we have considered are [54]: 1. Spin-temperature: Neglecting all correlations σ + i σ − j ≈ 0 for i = j, Eq. (D3) leads to an independent EoM for each nucleus Note that within this approximation the polarization of each nucleus in the first ensemble should be computed as σ z i = 3 − 2γ ii .For the second ensemble we substitute σ − i → a i and σ The polarization of the nuclei in the second ensemble can be computed using the usual formula σ z i = 2γ ii − 1.For the rest of schemes, we use Eqs.(D8) and (D9), approximating 3. Hartree: We may approximate as well the zz-correlations σ z i σ z j = 2 σ + i σ z j σ − i − σ z j , or use its own equation of motion, Eq. (D10).
In Fig. D1 we show the polarization dynamics of the first ensemble, obtained using the different approximation schemes 1-5.The Hartree, Wick, and Spin approximations are implemented without including the EoMs for the zz-correlations.As can be appreciated, these three approaches result in different dynamics that eventually converge to the same steady-state value as that obtained in the Spin-temperature description.If we include the EoMs for the zz-correlations, then the results are much better, see Fig. D2.Curiously enough, the best approximation for the first ensemble is the Spin approximation, while the best for the second is the Wick approximation.This can be understood as follows: In a rotated frame, we already established that the equation of motion of σz i2 is the same as that of σ z i1 .Applying the Spin approximation in this rotated frame amounts to Rotating back to the original frame we have Rearranging the operators, we then have which is nothing but the Wick approximation.We also note that the Boson approximation seems to be fine.However, as we show next, it is not so good for the variances.

D.2. Variance dynamics
Due to the symmetries of the protocol (we only send electrons polarized in the z-direction), we expect var(J x 1 + J x 2 ) = var(J y 1 + J y 2 ).The variance of the total nuclear spin along the x-direction is ij * (z * denotes the complex conjugate of z).However, since the Liouvillian superoperator and the initial nuclei's reduced density matrix are both real, the covariance matrix is symmetric at all times, σ + ij σ − i j = σ + i j σ − ij .This property of the covariance matrix is also maintained by the EoMs.Also, it can be shown that if σ + ij and σ + ij σ + i j are zero initially, they remain equal to zero throughout time evolution.With this assumption, noting that {σ + ij , σ − ij } = 1, we have We have worked out already the EoMs of γ 11 ij and γ 22 ij in the previous section, and found the best approximations to compute them, namely  12) and ( 13) in the main text].The coupling constants are all different and they have a relative standard deviation of 32% for the first ensemble, and 51% for the second; they are the same for all different initial states considered.The remaining parameters have values µ = 0.5/ √ N1 and ν = 0.8µ.

E. SHELL MODEL
We define the "shell" S j α as the set of nuclei i ∈ QDj with the same coupling constant g ij = g αj .We will use Greek indices to label different shells.As explained in the main text, the correlations are the same for all the nuclei inside a given shell: γ jj ii = x j α for all i ∈ S j α ; γ jj ii = y jj αβ for all i ∈ S j α , i ∈ S j β , (i, j) = (i , j ); and λ jj ii = z jj αβ for all i ∈ S j α , i ∈ S j β , (i, j) = (i , j ).Unless the coupling constants in a given QD are all different, the number of shells is smaller than the number of nuclei in that QD (α ≤ N j ).Thus, we can reduce the number of variables we need to consider in order to compute the polarizations and transverse spin variance.From Eqs. (D8-D10) (and the equivalent EoM for QD2), using the Z-rules approximation, we obtain the following set of equations: ∆x j α = −(−1) j κg αj 2S j α + g αj − ηg 2 αj 2x j α − 1 , (E1) ∆y jj αβ = (−1) j κ g βj 2x j β − 1 S j αβ + g αj 2x j α − 1 S j βα + g αj g βj x j α + x j β − 1 + 2κ g αj S j αβ + g βj S j βα y jj αβ − η g 2 αj + g 2 βj y jj αβ + ηg αj g βj z j αβ , (E2) ∆z jj αβ = −(−1) j 2κ g αj 2x j β − 1 2S j αβ + g αj + g βj 2x j α − 1 2S j βα + g βj − 8κ g βj S j αβ + g αj S j βα y jj αβ + 8ηg αj g βj y jj αβ − 2η g 2 αj + g 2 βj z jj αβ . (E3) Here, we have defined S j α = γ g γj N j γ y jj αγ − g αj y jj αα , and S j αβ = S j α − g βj y jj αβ ; N j α denotes the number of nuclei in shell S j α .Note that it only makes sense to consider y jj αα and z jj αα if N j α > 1, otherwise we should set y jj αα = z jj αα = 0.As for variables involving shells in different QDs, we obtain form Eqs. (D20, D21), using the Z-rules approximation, the following equations: The steady state is a solution of ∆x j α = 0, ∆y jj αβ = 0, and ∆z jj αβ = 0. Looking at this set of equations, we realize that ∆z jj αα = −4∆y jj αα , so the system is underdetermined.Nonetheless, at any given step n, z jj αα (n + 1) = z jj αα (n) + ∆z jj αα (n) = z jj αα (n) − 4∆y jj αα (n) , (E6) so the steady state solution satisfies z jj αα = z jj αα (0) − 4 y jj αα − y jj αα (0) , where y jj αα (0) and z jj αα (0) are the initial values of the corresponding variables.In this manner, the steady state actually depends on the initial state.Furthermore, we find numerically that for our particular choice of initial state, the steady state is such that y jj αβ = 0, for α = β.Substituting in back in Eqs.(E1-E3), we obtain Thus, the steady-state polarization of each shell is independent of the others, and it can be found solving a single quadratic equation.It depends on the values of µ, ν, and the number of nuclei in the shell N j α , but not on its coupling constant g αj .
Having found the steady state values of x j α and y jj αβ , we can substitute them back in Eqs.(E4, E5) and use a numerical algorithm to find the steady state values of y 12  αβ (we used a version of the Levenberg-Marquardt algorithm).As an example, in Fig. E1 we show the dynamics of a system where we model the electronic wavefunction with a Gaussian.The two ensembles consist of N 1 = 21 and N 2 = 25 nuclei arranged in a square lattice.In one of the cases considered, the electronic wavefunction in each QD is centered with respect to the nuclear lattice, so that many nuclei have the same coupling constant.Therefore, this case can be described with few shells containing several nuclei each.In the other case, the electronic wavefunction in QD1 has the same shape, but it has been displaced a little bit away from the center of the nuclear lattice, resulting in a system with many more shells, and with fewer nuclei per shell.The former case features a higher value of the nuclear polarizations, a lower value of the transverse spin variances, and a value of ∆ EPR < 1 in the steady state, which is reached within the computed time span.The second case, on the other hand, does not present an entangled steady state but metastable entanglement for many cycles.In Fig. E2, we compare the two entanglement witnesses ∆ and ∆ EPR for the system analyzed in Fig. 4 in the main text.We also show the evolution of some of the quantitites needed to compute them, which present interesting

FIG. 2 .
FIG. 2. (a) Contour plot of ∆ HPEPR , Eq. (11), for the subspace with maximal total angular momenta Jj = NjI (rj = 1).(b) Characterization of the steady state in the subspace with maximum total angular momenta as a function of the number of nuclei, for two identical (N1 = N2 ≡ N ), homogeneous ensembles of spin-1/2 nuclei.We compare the exact dynamics (circles) with the result given by the second-order master equation (squares), which does not depend on the individual values of µ and ν but only on their ratio.For all the curves plotted this ratio is the same, ν/µ = 0.8.The dashed horizontal line marks the value obtained by the Holstein-Primakoff approximation in the limit Jj → ∞, ∆ HP EPR .(c) Steady-state ∆EPR for an initial mixed state: blue circles and dots correspond to an initial equal-weight mixture of all states with J1, J2 ≥ Jmin, while orange squares and crosses correspond to an initial state with a Gaussian probability distribution, pj(J) ∝ exp[−(J − J) 2 /(2w 2 )], of the total angular momenta with mean J and variance w 2 = 3.The system is the same as in panel (b), with N = 40 spins.The ratio ν/µ = 0.8 is the same for all the points.The results have been obtained using the Holstein-Primakoff approximation and the 2nd-order master equation.

FIG. 4 .r max 1 = 2 .5 and r max 2 = 3 ,
FIG. 4. (a) Dynamics of a shell model with different inhomogeneity degrees (color-coded).The couplings follow a Gaussian function gij ∝ exp(− rij 2 /c), where {rij} denote the positions of the nuclei, which are assumed to be forming a 2D square lattice, rij ∈ Z 2 .We only consider the nuclei within a given distance from the origin rij < r max j , r max 1 = 2.5 and r max 2 = 3, resulting in a total of N1 = 21 and N2 = 25 nuclei in each ensemble.The inhomogeneity in each QD is determined by the ratio c/r max j , such that larger ratios correspond to more homogeneous systems.Parameters: µ = 0.5/ √ N1, ν = 0.8µ, k = 1.The initial state is |ΨFP .Dashed vertical lines roughly indicate the two different timescales present in the dynamics.Grey curves show the dynamics for an equivalent homogeneous system.(b) Stabilization times for the generalized polarization in QD1 vs the inverse squared smallest non-zero difference between the couplings in QD1.For very inhomogeneous systems the points obtained (not shown) do not follow this linear dependence.The color code used to denote different inhomogeneity degrees is the same in both panels.

1 = d 1 = 1 .
A13) where the coefficients obey a simple recurrence relation c m = c m−1 + 1 and d m = d m−1 + 1 + 2c m−1 for m > 1, and c Therefore, the concatenation of k such steps leads to

µ 1 = µ 2 , ν 1 = ν 2 1 1 N 10 3 10 4 T
FIG. C1.Dependence of the EPR uncertainty as given by the HP approximation on various parameters.In all the plots shown, k = 1, r1 = 0.9 and r2 = 0.8.Other fixed parameters have the values indicated in each plot.The orange dots in the panels of the left column correspond to the fixed values of N1/N2 and ν2/µ1 used in the other two panels.The orange squares in the panels of the right column correspond to the "ideal" choice µ1 = µ2 and ν1 = ν2.

FIG. D1 .
FIG.D1.Polarization dynamics for the first ensemble as given by the different approximation schemes.The plot on the right is just a zoom in the short-term dynamics.The results correspond to an homogeneous system with parameters: N1 = 50, µ = 0.5/ √ N1, and ν = 0.8µ.The Hartree, Wick, and Spin approximations have been used without including the EoMs for the zz-correlations.
FIG.D4.Evolution of the polarizations of each ensemble and variance of the total spin in the transverse plane for an inhomogeneous system with N1 = N2 = 5.Each column corresponds to a different initial state (ρ0 = |Ψ0 Ψ0|).Continuous lines show the results obtained using the 2nd-order master equation [Eq.(3) in the main text], while dashed lines show the results given by the EoMs [Eqs.(12) and (13) in the main text].The coupling constants are all different and they have a relative standard deviation of 32% for the first ensemble, and 51% for the second; they are the same for all different initial states considered.The remaining parameters have values µ = 0.5/ √ N1 and ν = 0.8µ.

− 1 ) − η y jj αα 2 +
2κ 2 (2N j α − 3) + 6η 2 y jj αα + κ 2 − η 2 = 0 .(E8) FIG. E1.(a) Schematics of the system under consideration.Each ensemble has a total of N1 = 21 and N2 = 25 nuclei, arranged in a lattice as shown by the black dots, with coupling constants proportional to a Gaussian function (colormap).The center of the Gaussian is indicated by a black cross.The case shown corresponds to the "centered" sytem.(b) Comparison between the dynamics of such system in case the Gaussian function is centered with respect to the lattice of nuclei, r0,j = (0, 0) (j = 1, 2) and in case it is shifted a little bit away from the center in QD1, r0,1 = (0.2, 0.1), r0,2 = (0, 0).In both cases the spread of the function is the same c = 10.The stars on the right side of the plot mark the steady state values of the corresponding quantities.The rest of parameters are in both cases µ = 0.5/ √ N1 and ν = 0.8µ, and the initial state is ρ0 = |ΨFP ΨFP|.