Quantum Coherence from Commensurate Driving with Laser Pulses and Decay

Non-equilibrium physics is a particularly fascinating field of current research. Generically, driven systems are gradually heated up so that quantum effects die out. In contrast, we show that a driven central spin model including controlled dissipation in a highly excited state allows us to distill quantum coherent states, indicated by a substantial reduction of entropy; the key resource is the commensurability between the periodicity of the pump pulses and the internal processes. The model is experimentally accessible in purified quantum dots or molecules with unpaired electrons. The potential of preparing and manipulating coherent states by designed driving potentials is pointed out.


Introduction
Controlling a quantum mechanical system in a coherent way is one of the long-standing goals in physics. Obviously, coherent control is a major ingredient for handling quantum information. In parallel, non-equilibrium physics of quantum systems is continuing to attract significant interest. A key issue in this field is to manipulate systems in time such that their properties can be tuned and changed at will. Ideally, they display properties qualitatively different from what can be observed in equilibrium systems. These current developments illustrate the interest in understanding the dynamics induced by time-dependent Hamiltonians H(t).
The unitary time evolution operator U (t 2 , t 1 ) induced by H(t) is formally given by where T is the time ordering operator. While the explicit calculation of U (t 2 , t 1 ) can be extremely difficult it is obvious that the dynamics induced by a time-dependent Hamiltonian maps quantum states at t 1 to quantum states at t 2 bijectively and conserves the mutual scalar products. Hence, if initially the system is in a mixed state with high entropy S > 0 it stays in a mixed state for ever with exactly the same entropy. No coherence can be generated in this way even for a complete and ideal control of H(t) in time. Hence, one has to consider open systems. The standard way to generate a single state is to bring the system of interest into thermal contact with a cold system. Generically, this is an extremely slow process. The targeted quantum states have to be ground states of some given system. Alternatively, optical pumping in general and laser cooling in particular [1] are well established techniques to lower the entropy of microscopic systems using resonant pumping and spontaneous decay. Quite recently, engineered dissipation has been recognized as a means to generate targeted entangled quantum states in small [2,3] and extended systems [4,5]. Experimentally, entanglement has been shown for two quantum bits [6,7] and for two trapped mesoscopic cesium clouds [8].
In this article, we show that periodic driving can have a quantum system converge to coherent quantum states if an intermediate, highly excited and decaying state is involved.
The key aspect is the commensurability of the periodic pump pulses to the internal process. This distinguishes our proposal from established optical pumping protocols. The completely disordered initial mixture can be made almost coherent. The final mixture only has an entropy S ≈ k B ln 2 corresponding to a mixture of two states. An appealing asset is that once the driving is switched off the Lindbladian decay does not matter anymore and the system is governed by Hamiltonian dynamics only.
The focus of the present work is to exemplarily demonstrate the substantial reduction of entropy in a small spin system subject to periodic laser pulses. The choice of system is motivated by experiments on the electronic spin in quantum dots interacting with nuclear spins [9][10][11][12][13][14][15][16]. The model studied is also applicable to the electronic spin in molecular radicals [17] or to molecular magnets, see Refs. [18][19][20]. In organic molecules the spin bath is given by the nuclear spins of the hydrogen nuclei in organic ligands.

Model
The model comprises a central, electronic spin S = 1/2 which is coupled to nuclear spins where H eZ = hS x is the electronic Zeeman term with h = gµ B B ( is set to unity) and the external magnetic field B in x-direction; H nZ = zh N i=1 I x i is the Zeeman term acting on the nuclear spins taken to be I = 1/2. Due to the large nuclear mass, the factor z is of the order of 10 −3 , but in principle other z-values can be studied as well. In the central spin part H CS = S · A the Overhauser field A results from the hyperfine interactions J i between the nuclei and the central spin We assume that the couplings J i are distributed evenly within a certain interval. Besides the spin system there is a single trion state |T polarised in z-direction at the very high energy ε (≈ 1 eV) so that the total Hamiltonian reads The laser pulse is taken to be very short as in experiment where it is of the order of picoseconds. Hence, we describe it as instantaneous unitary U puls which takes the | ↑ of the central spin to the trion state and vice versa where c := | ↑ T| and c † := |T ↑ |. Such pulses are applied in long periodic trains lasting seconds and minutes. The repetition time between two consecutive pulses is T rep of the order of 10 ns. The decay of the trion is described by the Lindblad equation for the density matrix ρ where the prefactor γ > 0 of the dissipator term [21] defines the decay rate. The corresponding process with c and c † swapped needs not be included because its decay rate is smaller by exp(−βε), i.e., it vanishes for all physical purposes.

Mathematical Properties of Time Evolution
The key observation is that the dynamics from just before the nth pulse at t = nT rep − to just before the n + 1st pulse at t = (n + 1)T rep − is a linear mapping M : ρ(nT rep −) → ρ((n + 1)T rep −) which does not depend on n. Since it is acting on operators one may call it a superoperator. Its matrix form is derived explicitly in Appendix A. If no dissipation took place (γ = 0) the mapping M would be unitary. But in presence of the dissipative trion decay it is a general matrix with the following properties: 1. The matrix M has an eigenvalue 1 which may be degenerate. If the dynamics of the system takes place in n separate subspaces without transitions between them the degeneracy is at least n.
2. All eigenoperators to eigenvalues different from 1 are traceless.
3. At least one eigenoperator to eigenvalue 1 has a finite trace.

4.
The absolute values of all eigenvalues of M are not larger than 1.

5.
If there is a non-real eigenvalue λ with eigenoperator C, the complex conjugate λ * is also an eigenvalue with eigenoperator C † .
6. The eigenoperators to eigenvalues 1 can be scaled to be hermitian.
While the above properties can be shown rigorously, see Appendix B, for any Lindblad evolution, the following ones are observed numerically in the analysis of the particular model (6) under study here: (a) The matrix M is diagonalizable; it does not require a Jordan normal form.
(c) The eigenoperators to eigenvalue 1 can be scaled to be hermitian and non-negative.
In the generic, non-degenerate case we denote the properly scaled eigenoperator V 0 with Tr(V 0 ) = 1.
(d) No eigenvalue = 1, but with absolute value 1, occurs, i.e., all eigenvalues different from 1 are smaller than 1 in absolute value.
(e) Complex eigenvalues and complex eigenoperators do occur.
The above properties allow us to understand what happens in experiment upon application of long trains of pulses corresponding to 10 10 and more applications of M . Then it is safe to conclude that all contributions from eigenoperators to eigenvalues smaller than 1 have died out completely. Only the (generically) single eigenoperator V 0 to eigenvalue 1 is left such that lim The quasi-stationary state after long trains of pulses is given by V 0 1 . This observation simplifies the calculation of the long-time limit greatly compared to previous quantum mechanical studies [12,13,16,22]. One has to compute the eigenoperator of M to the eigenvalue 1. Below this is performed by diagonalization of M which is a reliable approach, but restricted to small systems N 6. We stress that no complete diagonalization is required to know V 0 because only the eigenoperator to the eigenvalue 1 is needed. Hence we are optimistic that further computational improvements are possible. If, however, the speed of convergence is of interest more information on the spectrum and the eigenoperators of M is needed, see also Sect. 5.

Results on Entropy
It is known that in pulsed quantum dots nuclear frequency focusing occurs (NFF) [9,10,23] which can be explained by a significant change in the distribution of the Overhauser field [11-13, 15, 16, 22] which is Gaussian initially. This distribution develops a comb structure with equidistant spikes. The difference ∆A x between consecutive spikes is such that it corresponds to a full additional revolution of the central spin T rep ∆A x = 2π. A comb-like probability distribution is more structured and contains more information than the initial featureless Gaussian. For instance, the entropy reduction of the Overhauser field distributions computed in Ref. [16], Fig. 12, relative to the initial Gaussians is ∆S = −0.202k B at B = 0.93T and ∆S = −0.018k B at B = 3.71T. Hence, NFF decreases the entropy, but only slightly for large spin baths. This observation inspires us to ask to which extent continued pulsing can reduce entropy and which characteristics the final state has.
Inspired by the laser experiments on quantum dots [9,10,23] we choose an (arbitrary) energy unit J Q and thus /J Q as time unit which can be assumed to be of the order of 1ns. The repetition time T rep is set to 4π /J Q which is on the one hand close to the experimental values where T rep = 13.2ns and on the other hand makes it easy to recognize resonances, see below. The trion decay rate is set to 2γ = 2.5J Q to reflect a trion life time of ≈ 0.4ps. The bath size is restricted to N ∈ {1, 2, . . . , 6}, but still allows us to draw fundamental conclusions and to describe electronic spins coupled to hydrogen nuclear spins in small molecules [17][18][19][20]. The individual couplings J i are chosen to be distributed according to which is a uniform distribution between J min and J max with √ 5 inserted to avoid accidental commensurabilities. Results for a frequently used exponential parameterization [24] with α ∈ {0.5, 1} and for a Gaussian parametrization motivated by the electronic wave function in quantum dots [25] are given in the next section and in Appendix D. Figure 1 displays a generic dependence on the external magnetic field h = gµ B B x of the entropy of the limiting density matrix V 0 obtained after infinite number of pulses. Two nested resonances of the Larmor precessions are discernible: the central electronic spin resonates for hT rep = 2πn (n ∈ Z) while the nuclear bath spin resonates for zhT rep = 2πn . These conditions, however, apply only without pulsing or interactions. The driven system displays important shifts. The nuclear resonance appears to be shifted by z∆h ≈ ±J max /2, see right panel of Fig. 1(a). The explanation is that the dynamics of the central spin S = 1/2 creates an additional magnetic field acting on each nuclear spin of the order of J i /2 which is estimated by J max /2. Further support of this explanation is given in Appendix C.
The electronic resonance is shifted by where A max is the maximum Overhauser field given by A max := (1/2) N i=1 J i for maximally polarized bath spins. This is shown in the right panel of Fig. 1(b). The commensurability of these resonances is crucial for the advocated mechanism.
For the minimum entropy, it is essential not to look at the bare resonances, but to take the two above mentioned shifts into account. We posed the question to which extent the initial entropy of complete disorder S init = k B (N + 1) ln 2 (in the figures and henceforth k B is set to unity) can be reduced by commensurate periodic pumping. The results in Fig.  1 clearly show that remarkably low values of entropy can be reached. The residual value of S ≈ 0.5k B in the minima of the right panel of Fig. 1(b) corresponds to a contribution of less than two states (S = ln 2k B ≈ 0.7k B ) while initially 16 states were mixed for N = 3 so that the initial entropy is S init = 4 ln 2k B ≈ 2.77k B . This represents a remarkable distillation of coherence.
Hence, we focus on the minima and in particular on the left minimum. We address the question whether the distillation of coherence still works for larger systems. Unfortunately, the numerical analysis cannot be extended easily due to the dramatically increasing dimension D = 2 2(N +1) because we are dealing with the Hilbert space of density matrices of the spin bath and the central spin. Yet a trend can be deduced from results up to N = 6 displayed in Fig. 2(a). The entropy reduction per N + 1 spins is −0.58k B for N = 3, −0.57k B for N = 4, −0.55k B for N = 5, and −0.52k B for N = 6. The reduction is substantial, but slowly decreases with system size. Presently, we cannot know the behavior for N → ∞. The finite value ≈ −0.2k B found in the semiclassical simulation [15,16] indicates that the effect persists for large baths. In Appendix D, results for the couplings defined in (9) or in (10) are given which corroborate our finding. The couplings may be rather close to each other, but not equal. It appears favorable that the spread of couplings is not too large.
Which state is reached in the minimum of the residual entropy? The decisive clue is provided by the lower panel Fig. 2(b) displaying the polarization of the spin bath. It is normalized such that its saturation value is unity. Clearly, the minimum of the residual entropy coincides with the maximum of the polarization. The latter is close to its saturation value though not quite with a minute decrease for increasing N . This tells us that the limiting density matrix V 0 essentially corresponds to the polarized spin bath. The central electronic spin is also almost perfectly polarized (not shown), but in z-direction. These observations clarify the state which can be retrieved by long trains of pulses. Additionally, Fig. 2(b) explains the shift of the electronic resonance. The polarized spin bath renormalizes the external magnetic field by (almost) ±A max . To the left of the resonance, it enhances the external field (+A max ) while the external field is effectively reduced (−A max ) to the right of the resonance. Note that an analogous direct explanation for the shift of the nuclear resonance in the right panel of Fig. 1 is not valid. The computed polarization of the central spin points in z-direction and thus does not shift the external field.

Results on Convergence
In order to assess the speed of convergence of the initially disordered density matrix ρ 0 = 1/Z to the limiting density matrix V 0 we proceed as follows. Let us assume that the matrices v i are the eigen matrices of M and that they are normalized ||v i || 2 := Tr(v † i v i ) = 1. Since the mapping M is not unitary, orthogonality of the eigenmatrices cannot be assumed. Note that the standard normalization generically implies that there is some factor between V 0 with Tr(V 0 ) = 1 and v 0 . The initial density matrix ρ 0 can be expanded in the {v i } After n pulses, the density matrix ρ n is given by where λ j are the corresponding eigenvalues of M and λ 0 = 1 by construction. We aim at ρ 0 being close to V 0 within p thresh , i.e., should hold for an appropriate n. A generic value of the threshold p thresh is 1%. To this end, the minimum n which fulfills (14) has to be estimated.  Such an estimate can be obtained by determining for j ∈ {1, 2, 3, . . . , D − 1}. The estimate of the required number of pulses is the maximum of these number, i.e., n puls := max 1≤j<D n j .
We checked exemplarily that the number determined in this way implies that the convergence condition (14) is fulfilled. This is not mathematically rigorous because it could be that there are very many slowly decreasing contributions which add up to a significant deviation from V 0 . But generically, this is not the case. In Fig. 3 the results are shown for various bath sizes and the parameters for which the data of the previous figures was computed. Since the entropy minima are located at the positions of the vertical dashed lines to good accuracy one can read off the required number of pulses at the intersections of the solid and the dashed lines. Clearly, about 2 · 10 12 pulses are necessary to approach the limiting, relatively pure density matrices V 0 . Interestingly, the number of required pulses does not depend much on the bath size, at least for the accessible bath sizes. This is a positive message in view of the scaling towards larger baths in experimental setups.    (9) and (10). Again, the range is about 3 · 10 12 . Still, there are relevant differences. The value n puls is higher for α = 1 (≈ 4 · 10 12 ) than for α = 1/2 ( 2 · 10 12 ). This indicates that the mechanism of distilling quantum states by commensurability with periodic external pulses works best if the couplings are similar, i.e., if their spread is small. The same qualitative result was obtained for the residual entropy, see Appendix D. Note that this argument also explains why the Gaussian parametrized couplings (10) require slightly less pulses than the exponential parametrized couplings (9). One could have thought that the cumulated couplings J i ≈ J max in the Gaussian case require longer pulsing in order to achieve a given degree of distillation because mathematically equal couplings J i = J i imply degeneracies preventing distillation, see the mathematical properties discussed in Sect. 3. But this is obviously not the case.
The total numbers of pulses is rather high. As many as 2 · 10 12 pulses for a repetition time T rep ≈ 10ns imply about six hours of pulsing. This can be achieved in the lab, but the risk that so far neglected decoherence mechanisms spoil the process is real. If, however, the pulses can be applied more frequently, for instance with T rep = 1ns, the required duration shrinks to about 30 minutes. The question arises why so many pulses are required. While a comprehensive study of this aspect is beyond the scope of the present article, first clue can be given.
It suggests itself that the slow dynamics in the bath is responsible for the large number of pulses required for convergence. This idea is corroborated by the results displayed in Fig. 5 where a larger maximum coupling and, importantly, a larger z factor is assumed. Recall that the z-factor is the ratio of the Larmor frequency of the bath spins to the Larmor frequency of the central spin. If it is increased, here by a factor of 100, the bath spins precess much quicker. Indeed, the range of the required number of pulses is much lower with 2 · 10 7 which is five orders of magnitude less than for the previous parameters. The former six hours then become fractions of seconds. Of course, the conventional g-factors of nuclear and electronic spins do not allow for z = 0.1. But the central spin model as such, built by a central spin and a bath of spins supplemented by a damped excitation  (8) (legend "equidist"), the exponential parametrization in (9) (legend "expo") and the Gaussian parametrization in (10) (legend "gaus"). The vertical dashed lines indicate the estimates for the entropy minima which are shifted from the resonances without interactions according to (11).
can also be realized in a different physical system.

Conclusion
Previous work has established dynamic nuclear polarization (DNP), for a review see Ref. [26]. But it must be stressed that the mechanism of this conventional DNP is fundamentally different from the one described here. Conventionally, the polarization of an electron is transferred to the nuclear spins, i.e., the polarization of the electrons induces polarization of the nuclei in the same direction. In contrast, in the setup studied here, the electron is polarized in z-direction while the nuclear spins are eventually polarized perpendicularly in x-direction. Hence, the mechanism is fundamentally different: it is NFF stemming essentially from commensurability. This is also the distinguishing feature compared to standard optical pumping. States in the initial mixture which do not allow for a time evolution commensurate with the repetition time T rep of the pulses are gradually suppressed. Eventually, only the particular state which allows for a dynamics commensurate with T rep persists. This mechanism can be used also for completely different physical systems, e.g., in ensembles of oscillators. The studied case of coupled spins extends the experimental and theoretical observations of NFF for large spin baths [9][10][11][12][13][14][15][16] where many values of the polarization of the Overhauser field can lead to commensurate dynamics. Hence, only a partial reduction of entropy occurred.
The above established DNP by NFF comprises the potential for a novel experimental technique for state preparation: laser pulses instead of microwave pulses as in standard NMR can be employed to prepare coherent states which can be used for further processing, either to perform certain quantum protocols or for analysis of the systems under study. The combination of optical and radio frequency pulsing appears promising because it enlarges the possibilities of experimental manipulations. Another interesting perspective is to employ the concept of state distillation by commensurability to physical systems other than localized spins, for instance to spin waves in quantum magnets. A first experimental observations of commensurability effects for spin waves in ferromagnets are already carried out [27].
In summary, we showed that dissipative dynamics of a highly excited state is sufficient to modify the dynamics of energetically low-lying spin degrees of freedom away from unitarity. The resulting dynamic map acts like a contraction converging towards a single density matrix upon iterated application. The crucial additional ingredient is commensurability which enables a substantial entropy reduction, almost to a single pure state. This has been explicitly shown for an exemplary small central spin model including electronic and nuclear Zeeman effect. This model served as proof-of-principle model to establish the mechanism of distillation by commensurability.
Such a model describes the electronic spin in quantum dots with diluted nuclear spin bath or the spin of unpaired electrons in molecules, hyperfine coupled to nuclear hydrogen spins. We stress that the mechanism of commensurability can also be put to use in other systems with periodic internal processes. The fascinating potential to create and to manipulate coherent quantum states by such approaches deserves further investigation.

Acknowledgements
The author thanks A. Greilich, J. Schnack, and O. P. Sushkov for useful discussions and the School of Physics of the University of New South Wales for its hospitality.
Funding information This work was supported by the Deutsche Forschungsgemeinschaft (DFG) and the Russian Foundation of Basic Research in TRR 160, by the DFG in project no. UH 90-13/1, and by the Heinrich-Hertz Foundation of Northrhine-Westfalia.

A Derivation of the Linear Mapping
The goal is to solve the time evolution of ρ(t) from just before a pulse until just before the next pulse. Since the pulse leads to a unitary time evolution which is linear (17) with U puls from (5) and the subsequent Lindblad dynamics defined by the linear differential equation (6) is linear as well the total propagation in time is given by a linear mapping M : ρ(nT rep −) → ρ((n + 1)T rep −). This mapping is derived here by an extension of the approach used in Ref. [16].
The total density matrix acts on the Hilbert space given by the direct product of the Hilbert space of the central spin comprising three states (| ↑ , | ↓ , |T ) and the Hilbert space of the spin bath. We focus on ρ TT := T|ρ|T which is a 2 N × 2 N dimensional density matrix for the spin bath alone because the central degree of freedom is traced out. By ρ S we denote the d × d dimensional density matrix of the spin bath and the central spin, i.e., d = 2 N +1 since no trion is present: ρ S |T = 0. The number of entries in the density matrix is D = d 2 , i.e., the mapping we are looking for can be represented by a D × D matrix.
The time interval T rep between two consecutive pulses is sufficiently long so that all excited trions have decayed before the next pulse arrives. In numbers, this means 2γT rep 1 and implies that ρ(nT rep −) = ρ S (nT rep −) and hence inserting the unitary of the pulse (5) yields where we used the standard ladder operators S ± of the central spin to express the projection | ↓ ↓ |. The equations (18) set the initial values for the subsequent Lindbladian dynamics which we derive next. For completeness, we point out that there are also nondiagonal contributions of the type T|ρ| ↑ , but they do not matter for M . Inserting ρ TT into the Lindblad equation (6) yields No other parts contribute. The solution of (19) reads By the argument 0+ we denote that the initial density matrix for the Lindbladian dynamics is the one just after the pulse. For ρ S , the Lindblad equation (6) implies Since we know the last term already from its solution in (20) we can treat it as given inhomogeneity in the otherwise homogeneous differential equation. With the definition U S (t) := exp(−iH spin t) we can write Integration leads to the explicit solution If we insert (20) into the above equation we encounter the expression where I x tot := S x + N i=1 I x i is the total momentum in x-direction. It is a conserved quantity commuting with H spin so that a joint eigenbasis with eigenvalues m α and E α exists. We determine such a basis {|α } by diagonalization in the d-dimensional Hilbert space (d = 2 N +1 ) of central spin and spin bath and convert (23) in terms of the matrix elements of the involved operators. For brevity, we write ρ αβ for the matrix elements of ρ S .
Elementary quantum mechanics tells us that with a := zht /2 which we need for the last row of equation (25). Replacing ρ TT (0+) by ↑ |ρ S (nT rep −)| ↑ according to (18b) and inserting (26) we obtain with the three d × d matrices In this derivation, we expressed ket-bra combinations by the spin ladder operators according to The final step consists in inserting (27b) into (25) and integrating the exponential time dependence straightforwardly from 0 to T rep . Since we assume that 2γT rep 1 so that no trions are present once the next pulse arrives the upper integration limit T rep can safely and consistently be replaced by ∞. This makes the expressions appear where τ ∈ {−1, 0, 1}. Finally, we use (18c) and summarize This provides the complete solution for the dynamics of d × d matrix ρ S from just before a pulse (t = 0−) till just before the next pulse for which we set t = T rep in (31). In order to set up the linear mapping M as D × D dimensional matrix with D = d 2 we denote the matrix elements M µ µ where µ is a combined index for the index pair αβ and µ for α β with α, β, α , β ∈ {1, 2 . . . , d}. For brevity, we introduce Then, (31) implies This concludes the explicit derivation of the matrix elements of M . Note that they are relatively simple in the sense that no sums over matrix indices are required on the right hand side of (33). This relative simplicity is achieved because we chose to work in the eigenbasis of H spin . Other choices of basis are possible, but render the explicit respresentation significantly more complicated.

B Properties of the Time Evolution
Preliminaries Here we state several mathematical properties of the mapping M which hold for any Lindblad dynamics over a given time interval which can be iterated arbitrarily many times. We assume that the underlying Hilbert space is d dimensional so that M acts on the D = d 2 dimensional Hilbert space of d × d matrices, i.e., M can be seen as D × D matrix. We denote the standard scalar product in the space of operators by where the trace refers to the d × d matrices A and B.
Since no state of the physical system vanishes in its temporal evolution M conserves the trace of any density matrix Tr(M ρ) = Tr(ρ).
This implies that M conserves the trace of any operator C. This can be seen by writing C = (C + C † )/2 + (C − C † )/2 = R + iG where R and G are hermitian operators. They can be diagonalized and split into their positive and their negative part R = p 1 − p 2 and G = p 3 − p 4 . Hence, each p i is a density matrix up to some real, positive scaling and we have where 1 d is the d×d-dimensional identity matrix and M † is the D ×D hermitian conjugate of M . From (38) we conclude which means that 1 d is an eigenoperator of M † with eigenvalue 1. Since the characteristic polynomial of M is the same as the one of M † up to complex conjugation we immediately see that 1 is also an eigenvalue of M . If the dynamics of the system takes place in n independent subspaces without transitions between them, the n different traces over these subspaces are conserved separately Such a separation occurs in case conserved symmetries split the Hilbert space, for instance the total spin is conserved in the dynamics given by (6) if all couplings are equal. Then, the above argument implies the existence of n different eigenoperators with eigenvalue 1. Hence the degeneracy is (at least) n which proves property 1. in the main text.
implies Tr(C) = 0, i.e., tracelessness as stated. Since all density matrices can be written as linear combinations of eigenoperators there must be at least one eigenoperator with finite trace. In view of property 2., this needs to be an eigenoperator with eigenvalue 1 proving property 3. The latter conclusion holds true if we assume that M cannot be diagonalized, but only has a Jordan normal form. If d J is the dimension of the largest Jordan block, the density matrix M d J −1 ρ will be a linear combination of eigenoperators while still having the trace 1.
Property 4. Next, we show that no eigenvalue λ can be larger than 1 in absolute value. The idea of the derivation is that the iterated application of M to the eigenoperator belonging to |λ| > 1 would make this term grow exponentially ∝ |λ| n beyond any bound which cannot be true. The formal proof is a bit intricate. First, we state that for any two density matrices ρ and ρ their scalar product is nonnegative 0 ≤ (ρ|ρ ) because it can be viewed as expectation value of one of them with respect to the other and both are positive operators. In addition, the Cauchy-Schwarz inequality implies Let C be the eigenoperator of M † belonging to λ; it may be represented as in (36) and scaled such that the maximum of the traces of the p i is 1. Without loss of generality this is the case for p 1 , i.e., Tr(p 1 ) = 1. Otherwise, C is simply rescaled: by C → −C to switch p 2 to p 1 , by C → −iC to switch p 3 to p 1 , or by C → iC to switch p 4 to p 1 . On the one hand, we have for any density matrix ρ n where the last inequality results form (41). On the other hand, we set ρ n := M n p 1 and obtain 2 ≥ |(C|ρ n )| = |((M † ) n C|p 1 )| = |λ * | n |(C|p 1 )| = |λ| n ( (C|p 1 )) 2 + ( (C|p 1 )) 2 (43a) ≥ |λ| n | (C|p 1 )| = |λ| n (p 1 |p 1 ) (43b) where we used (p 1 |p 2 ) = 0 in the last step; this holds because p 1 and p 2 result from the same diagonalization, but refer to eigenspaces with eigenvalues of different sign. In essence we derived 2 ≥ |λ| n (p 1 |p 1 ) which clearly implies a contradiction for n → ∞ because the right hand side increases to infinity for |λ| > 1. Hence there cannot be eigenvalues with modulus larger than 1.
where ρ 0 is an arbitrary initial density matrix which should contain contributions from all eigenspaces of M . For instance, a Gram-Schmidt algorithm applied to the Krylov basis generates an orthonormal basisρ n . Due to the fact, that all the operators ρ n from (45) are hermitian density matricesρ n =ρ † n , we know that all overlaps (ρ m |ρ n ) are real and hence the constructed orthonormal basisρ n consists of hermitian operators. Also, all matrix elements (ρ m |M ρ n ) = (ρ m |ρ n+1 ) are real so that the resulting representationM is a matrix with real coefficients whenceM by complex conjugation. Here c is a vector of complex numbers c n which define the corresponding eigenoperators by Thus, c and c * define C and C † , respectively.
Property 6. In view of the real representationM of M with respect to an orthonormal basis of hermitian operators derived in the previous paragraph the determination of the eigenoperators with eigenvalue 1 requires the computation of the kernel ofM −1 D . This is a linear algebra problem in R D with real solutions which correspond to hermitian operators by means of (47). This shows the stated property 6..

C Shift of the Nuclear Resonance
In the main text, the shift of the nuclear resonance due to the coupling of the nuclear spins to the central, electronic spin was shown in the right panel of Fig. 1(a). The effect can be estimated by z∆h ≈ ±J max /2. This relation is highly plausible, but it cannot be derived analytically because no indication for a polarization of the central, electronic spin in x-direction was found. Yet, the numerical data corroborates the validity of (48).
In Fig. 6, we show that the nuclear resonance without shift occurs for where n ∈ Z. But it is obvious that an additional shift occurs which is indeed captured by (48). In order to support (48) further, we also study various values of J max in Fig. 7. The estimate (48) captures the main trend of the data, but it is not completely quantitative because the position of the dashed lines relative to the minimum of the envelope of the resonances varies slightly for different values of J max . Hence, a more quantitative explanation is still called for.

D Entropy Reduction for Other Distributions of Couplings
In the main text, we analyzed a uniform distribution of couplings, see Eq. (8). In order to underline that our results are generic and not linked to a special distribution, we provide additional results for two distributions which are often considered in literature, namely an exponential parameterization as defined in (9) and a Gaussian parametrization as defined in (10).
The key difference between both parametrizations (9) and (10) is that due to the quadratic argument in (10) the large couplings in this parametrization are very close to each other, in particular for increasing N . Hence, one can study whether this feature is favorable of unfavorable for entropy reduction.
Additionally, the difference between α = 0.5 and α = 1 consists in a different spread of the couplings. For α = 1, one has J min /J max = 1/e in both parametrizations while one has J min /J max = 1/ √ e for α = 0.5, i.e., the spread is smaller. Figure 8 displays the results for the exponential parametrization (9) while Fig. 9 depicts the results for the Gaussian parametrization (10). Comparing both figures shows that the precise distribution of the couplings does not matter much. Exponential and Gaussian parametrization lead to very similar results. They also strongly ressemble the results shown in Fig. 2a in the main text for a uniform distribution of couplings. This is quite remarkable since the Gaussian parametrization leads to couplings which are very close to each other and to the maximum coupling. This effect does not appear to influence the achievable entropy reduction.
The ratio J min /J max between the smallest to the largest coupling appears to have an impact. If it is closer to unity, here for α = 0.5, the reduction of entropy works even better than for smaller ratios.