Bath-induced Zeno localization in driven many-body quantum systems
Thibaud Maimbourg, Denis M. Basko, Markus Holzmann, Alberto Rosso
BBath-induced Zeno localization in driven many-body quantum systems
Thibaud Maimbourg, ∗ Denis M. Basko, Markus Holzmann, and Alberto Rosso LPTMS, CNRS, Universit´e Paris-Saclay, 91405, Orsay, France. Universit´e Grenoble Alpes and CNRS, LPMMC, 25 rue des Martyrs, 38042 Grenoble, France
We study a quantum interacting spin system subject to an external drive and coupled to a thermalbath of spatially localized vibrational modes, serving as a model of Dynamic Nuclear Polarization.We show that even when the many-body eigenstates of the system are ergodic, a sufficiently strongcoupling to the bath may effectively localize the spins due to many-body quantum Zeno effect, asmanifested by the hole-burning shape of the electron paramagnetic resonance spectrum. Our resultsprovide an explanation of the breakdown of the thermal mixing regime experimentally observedabove 4 – 5 Kelvin.
The fate of an isolated many-body quantum system, i.e. thermalization or ergodicity breaking, depends onthe statistical properties of its eigenstates [1]. Indeed,the unitary dynamics projects, through dephasing, theinitial state on the eigenstates of the Hamiltonian. Ac-cording to the Eigenstate Thermalization Hypothesis(ETH), every eigenstate is representative of the ensem-ble; namely, the expectation value of physical observ-ables in an eigenstate depends only upon a few con-served quantities (such as the total energy) and coin-cides with the prediction of the microcanonical ensem-ble [2–4]. It was recently realized that ETH may breakdown due to many-body localization (MBL, see recentreviews [5–7] and references therein), when Anderson lo-calization in the many-body Hilbert space prevents theeigenstates from spreading ergodically over the availableportion of the Hilbert space [8, 9].In open systems, coupling to an external bath per-mits thermalization even when eigenstates are localized,as manifested by phonon-induced hopping transport oflocalized electrons [10]: the bath supplies or absorbsthe energy needed to allow hopping between localizedstates, thereby destroying the localization. In this Let-ter, we show how coupling to a local bath can instead in-duce localization in a quantum many-body system withergodic eigenstates, as revealed by the properties of thenonequilibrium steady state reached under an externaldrive. Such localization is not due to Anderson, butrather to the quantum Zeno effect [11–17] in the many-body Hilbert space. This scenario should not be con-fused with the quantum Zeno effect in quantum gaseswith localized particle losses [18–21] or dephasing [22],where the combined effect of local single-particle lossesand many-body interactions is mainly to renormalizesingle-particle quantities. It should also not be mistakenwith localization induced by subohmic baths at zerotemperature [23–26], essentially a polaronic effect de-stroyed by interactions or finite temperature [27]. Thelocalization studied here is more akin to the entan-glement transition between volume-law and area-lawphases found in schematic models with random localprojective measurements such as quantum circuits [28–36], free fermionic chains [37] and interacting bosonicchains [38–40].We use an out-of-equilibrium driven-dissipative pro-tocol: it has recently been emphasized [41–46] and ex-perimentally confirmed [47] that fingerprints of ther-malization or localization appear in the steady state ofsystems weakly coupled to a bath if an external drive is ∗ [email protected] FIG. 1. (a) Sketch of the system: N electron spins withdipolar interactions of strength U ij in a strong inhomoge-neous magnetic field ω e + ∆ i in contact with a thermostatare irradiated by microwaves at frequency ω MW (wavy ar-row). (b)(c) EPR spectrum f ( ω ) at Boltzmann equilibrium(blue curve) and under microwave driving displaying twodifferent shapes (orange curves): (b) the spin-temperatureprofile [48, 49] with a linear behavior (dashed red) of slope β s / ω = ω MW ) (c) the holeburning at microwave resonance. applied. In fact, if the transitions rates induced by thebath or drive are small with respect to dephasing, thedynamics still gets projected in the Hamiltonian eigen-states. For thermal systems, the steady state is thena mixture of ETH-satisfying eigenstates. Thus, evenin the presence of a drive, it is described by an effec-tive equilibrium with a unique temperature. On thecontrary, driven MBL systems end up in a non-thermalsteady state.Through a driven protocol we show that a bathof spatially-localized modes triggers spatially-localizedquantum jumps of the system’s state during its evo-lution which can induce localization in ETH systems.Indeed if the bath transition rates prevail over dephas-ing, the projection on eigenstates is hindered, prevent-ing the system from reaching a thermal stationary state.Consequently an effective equilibrium description of thesteady state breaks down. Note that the bath transitionrates are generally increasing with temperature; thus,curiously, localization happens in a high-temperaturephase. The phenomenon is analog to the quantumZeno effect, when the unitary evolution is impeded byinfinitely frequent measurements of a coupled probe.Here the spatially-localized quantum jumps play the a r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p FIG. 2. (left) A system of three levels | (cid:105) , | (cid:105) , | (cid:105) with ener-gies 0 , ω e ± ∆ and a coupling J between | (cid:105) and | (cid:105) . Amonochromatic drive at frequency ω MW couples | (cid:105) and | (cid:105) . Two independent zero-temperature baths lead to de-cay of the population from | (cid:105) and | (cid:105) to | (cid:105) with rate 2 γ .The largest energy scale is ω e . (right) Level populations ρ ii ( ω MW − ω e ) in three limits. In the Anderson limit thestrong imbalance between ρ and ρ is caused by the eigen-states’ localization on the respective levels. In the Zeno limitthe system eigenstates are uniformly spread over | (cid:105) and | (cid:105) ,but they do not have time to form because of the dissipation. role of the measurements. The essence of the differ-ence between the Zeno and the Anderson localizationin a driven-dissipative system can be illustrated by asimple example with just three levels, shown in Fig. 2.We investigate these ideas in the problem of Dy-namic Nuclear Polarization (DNP) [49] induced by thenonequilibrium steady state of N diluted electron spinsexposed to a strong magnetic field ω e 1 along the z axis.Their Hamiltonian reads (see Fig. 1(a))ˆ H S = N (cid:88) i =1 ( ω e + ∆ i ) ˆ S zi + ˆ H dip (1)Here ∆ i is a small disorder from the random orientationof the molecule where the spin lies and ˆ H dip stands forthe dipolar interaction. The large value of the magneticfield implies ˆ S z = (cid:80) i ˆ S zi is conserved, hence the dipolarHamiltonian gets truncated as [42, 49–51]ˆ H dip = (cid:88) i The electronspins are dilute, so we assume they are in contact withindependent spatially-localized vibration modes ˆ B µi :ˆ H int = (cid:88) i =1 ,...,Nµ = x,y,z ˆ S µi ⊗ ˆ B µi (3)The quantum dynamics of the spin and bath degrees offreedom amounts to an intractable unitary evolution ofthe full density matrix ρ S ⊗ B . Nonetheless, assuming thebath is equilibrated at temperature β − , one may traceout the ˆ B µi variables and write an effective evolutionfor the spin system density matrix ρ = Tr B ( ρ S ⊗ B ). Theensuing evolution is no longer unitary but must still pre-serve the trace and semi-positivity of ρ . The most gen-eral Markovian dynamics must then be of the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) form [16]:˙ ρ = − i [ ˆ H, ρ ] + (cid:88) α ˆ A α ρ ˆ A † α − (cid:110) ˆ A † α ˆ A α , ρ (cid:111) (4)where ˆ H is Hermitian and { ˆ A α } is a set of jumpoperators. To integrate the bath degrees of free-dom [SM,Sec.I], we consider the spin-bath couplingweak and perform a perturbative expansion of the fullunitary dynamics of ρ S ⊗ B at second order in ˆ H int .Within the Born-Markov approximation [16, 56, 57],one can turn the perturbative expansion into an effec-tive Markovian evolution for ρ . As the bath degrees See the Supplementary Material [SM] for further details. of freedom are independent, they enter the dynamicalequation only through a single equilibrium correlationfunction γ ( ω ) = ˆ ∞−∞ d τ e iωτ (cid:68) ˆ B µi ( τ ) ˆ B µi (0) (cid:69) B = h ( ω ) T ( | ω | ) . (5)The hypothesis of an equilibrated bath implies the lastequality, in which h ( ω ) = (1+ e − βω ) − enforces detailedbalance at inverse temperature β , while T ( | ω | ) is thetimescale of the exchange of energy ω with the spins.The way to implement the Markovian approximationis not unique, yielding different dynamical equationsin general not of the GKSL form (4). A Markovianprescription was recently proposed [58–60], which leadsto a GKSL form by setting the unitary part ˆ H = ˆ H S 3 and the jump operators of the dissipative partˆ A α = (cid:88) n,m (cid:112) γ ( ω nm ) (cid:104) m | ˆ S µi | n (cid:105) | m (cid:105)(cid:104) n | (6)with α = ( i, µ ) and ω nm = ε n − ε m are energy gaps ofˆ H S . We have three characteristic timescales T ( | ω | ) cor-responding to three energy scales: (i) T for transitionsof energy gap ± ω e , providing the jump operatorsˆ A xi = (cid:115) h ( ω e )2 T (cid:16) ˆ S − i + e − βω e / ˆ S + i , (cid:17) ˆ A yi = i (cid:115) h ( ω e )2 T (cid:16) ˆ S − i − e − βω e / ˆ S + i (cid:17) , (7) (ii) T ∗ for transitions of finite energy | ω | (cid:28) ω e , and (iii) T (0) for zero-energy transitions within the sameeigenstate, givingˆ A zi = ˆ S zi √ T ∗ + (cid:32) (cid:112) T (0) − √ T ∗ (cid:33) (cid:88) n (cid:104) n | ˆ S zi | n (cid:105) | n (cid:105)(cid:104) n | . (8)The operators (7) and (8) (when T (0) ≈ T ∗ ) are welllocalized in space and called nonsecular jump operators.The quantum trajectories thus result from a competi-tion between the unitary dynamics which, through de-phasing, projects the system’s state on thermal eigen-states, and repeated measurements performed by thenonsecular jump operators. If the jump rates domi-nate the dephasing, the system’s projection on thermaleigenstates is hence hampered in a way reminiscent ofthe quantum Zeno effect.This choice of jump operators contrasts with theMarkovian prescription employed in the usual weak-coupling scheme [16, 56, 57]. Within this approach, theGKSL form of master equation (4) is actually recoveredthrough an additional secular approximation where thejump operators select only a given transition of energy ω nm between eigenstates | m (cid:105) and | n (cid:105) :ˆ A sec α ( ω nm ) = (cid:112) γ ( ω nm ) (cid:104) m | ˆ S µi | n (cid:105) | m (cid:105)(cid:104) n | (9)Note that the secular jump operators (9) are projectedon the eigenbasis of ˆ H S ; if the eigenstates satisfy ETH,the jump operators are delocalized in space. Inserting A Hermitian Lamb-shift operator should be added to ˆ H S . Weshow in [SM] that it is negligible for our system. them in Eq. (4) yields an exponential decay of the off-diagonal elements (coherences) in the eigenbasis:˙ ρ nm = − (cid:18) iω nm + 1 T nm (cid:19) ρ nm (10)where iω nm is the dephasing due to the unitary evolu-tion and T nm > Hilbert approximation where the dynamics is pro-jected on the diagonal elements only: it amounts totransit from an eigenstate to the other with rates givenin [SM,Sec.I.E]. The nonsecular dynamics (6) adds tothe right-hand side of Eq. (S48) other entries than ρ nm ,with associated bath rates [SM,Sec.I.F], allowing the ex-istence of coherences in the steady-state solution. TheHilbert approximation is thus retrieved within the non-secular dynamics when the bath timescales are longwith respect to dephasing. Microwave drive – In equilibrium, the steady stateis described by the Boltzmann distribution with eitherchoice of jump operators. The nonsecular evolutionbrings drastic changes out of equilibrium: in a DNP pro-tocol the system is irradiated by microwaves describedby ˆ H MW ( t ) = ω (cid:104) ˆ S x cos( ω MW t ) + ˆ S y sin( ω MW t ) (cid:105) . Thetime dependence is handled using a rotating frameat frequency ω MW where the microwaves are station-ary. The dynamics of the rotated density matrix e iω MW t ˆ S z ρ ( t ) e − iω MW t ˆ S z → ρ ( t ) remains given by Eq. (4)with the shift ˆ H = ˆ H S − ω MW ˆ S z + ω ˆ S x [SM,Sec.I.D]. ω (2 π GHz)FIG. 3. EPR spectra : dots represent numerical pro-files for Hilbert (red) and nonsecular (black) evolutions.Dashed lines are calculated through a spin-temperatureansatz. (top) β − = 1 . β − = 12 K. Thebath timescales are short and the full (nonsecular) dynam-ics gets localized and displays the hole burning. Here theHilbert approximation fails, predicting a spin-temperaturebehavior. Averages are done over 1000 realizations. Numerical computation of the EPR spectrum – Weare interested in comparing the stationary states pre-dicted by the Hilbert dynamics with the ones obtainedby the nonsecular evolution Eq. (4) with jumps (7),(8).In our numerics, the disorder ∆ i is drawn uniformly inthe interval [ − ∆ ω e , ∆ ω e ] and the dipolar couplings U ij are mimicked by independent Gaussian distribu-tions with zero mean and variance U /N , with U ofthe same order of magnitude as ∆ ω e . We fix the disor-der ∆ ω e = 5 · π MHz and U = 0 . · π MHz (note that ω e = 93 . · π GHz) where ˆ H S has ETH statistics. Weconsider two temperatures ; at high temperature thecharacteristic times of the bath are short (see Table I). β − (K) T (0)( µ s) T ∗ ( µ s) T ( µ s) ω (2 π MHz) ω MW (2 π GHz)1.2 1.6 80 160 0.628 93.898812 0 . 16 0 . 16 1 . We compute numerically the steady-state density ma-trix ρ stat [62]. The Hilbert case amounts to a 2 N × N linear system, which for N = 10 spins is sufficientlysmall to be treated by exact diagonalization. The non-secular dynamics Eq. (4) is instead a 4 N × N linearsystem which requires the use of Krylov subspace meth-ods (namely the biconjugate gradient stabilized algo-rithm) [63]. To probe the stationary state we focuson the EPR spectrum. The EPR experiment starts attime τ = 0 by a π/ i on the y axis: ρ π/ = e i π ˆ S xi ρ stat e − i π ˆ S xi . For short times afterthe pulse, the evolution is unitary and the polarizationin the ( x, y ) plane is encoded in g i ( τ ) = − i Tr (cid:104) ˆ S + i ( τ ) ρ π/ (cid:105) (11)where ˆ S + i ( τ ) = e i ˆ H S τ ˆ S + i e − i ˆ H S τ . The EPR spec-trum is then defined by the Fourier transform f ( ω ) = N (cid:80) i Re (cid:104) ´ ∞ τπ g i ( τ ) e − iωτ (cid:105) , explicitly calcu-lated in [SM,Sec.II].The EPR spectra are shown in Fig. 3. At low tem-perature, we observe a spin-temperature curve for bothdynamics. Here the bath timescales are long with re-spect to dephasing. Consequently, the density ma-trix still gets projected in the eigenstate basis, yield-ing a similar situation as in the Hilbert approxima-tion. At higher temperature, the EPR spectrum isspin-temperature-like for Hilbert dynamics, whereas ithas a hole-burning shape in the nonsecular evolution.The spin-temperature behavior observed in the Hilbertapproximation is expected [41–43] : due to the ETH,the jump operators projected on the eigenstates inducechanges in energy and polarization without any otherinformation such as the spatial location of the spins.On the contrary, in the nonsecular equation the jumpoperators are well localized in space and compete with The last dipolar term ∝ ˆ S zi ˆ S zj in Eq. (2) is dropped as itcommutes with all the operators ˆ S zi . n FIG. 4. (main) Entanglement entropy of the eigenstatesof ρ stat (red : Hilbert, black : nonsecular) classified byincreasing energy. To compute it, we performed a partialtrace of each eigenstate | n (cid:105) (cid:104) n | over spins 5 to 10. Verticalgray lines delimit sectors of constant polarization s zn . (inset)Smoothed histogram of S zi expectation values ( i = 1 , . . . , N )for 100 eigenstates in the middle of the s zn = 0 sector. Forthe Hamiltonian eigenstates, the local observable is peakedaround its thermal value ( S zi (cid:39) 0) as required by ETH.The nonsecular eigenstates present instead strong fluctua-tions, as for MBL states. Parameters correspond to those ofFig. 3(bottom). the dephasing, which is not able to project the systemon the eigenstates. The EPR spectrum becomes similarto the MBL case, although the eigenstates of ˆ H S are er-godic for the chosen parameters. The spin-temperaturepicture no longer holds and one can regard this break-down as a type of Zeno localization induced by the bath.This breakdown is confirmed by comparing the EPRprofiles with the ones (dashed lines in Fig. 3) obtainedthrough a spin-temperature ansatz for the steady statedensity matrix ρ ans nn ( β s , h ) ∝ e − β s ( ε n − hs zn ) where thespin temperature β s is conjugated to the energy andthe magnetic field h is conjugated to the polarization( s zn are the eigenvalues of the conserved ˆ S z ). These pa-rameters are computed by a fitting method introducedin [42] and recalled in [SM,Sec.III.B]. Spectral properties – The localization phenomenonexhibited by the experimentally relevant EPR spectracan be revealed through other observables. For in-stance, in [SM,Sec.III] we display the polarization pro-files. In Fig. 4, we focus instead on the spectral prop-erties of ρ stat . We compare the entanglement entropyof the ρ stat eigenstates in the Hilbert and nonsecularcases. The latter case is much less entangled and simi-lar to MBL eigenstates [SM,Sec.III.D]. The present sce-nario is analog to the entanglement transition betweenvolume-law and area-law phases analyzed in Refs. [28–40]. However we stress that the bath-induced local-ization unveiled here is very different from the local-ization studied in Refs. [41–46] and experimentally ob-served in [47]. In that case it is caused by MBL eigen-states that stem from a strong disorder (or weak in-teraction). In [SM,Sec.III.C] we show that in presenceof MBL eigenstates, the nonsecular terms are negligi-ble irrespective of temperature : Hilbert dynamics istherefore sufficiently accurate. Conclusion – We have shown that bath transitionsinduced by localized vibrations compete with the de-phasing (delocalizing long-range interactions) and maybe at the origin of a localization effect in driven sys-tems. To capture this phenomenon it is crucial togo beyond the conventional weak-coupling secular ap-proximation [16] and solve a more general GKSL equa-tion [58–60]. In the context of DNP, our work providesan explanation for the breakdown of thermal mixingupon increasing the temperature of the sample. Its ori-gin lies in the enhanced dynamics of the vibrationalmodes. In our model we used heuristic values of themicroscopic timescales consistent with the experimen-tal values of the standard T and T relaxation times.The present analysis calls for a thorough experimen-tal test of the influence of temperature on the differentregimes of hyperpolarization. Acknowledgments – We thank Fabien Alet, Leti-cia Cugliandolo, Andrea De Luca, Laura Foini, Nico-las Laflorencie, Leonardo Mazza and Marco Schir´o forvaluable discussions around different facets of this work.We also are grateful to Filippo Vicentini for advices onthe numerical methods. We acknowledge funding fromthe ANR-16-CE30-0023-01 (THERMOLOC) grant andthe Galileo Galilei Institute for hospitality during partof this project. Supplementary Material In the supplementary material we provide some technical details and additional numerics. In Sec. I we reviewthe derivation of the different master equations analyzed in the main text (secular, nonsecular and Hilbert masterequations). In Sec. II we derive a numerically convenient formula to compute the EPR spectrum. In Sec. III wedisplay additional plots of EPR and polarization profiles, for an ETH and a MBL Hamiltonian. We explain thefitting procedure introduced in [42] to a spin-temperature ansatz and exhibit extra data concerning the spectralproperties of the bath-induced localized steady state. We also discuss the effect of the different bath timescalesin the DNP model. The numerics are done for the N = 10 spins setup of the main text, but for completenesswe analyze as well a case of full localization with N = 12 spins. Finally in Sec. V we explain the choice of ourDNP Hamiltonian and show a simple model of coupling to the bath that allow to compute the microscopic bathtimescales T ( | ω | ) as a function of frequency and temperature. These timescales are decreasing with temperature,an essential feature for the Zeno localization discussed in the main text. CONTENTS I. Open and driven system: Weak-coupling nonsecular master equation 7A. Born-Markov approximation at weak coupling 7B. Secular approximation and the Gorini-Kossakowski-Sudarshan-Lindblad equation 8C. The nonsecular quantum master equation 91. Derivation of the nonsecular equation 92. Lamb shift and nonsecular jump operators 11D. Turning the microwaves on 11E. The Hilbert approximation 12F. Comparison between secular and nonsecular GKSL equations 14II. Electron Paramagnetic Resonance spectrum 14III. Bath-induced localization by temperature variation: additional numerical results 16A. Thermal mixing and spin temperature behavior 16B. Estimating the spin-temperature ansatz parameters from simulation data 18C. Hole burning 18D. Spectral properties of the stationary state 20IV. A toy model for Zeno localisation 22V. Microscopic model of the electronic spins and bath 23A. Coupling between spins and the magnetic field 23B. Vibrational modes of the embedding material 241. Direct process 242. Two-phonon processes 243. Discussion of the one- and two-phonon bath timescales 25References 26 I. OPEN AND DRIVEN SYSTEM: WEAK-COUPLING NONSECULAR MASTER EQUATION The electronic spin system is in contact with a thermal reservoir and irradiated by microwaves. The totalHamiltonian reads H tot = H S + H MW ( t ) + H int + H B (S1)In this section we will carefully study the system-bath interaction to derive an effective evolution of the spin system.We thus switch the microwaves off H MW ( t ) = 0 and shall reinstate them in Sec. I D. The term H int describes thecoupling of the bath to the spins and H B the Hamiltonian of the bath alone, assumed to remain in equilibrium atinverse temperature β . The system-bath interaction is written as H int = (cid:88) i =1 ,...,Nµ = x,y,z S µi ⊗ B µi (S2)where B µi are Hermitian operators acting on the bath’s Hilbert space, representing vibrational modes of the glassymedium. This choice shall be justified for diluted spins in Sec. V, where it also turns out the bath degrees offreedom can be considered independent : (cid:68) B µi ( t ) B νj ( t (cid:48) ) (cid:69) B = 0 if ( i, µ ) (cid:54) = ( j, ν ) (S3)where (cid:104)•(cid:105) B = Tr( • ρ B ) are bath averages with ρ B = e − βH B / Tr e − βH B , and B µi ( t ) = e itH B B µi e − itH B . Besides thebath degrees of freedom are fluctuations that satisfy (cid:10) B µi (cid:11) B = 0.In the following we first summarize the steps in [16, Sec.3.3] to treat perturbatively the effect of the bath onthe system evolution. We then review the standard secular approximation that allow to turn this expansion intoa Gorini-Kossakowski-Sudarshan-Lindblad master equation. Next we summarize the steps that provide a GKLSform without throwing away nonsecular terms, introduced in Ref. [58], by modifying the Markov approximation.Finally we show how to include the microwave drive in this approach and give the main formal difference betweensecular and nonsecular master equations that account for the effects analyzed in this article. A. Born-Markov approximation at weak coupling We start defining the non-interacting Hamiltonian H = H S + H B and the corresponding unitary evolutionoperator U ( t ) = e − itH . The evolution of the whole system’s density matrix is ruled by the Liouville-von Neumannequation which reads ˙ ρ tot = − i [ H, ρ tot ] in Schr¨odinger’s picture. Since we treat H int as a perturbation it is moreconvenient to start from the interaction picture version of it,˙ ρ ( t ) = − i [ (cid:101) H int ( t ) , ρ ( t )] (S4)where by definition ρ ( t ) = U † ( t ) ρ tot ( t ) U ( t ) and (cid:101) H int ( t ) = U † ( t ) H int U ( t ).We aim at getting an equation for the reduced density matrix describing the spin system ρ S = Tr B ρ . At weakcoupling one (i) solves perturbatively Eq. (S4) by integrating it once, ρ ( t ) = ρ (0) − i ´ t d s [ (cid:101) H int ( s ) , ρ ( s )], andplugging the solution back in Eq. (S4) – note that there is no approximation made in this step – (ii) makes the Born approximation ρ ( t ) (cid:39) ρ S ( t ) ⊗ ρ B . One obtains after tracing over the bath˙ ρ S ( t ) = − ˆ t d s Tr B (cid:104) (cid:101) H int ( t ) , [ (cid:101) H int ( s ) , ρ S ( s ) ⊗ ρ B ] (cid:105) = τ = t − s − ˆ t d τ Tr B (cid:104) (cid:101) H int ( t ) , [ (cid:101) H int ( t − τ ) , ρ S ( t − τ ) ⊗ ρ B ] (cid:105) (S5)The latter equation is not affected by the additional term coming from the initial condition when integratingEq. (S4), as with S µi ( t ) = e itH S S µi e − itH S ,Tr B (cid:104) (cid:101) H int ( t ) , ρ (0) (cid:105) = (cid:88) i,µ (cid:2) S µi ( t ) , ρ S (0) (cid:3) Tr (cid:0) B µi ρ B (cid:1) = 0 (S6)The second approximation, after the weak coupling, is a Markov approximation :˙ ρ S ( t ) = − ˆ ∞ d τ Tr B (cid:104) (cid:101) H int ( t ) , [ (cid:101) H int ( t − τ ) , ρ S ( t ) ⊗ ρ B ] (cid:105) (S7)The latter approximation amounts to say that products like (cid:101) H int ( t ) (cid:101) H int ( t − τ ) decay very rapidly to zero with τ compared to the relaxation time τ R of the system (the relaxation time of ρ S ( t )), which is the case if the bathcorrelation functions (whose relaxation time is τ B ) decay very fast: τ B (cid:28) τ R . The evolution is now Markovian.Note that there is not a unique way to perform this approximation, as there is a wide freedom in the choice of thetime substitution made above t − τ → t : for instance any time t (cid:48) such that | ( t − τ ) − t (cid:48) | (cid:46) τ B could be chosen. It can be shown through the Nakajima-Zwanzig projection op-erator formalism [16, Sec.9.1] that Eq. (S5) is actually exactat second order in H int (provided Eq. (S6)). Alternatively, asimpler cumulant expansion has been devised by Alicki andcollaborators in [56, 57] for the evolution operator Λ such that ρ S ( t ) = Tr B ρ ( t ) = Λ( t, ρ S (0). Taking the time derivative ofthe latter equation, one recovers the same result at second or-der in H int . This result is guessed in the present perturbativeself-consistent derivation. B. Secular approximation and the Gorini-Kossakowski-Sudarshan-Lindblad equation We shall now see that Eq. (S7) can be simplified further using the eigendecomposition of H S . We define theprojector Π( ε ) on the eigenspace associated to the energy ε of H S , and for each gap ω of the spin Hamiltonian S µi ( ω ) ≡ (cid:88) ε (cid:48) − ε = ω Π( ε ) S µi Π( ε (cid:48) ) ⇒ H int = (cid:88) i,µ,ω S µi ( ω ) ⊗ B µi (S10)This allows to simply write (cid:101) H int ( t ) = (cid:88) i,µ,ω e − iωt S µi ( ω ) ⊗ B µi ( t ) = (cid:88) i,µ,ω e iωt S µi ( ω ) † ⊗ B µ † i ( t ) (S11)with S µi ( ω ) † = S µi ( − ω ) from Eq. (S10). Plugging it in Eq. (S7) using the Hermitian conjugated form for (cid:101) H int ( t )and the direct one for (cid:101) H int ( t − τ ), we get˙ ρ S ( t ) = (cid:88) i,jµ,ν (cid:88) ω,ω (cid:48) Γ µνij ( ω ) e i ( ω (cid:48) − ω ) t (cid:16) S νj ( ω ) ρ S ( t ) S µi ( ω (cid:48) ) † − S µi ( ω (cid:48) ) † S νj ( ω ) ρ S ( t ) (cid:17) + H.c. (S12)The equilibrium bath correlation functions Γ µνij ( ω ) satisfy time-translational invariance owing to the commutation[ H B , ρ B ] = 0, so that using (cid:68) B µi ( t ) B νj ( t − τ ) (cid:69) = (cid:68) B µi ( τ ) B νj (0) (cid:69) , they are defined asΓ µνij ( ω ) = ˆ ∞ d τ e iωτ (cid:68) B µi ( τ ) B νj (0) (cid:69) B = δ ij δ µν Γ( ω ) (S13)We assumed as previously mentioned that different degrees of freedom of the bath decouple. Let us then writeEq.(S12) back to the Schr¨odinger picture (see Eq. (S8)):˙ ρ s S ( t ) = − i (cid:2) H S , ρ s S ( t ) (cid:3) + (cid:88) i,µ (cid:88) ω,ω (cid:48) (cid:26) Γ( ω ) (cid:16) S µi ( ω ) ρ s S ( t ) S µi ( ω (cid:48) ) † − S µi ( ω (cid:48) ) † S µi ( ω ) ρ s S ( t ) (cid:17) + H.c. (cid:27) (S14)where we note that the phases disappear compared to the Heisenberg picture equation (S12). Eq.(S14) is writtenusing an eigenbasis of H S . This is a convenient choice, but we could have chosen a different basis, yielding a morecomplicated equation, and at this level of approximation the time evolution of ρ s S would be identical.Nonetheless our Markovian quantum master equation must preserve the defining properties of a density matrix(such as positivity and trace), i.e. be of the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) form ˙ ρ = L ρ wherethe superoperator L is a generator of a quantum dynamical semi-group [16, Sec.3.2]. A more explicit definitionis given by Eq. (4) of the main text. This is not guaranteed by Eq. (S14). To ensure so in this weak-couplingapproach, the standard prediction resorts to the secular approximation . It amounts to neglect the terms for which ω (cid:54) = ω (cid:48) . It is then useful to define the real and imaginary parts:Γ( ω ) = γ ( ω )2 + iS ( ω ) , S ( ω ) = Im Γ( ω ) γ ( ω ) =2Re Γ( ω ) = ˆ ∞−∞ d τ e iωτ (cid:10) B µi ( τ ) B µi (0) (cid:11) B (S15)As ρ B describes the Boltzmann-Gibbs distribution at inverse temperature β , one can easily prove the Kubo-Martin-Schwinger condition (cid:68) B µ † i ( t ) B νj (0) (cid:69) B = (cid:68) B νj (0) B µ † i ( t + iβ ) (cid:69) B , from which the detailed balance condition γ ( ω ) /γ ( − ω ) = e βω follows [64, Sec.12]. We can thus equivalently define γ ( ω ) = h ( ω ) T ( | ω | ) with h ( ω ) = e βω e βω (S16) Let us note ρ s S = Tr B ρ tot the Schr¨odinger picture density ma-trix. We have ρ S ( t ) =Tr B (cid:16) U † ( t ) ρ tot ( t ) U ( t ) (cid:17) = e itH S Tr B (cid:16) e itH B ρ tot ( t ) e − itH B (cid:17) e − itH S = e itH S ρ s S ( t ) e − itH S (S8) i.e. the standard relationship between Heisenberg andSchr¨odinger pictures for the reduced system S. If we switchto the Schr¨odinger picture, Eq. (S7) becomes ˙ ρ s S ( t ) = − i (cid:2) H S , ρ s S ( t ) (cid:3) + (cid:88) i,jµ,ν ˆ ∞ d τ (cid:16) S νj ( − τ ) ρ s S ( t ) S µi − S µi S νj ( − τ ) ρ s S ( t ) (cid:17) (cid:68) B µi ( τ ) B νj (cid:69) B + H.c. (S9)We used time-translation invariance of the equilibrium bathcorrelation function. One sees that the extra difficulty of evolv-ing in time the interaction Hamiltonian remains even in thisformulation, which calls for the eigendecomposition of H S . T ( | ω | ) defines relaxation times of the system for transitions at energy ω while h enforces the detailed balancecondition. Inserting the above definitions in Eq. (S14) where the ω (cid:54) = ω (cid:48) terms are neglected, one arrives at aGKSL master equation: ˙ ρ s S ( t ) = − i (cid:2) H S + H LS , ρ s S ( t ) (cid:3) + D ρ s S ( t ) H LS = (cid:88) i,µ,ω S ( ω ) S µ † i ( ω ) S µi ( ω ) D ρ = (cid:88) i,µ,ω h ( ω ) T ( | ω | ) (cid:18) S µi ( ω ) ρS µ † i ( ω ) − (cid:110) S µ † i ( ω ) S µi ( ω ) , ρ (cid:111)(cid:19) (S17)The Lamb-shift Hamiltonian H LS being Hermitian and [ H S , H LS ] = 0, it can be absorbed into the system Hamil-tonian H S in Eq. (S17) just by shifting the energy levels.For our many-body electron system described by the Hamiltonian H S , it is safe to consider the energy gaps asnon degenerate. Eq. (S17) takes then a simple form in the eigenbasis {| n (cid:105)} ˙ ρ nn = (cid:88) m (cid:54) = n W B mn ρ mm − W B nm ρ nn ˙ ρ m (cid:54) = n = − (cid:18) iω mn + 1 T mn (cid:19) ρ m (cid:54) = n (S18)with respectively the bath and decoherence rates W B nm = h ( ω nm ) T ( | ω nm | ) (cid:88) i,µ (cid:12)(cid:12) (cid:104) n | S µi | m (cid:105) (cid:12)(cid:12) T nm = 12 (cid:88) k (cid:54) = n h ( ω nk ) W B nk + 12 (cid:88) k (cid:54) = m h ( ω mk ) W B mk + 14 T (0) (cid:88) i (cid:104) n | S zi | n (cid:105) + (cid:104) m | S zi | m (cid:105) − (cid:104) n | S zi | n (cid:105) (cid:104) m | S zi | m (cid:105) (S19)The first line in Eq. (S18) is a master equation for the diagonal elements only, with rates satisfying the detailedbalance condition W B mn /W B nm = e βω mn , and ω mn = ε m − ε n the gap between levels m and n of the Hamiltonian H S + H LS . The second equation describes the evolution of the off-diagonal terms: the first term iω mn represents thedephasing from the unitary evolution while the second term T − mn describes the decoherence due to the bath. Underthis evolution, the coherences vanish exponentially while the populations tend to their thermal values imposed bythe Boltzmann distribution, emerging from the detailed balance at temperature β − .Note that after the secular approximation, the time evolution of ρ s S depends on the choice of the eigenbasis. Inpresence of a drive, different steady states can be reached by changing this choice of basis, which is an inconsistencyof the secular approximation. C. The nonsecular quantum master equation A usual justification [16] for the secular approximation, on top of allowing to preserve the properties of a densitymatrix, is that in Heisenberg representation the terms ω (cid:54) = ω (cid:48) produce rapidly oscillating phases (see Eq. (S12))which cancel. Yet for a large many-body system there exist a lot of quasi-degenerate gaps for which the argumentdoes not hold. There is a way to retain all the terms from Eq. (S12) while ensuring a GKSL form by modifying theMarkovian approximation. Below we summarize the derivation as given by Nathan and Rudner’s recent paper [58]. 1. Derivation of the nonsecular equation We start back from Eqs. (S5)-(S7). Coming back to the variable s we see that the initial condition is in the endpreponed to −∞ by the Markov approximation performed there, and we have˙ ρ S ( t ) = − ˆ t −∞ d s Tr B (cid:104) (cid:101) H int ( t ) , [ (cid:101) H int ( s ) , ρ S ( t ) ⊗ ρ B ] (cid:105) = − (cid:88) i,µ ˆ t −∞ d t (cid:48) Γ( t − t (cid:48) ) (cid:2) S µi ( t ) , S µi ( t (cid:48) ) ρ S ( t ) (cid:3) + H.c. (S20)A first progress is to take the square root of the bath rate in Fourier space, i.e. to define g such thatΓ( t − t (cid:48) ) = ˆ R d s g ( t − s ) g ( s − t (cid:48) ) (S21)Eq. (S20) becomes ˙ ρ S ( t ) = ˆ R d t (cid:48) d s F ( t, s, t (cid:48) )[ ρ S ( t )] F ( t, s, t (cid:48) )[ ρ ] = (cid:88) i,µ θ ( t − t (cid:48) ) g ( t − s ) g ( s − t (cid:48) ) (cid:2) S µi ( t ) , S µi ( t (cid:48) ) ρ (cid:3) + H.c. (S22)0Integrating on [ t , t ] with t − t (cid:29) τ B : ρ S ( t ) − ρ S ( t ) = ˆ t t d t ˆ R d t (cid:48) d s F ( t, s, t (cid:48) )[ ρ S ( t )] (S23)Next, due to t − t (cid:29) τ B , most contributions in the integrals of Eq. (S23) arise for ( t, s, t (cid:48) ) ∈ [ t , t ] , thus they areapproximately unaffected if we change the boundaries as {−∞ < s, t (cid:48) < ∞ , t (cid:54) t (cid:54) t } → {−∞ < t, t (cid:48) < ∞ , t (cid:54) s (cid:54) t } : ρ S ( t ) − ρ S ( t ) (cid:39) ˆ t t d s ˆ R d t d t (cid:48) F ( t, s, t (cid:48) )[ ρ S ( s )] (S24)where, similarly to Eq. (S7) one makes another Markovian type of approximation for the density matrix timedependence. Then we derive Eq. (S24) with respect to t and set t = t :˙ ρ S ( t ) = L ( t )[ ρ S ( t )] , L ( t ) = ˆ R d s d s (cid:48) F ( s, t, s (cid:48) ) (S25)and using θ ( t ) = 1 / t ) / ρ S ( t ) = − i (cid:2) H LS ( t ) , ρ S ( t ) (cid:3) + (cid:88) i,µ A µi ( t ) ρ S ( t ) A µi ( t ) † − (cid:110) A µi ( t ) † A µi ( t ) , ρ S ( t ) (cid:111) A µi ( t ) = ˆ R d s g ( t − s ) S µi ( s ) , H LS ( t ) = 12 i (cid:88) i,µ ˆ R d s d s (cid:48) S µi ( s ) g ( t − s ) g ( t − s (cid:48) ) S µi ( s (cid:48) )sign( s − s (cid:48) ) (S26)By Fourier transform and the definition (S21) we can relate g ( t ) and γ ( ω ):˜ g ( ω ) = ˆ R d t π e iωt g ( t ) ⇒ ˜ g ( ω ) = (cid:112) γ ( ω )2 π (S27)We now seek to come back to Schr¨odinger representation and decompose the operators on the eigenbasis of thespin Hamiltonian through S µi ( t ) = (cid:80) ω S µi ( ω ) e − iωt . For the Lamb-shift Hamiltonian we show through the Fourierdecompositions H LS ( t ) = (cid:88) i,µ,ω,ω (cid:48) S µi ( ω ) S µi ( ω (cid:48) ) e − it ( ω + ω (cid:48) ) ˆ R dΩdΩ (cid:48) ˜ g (Ω)˜ g (Ω (cid:48) ) k ( ω + Ω , ω (cid:48) − Ω (cid:48) ) k ( p, q ) = 12 i ˆ R d s d s (cid:48) sign( s − s (cid:48) ) e − i ( ps + qs (cid:48) ) (S28)Through a direct computation and usual regularizations we get k ( p, q ) = 2 πδ ( p + q )Im 1 ip + 0 + (S29)which yields H LS ( t ) = (cid:88) i,µ,ω,ω (cid:48) S µi ( ω ) S µi ( ω (cid:48) ) e − it ( ω + ω (cid:48) ) π P ˆ R dΩ ˜ g (Ω − ω )˜ g (Ω + ω (cid:48) )Ω (S30)where P stands for Cauchy’s principal value. It is then straightforward to pass from Heisenberg to Schr¨odingerrepresentations, which cancels the phases. This allows to define nonsecular jump operators A µi = (cid:88) ω (cid:112) γ ( ω ) S µi ( ω ) (S31)In Schr¨odinger representation, the nonsecular version of Eq. (S17) thus reads:˙ ρ s S ( t ) = − i (cid:2) H S + H LS , ρ s S ( t ) (cid:3) + (cid:88) i,µ A µi ρ s S ( t ) A µ † i − (cid:110) A µ † i A µi , ρ s S ( t ) (cid:111) H LS = (cid:88) i,µ,ω,ω (cid:48) S µi ( ω ) S µi ( ω (cid:48) ) f ( ω, ω (cid:48) ) , f ( ω, ω (cid:48) ) = P ˆ R dΩ2 π (cid:112) γ (Ω − ω ) γ (Ω + ω (cid:48) )Ω (S32)The Lamb-shift H LS is Hermitian and the evolution is thus manifestly in GKSL form. The secular equation isrecovered again by throwing away all gaps ω (cid:54) = ω (cid:48) in the sums.1 2. Lamb shift and nonsecular jump operators For our system the Lamb-shift H LS becomes negligible due to the wide difference of characteristic energy scales.Indeed considering that f ( ω, ω (cid:48) ) is smooth in both arguments, f ( ω, ω (cid:48) ) (cid:39) f (0 , 0) for S zi gaps while for S x,yi we canuse S ± i = S xi ± iS yi operators and get: H LS (cid:39) f (0 , (cid:88) i ( S zi ) + (cid:88) i,ω,ω (cid:48) (cid:16) S + i ( ω ) + S − i ( ω ) (cid:17) (cid:16) S + i ( ω (cid:48) ) + S − i ( ω (cid:48) ) (cid:17) f ( ω, ω (cid:48) ) − (cid:88) i,ω,ω (cid:48) (cid:16) S − i ( ω ) − S + i ( ω ) (cid:17) (cid:16) S − i ( ω (cid:48) ) − S + i ( ω (cid:48) ) (cid:17) f ( ω, ω (cid:48) )= N f (0 , 0) + 12 (cid:88) i,ω,ω (cid:48) f ( ω, ω (cid:48) ) (cid:104) S + i ( ω ) S − i ( ω (cid:48) ) + S − i ( ω ) S + i ( ω (cid:48) ) (cid:105) (cid:39) N f (0 , 0) + 12 (cid:88) i,ω,ω (cid:48) (cid:104) f ( − ω e , ω e ) S + i ( ω ) S − i ( ω (cid:48) ) + f ( ω e , − ω e ) S − i ( ω ) S + i ( ω (cid:48) ) (cid:105) = N f (0 , 0) + f ( ω e , − ω e )2 (cid:88) i (cid:16) S + i S − i + S − i S + i (cid:17) = N f (0 , 0) + 2 f ( ω e , − ω e )) (S33)where we have used the symmetry f ( ω e , − ω e ) = f ( − ω e , ω e ). In conclusion the Lamb-shift Hamiltonian is propor-tional to the identity up to very small corrections and thus we can take (cid:2) H LS , ρ s S ( t ) (cid:3) = 0 in Eq. (S32).With very similar considerations one obtains from the definition (S31) the form of the nonsecular jump operatorsgiven by Eq. (7) of the main text. D. Turning the microwaves on In the following we set H LS = 0 and incorporate the microwave field H MW ( t ) = ω (cid:2) S x cos( ω MW t ) + S y sin( ω MW t ) (cid:3) (S34)The interaction picture dynamics Eq. (S4) becomes˙ ρ ( t ) = − i [ (cid:101) H MW ( t ) , ρ ( t )] − i [ (cid:101) H int ( t ) , ρ ( t )] (S35)with (cid:101) H MW ( t ) = U † ( t ) H MW ( t ) U ( t ) = e iH S t H MW ( t ) e − iH S t . Then we treat the last term as before by integratingonce the equation and obtain˙ ρ ( t ) = − i [ (cid:101) H MW ( t ) , ρ ( t )] − ˆ t d s (cid:104) (cid:101) H int ( t ) , [ (cid:101) H int ( s ) , ρ ( s )] (cid:105) − ˆ t d s (cid:104) (cid:101) H int ( t ) , [ (cid:101) H MW ( s ) , ρ ( s )] (cid:105) − i [ (cid:101) H int ( t ) , ρ (0)] (S36)Then if we make the Born approximation ρ ( t ) (cid:39) ρ S ( t ) ⊗ ρ B and trace over the bath, since (cid:101) H MW ( s ) acts only onthe spin degrees of freedom, the last two terms in Eq. (S36) are zero if the bath degrees of freedom are of averagezero, (cid:10) B µi ( t ) (cid:11) B = 0. Then the quadratic term in H int can be dealt with as in the previous sections.In Schr¨odinger picture we therefore get the equation, similar to Eq. (S32):˙ ρ s S ( t ) = − i (cid:2) H MW ( t ) , ρ s S ( t ) (cid:3) − i (cid:2) H S , ρ s S ( t ) (cid:3) + D ρ s S ( t ) D ρ = − i [ H LS , ρ ] + (cid:88) i,µ A µi ρA µ † i − (cid:110) A µ † i A µi , ρ (cid:111) (S37)Note that H MW ( t ) is time dependent ; a usual trick for NMR studies is to place ourselves in a rotated frame given bythe Larmor precession [65]. We shall now rotate accordingly the reference frame, which, if done at frequency ω MW ,makes the microwaves field effectively stationary, resulting in a time-independent dynamical semigroup generatorfor Eq. (S37). The rotation is implemented as U MW ( t ) = e − iω MW tS z , ρ r ( t ) = U † MW ( t ) ρ s S ( t ) U MW ( t ) (S38)Since [ S z , H S ] = 0, the unperturbed spin Hamiltonian is unchanged. We then need to express the rotated microwaveHamiltonian U † MW ( t ) H MW ( t ) U MW ( t ). To achieve this one can write the microwave Hamiltonian as H MW ( t ) = ω (cid:16) S − e iω MW t + S + e − iω MW t (cid:17) (S39) If the field is oscillating along a single direction, e.g. H MW ( t ) = ω S x cos( ω MW t ), the conclusion of this section re- mains valid at small drive through the rotating wave approxi-mation [41, 42]. S ± ( t ) = U † MW ( t ) S ± U MW ( t ) verifying˙ S ± ( t ) = U † MW ( t ) iω MW (cid:2) S z , S ± (cid:3) U MW ( t ) = ± iω MW S ± ( t ) ⇒ S ± ( t ) = e ± iω MW t S ± (S40)where we used (cid:2) S z , S ± (cid:3) = ± S ± . We thus conclude that in the rotating frame the microwave Hamiltonian becomesstationary: U † MW ( t ) H MW ( t ) U MW ( t ) = ω S x (S41)and the GKSL equation (S37) becomes˙ ρ r ( t ) = − i (cid:2) H S + ω S x − ω MW S z , ρ r ( t ) (cid:3) + U † MW ( t ) D (cid:16) U MW ( t ) ρ r ( t ) U † MW ( t ) (cid:17) U MW ( t ) (S42)Let us define as before H S | n (cid:105) = ε n | n (cid:105) the eigendecomposition. As [ S z , H S ] = 0, the total polarization along z isconserved, and we note S z | n (cid:105) = s zn | n (cid:105) . We consider as well for large enough N that the finite gaps ω nm = ε n − ε m are non degenerate, implying that the jump operators read: S zi (0) = (cid:88) n (cid:104) n | S zi | n (cid:105) | n (cid:105)(cid:104) n | , S µi ( ω nm ) = | m (cid:105)(cid:104) m | S µi | n (cid:105)(cid:104) n | (S43)To analyze the last term of Eq. (S42), let us go back to the original dissipative part of Eq. (S37) (second line). Itconsists in sums over pairs of gaps ( ω, ω (cid:48) ) of the jumps operators (S43). Each term in this sum contains a productof the form (cid:104) n | S µi ( ω = ω mn ) | m (cid:105) (cid:104) m (cid:48) | S µi ( ω (cid:48) = ω m (cid:48) n (cid:48) ) † | n (cid:48) (cid:105) (S44)The bath operators are of two kinds: S xi , S yi which induce flips of the spin i with energy jump ≈ ± ω e and S zi whosetransitions have gaps ω (cid:28) ω e which conserve the total polarization along z . For the S xi , S yi operators, there are 4possibilities: • ω ≈ ω e and ω (cid:48) ≈ − ω e : both factors imply a flip |−(cid:105) i → | + (cid:105) i , so that i (cid:104) n | S xi | m (cid:105) = (cid:104) n | S yi | m (cid:105) and i (cid:104) m (cid:48) | S xi | n (cid:48) (cid:105) = (cid:104) m (cid:48) | S yi | n (cid:48) (cid:105) . Consequently the terms generated by the operator S yi have the opposite value tothe one generated by S xi : and their contribution vanishes. • ω ≈ − ω e and ω (cid:48) ≈ ω e : both factors imply a flip | + (cid:105) i → |−(cid:105) i , so that − i (cid:104) n | S xi | m (cid:105) = (cid:104) n | S yi | m (cid:105) and − i (cid:104) m (cid:48) | S xi | n (cid:48) (cid:105) = (cid:104) m (cid:48) | S yi | n (cid:48) (cid:105) . Once again their combined contribution vanishes. • ω ≈ ω (cid:48) ( ≈ ± ω e ): contrary to the above cases those terms are allowed.For all allowed transitions the total polarization jump for the term (S44) ( i.e. from | n (cid:105) to | n (cid:105) (cid:48) ) is s zn − s zm − ( s zn (cid:48) − s zm (cid:48) ) = 0. This holds for the S xi , S yi operators (for which s zn − s zm = s zn (cid:48) − s zm (cid:48) = ± S zi transitionswhere s zn = s zm and s zn (cid:48) = s zm (cid:48) .When shifting to the rotating frame, the rotation (S38) produces oscillating phases in the dissipative term ofEq. (S42) according to the following relation U MW ( t ) S µi ( ω mn ) U † MW ( t ) = e iω MW t ( s zn − s zm ) S µi ( ω mn ) (S45)Therefore each term writing as Eq.(S44) receives a phase e iω MW t [ s zn − s zm − ( s zn (cid:48) − s zm (cid:48) ) ]. As emphasized above, this phaseis 1 for all allowed transitions. Consequently U † MW ( t ) D (cid:16) U MW ( t ) ρ r ( t ) U † MW ( t ) (cid:17) U MW ( t ) = D ρ r ( t ), i.e. the quantummaster equation in the rotating frame is˙ ρ r ( t ) = − i (cid:2) H S + H LS + ω S x − ω MW S z , ρ r ( t ) (cid:3) + (cid:88) i,µ A µi ρ r ( t ) A µ † i − (cid:110) A µ † i A µi , ρ r ( t ) (cid:111) (S46)Note that if for example only the spin j is irradiated, one has to replace S x → S xj . E. The Hilbert approximation In this section we provide the details of the Hilbert approximation for the secular GKSL equation given byEq. (S46) where nonsecular terms ω (cid:54) = ω (cid:48) are thrown away (as in Eq. (S17)), which amounts to take A µi → (cid:112) γ ( ω ) S µi ( ω ) in Eq. (S46). Below we drop the rotating frame index r and consider H LS = 0 as explained in theprevious sections:˙ ρ = − i [ H S − ω MW S z + ω S x , ρ ] + (cid:88) ω,i,µ h ( ω ) T ( | ω | ) (cid:18) S µi ( ω ) ρS µi ( ω ) † − (cid:110) S µi ( ω ) † S µi ( ω ) , ρ (cid:111)(cid:19) (S47)3 FIG. S5. Effective quantum dynamics of the electron spins. (a) Energy levels of H S are drawn. They are arranged into N + 1 polarization sectors, each containing many levels, separated by multiples of ω e . In the Hilbert approximation due tofast dephasing, the jump operators (arrows) select transitions between eigenstates with corresponding timescales, and thesystem is a mixture of H S eigenstates. (b) Example of a nonsecular quantum trajectory in its average energy vs time plane:the unitary dynamics (wavy arrows) gets projected with a characteristic rate by spatially-localized jump operators. TheHilbert approximation is retrieved for slow jump rates for which the unitary evolution is able to project the system ontoeigenstates. Projecting on the diagonal and off-diagonal terms ( n (cid:54) = m ) leads to˙ ρ nn = (cid:88) k (cid:54) = n W B kn ρ kk − W B nk ρ nn − iω (cid:0) (cid:104) n | S x | k (cid:105) ρ kn − c.c. (cid:1) ˙ ρ nm = − (cid:18) i ∆ ω nm + 1 T nm (cid:19) ρ nm − iω (cid:88) k (cid:0) (cid:104) n | S x | k (cid:105) ρ km − (cid:104) k | S x | m (cid:105) ρ nk (cid:1) (S48)where rates have been defined in Eq. (S19) and the dephasing is now∆ ω nm = ε n − ε m − ( s zn − s zm ) ω MW (S49)If the term i ∆ ω nm + 1 /T nm dominates the coherences’ evolution – making them vanish exponentially fast – thedynamics is projected on the diagonal elements. This is the spirit of the Hilbert approximation . It is achieved byconsidering T nm and 1 / ∆ ω nm go to zero with ∆ ω nm T nm constant and solving perturbatively the master equation.We thus write the evolution as ˙ ρ = L ρ + L ρ with L the dominant contribution to the generator, defined by L e nm = n = m − (cid:18) i ∆ ω nm + 1 T nm (cid:19) e nm n (cid:54) = m (S50)with { e nm } a basis for our 2 N × N density matrices: ( e nm ) ab = δ an δ bm . At dominant order the steady-state densitymatrices fall in the subspace of diagonal matrices. Treating L as a perturbation, the Hilbert approximation projectsthe dynamical evolution in the diagonal subspace, leading to an approximate master equation for the populations.Details of the (Schrieffer-Wolff) perturbation theory are described in [41, 66, 67]. One has to go to second order in L to get the first terms dependent on ∆ ω nm and T nm . One finds the master equation for the populations˙ ρ nn = (cid:88) n (cid:48) ( L ) nn,n (cid:48) n (cid:48) + (cid:88) m (cid:54) = m (cid:48) ( L ) nn,mm (cid:48) ( L ) mm (cid:48) ,n (cid:48) n (cid:48) i ∆ ω mm (cid:48) + 1 /T mm (cid:48) ρ n (cid:48) n (cid:48) (S51)with the basis decomposition ( L ρ ) nm = (cid:80) ij ( L ) nm,ij ρ ij . Concretely the Hilbert rate equation (S51) reads in thepresent secular case ˙ ρ nn = (cid:88) m (cid:54) = n (cid:104) W B mn + W MW nm (cid:105) ρ mm − (cid:104) W B nm + W MW nm (cid:105) ρ nn (S52)where the microwave rates are W MW nm = 2 ω T nm T nm ∆ ω nm ) (cid:12)(cid:12) (cid:104) n | S x | m (cid:105) (cid:12)(cid:12) (S53)The physical interpretation of Hilbert versus nonsecular dynamics is sketched in Fig. S5.4 F. Comparison between secular and nonsecular GKSL equations In this section we highlight the additional terms contained in the nonsecular GKSL equation (S46) comparedto the secular one (S47) (with H LS = 0). We write the off-diagonal elements of the nonsecular equation in theeigenbasis, as these contain the main new terms to account for the bath-induced localization. The first line of eachequation corresponds to the secular part given in Eq. (S48).˙ ρ nm = − (cid:18) i ∆ ω nm + 1 T nm (cid:19) ρ nm − iω (cid:88) k (cid:0) (cid:104) n | S x | k (cid:105) ρ km − (cid:104) k | S x | m (cid:105) ρ nk (cid:1) + (cid:88) iµ = x,y,z (cid:88) k (cid:54) = nk (cid:48) (cid:54) = m (cid:112) γ ( ω kn ) γ ( ω k (cid:48) m ) (cid:104) n | S µi | k (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) ω = ω kn ρ kk (cid:48) (cid:104) k (cid:48) | S µi | m (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) ω (cid:48) = ω k (cid:48) m − (cid:88) kk (cid:48) (cid:54) = n (cid:112) γ ( ω nk ) γ ( ω k (cid:48) k ) (cid:104) n | S µi | k (cid:105)(cid:104) k | S µi | k (cid:48) (cid:105) ρ k (cid:48) m (cid:124) (cid:123)(cid:122) (cid:125) ω (cid:48) = ω nk , ω = ω k (cid:48) k − (cid:88) kk (cid:48) (cid:54) = m (cid:112) γ ( ω k (cid:48) k ) γ ( ω mk ) ρ nk (cid:48) (cid:104) k (cid:48) | S µi | k (cid:105)(cid:104) k | S µi | m (cid:105) (cid:124) (cid:123)(cid:122) (cid:125) ω (cid:48) = ω k (cid:48) k , ω = ω mk (S54)Some remarks:1. The constraints over eigenstates indices below the sum signs correspond to secular terms already taken intoaccount in the first line of the equation. For instance if the constraint is k (cid:54) = n this means that the sameterm with k = n is a secular term contained in the first line.2. For S x,yi terms, we have explicitly written the pair of gaps ( ω, ω (cid:48) ) involved. This emphasizes, as remarked inSec. I D, that terms not satisfying ω ≈ ω (cid:48) actually vanish, although we have not mentioned it in the sums tolighten the notation. This is not true for S zi terms.3. In the secular equation at small drive allows only a steady-state solution where all coherences ρ nm = 0.The nonsecular terms render possible non-zero coherences in the stationary state: the system does not getprojected anymore onto the Hamiltonian eigenstates. This can occur if timescale T ∗ gets lowered so that thenonsecular terms in Eq. (S54) do have an impact. This timescale comes from S zi transitions. All nonsecular S zi transitions in Eq. (S54) occur between different eigenstates ; as a consequence, in the limit of strongdisorder or non-interacting spins, such terms vanish. This explains why the drastic effect of the bath seen foran ETH Hamiltonian is not present in the case of a MBL Hamiltonian, see Fig S7. II. ELECTRON PARAMAGNETIC RESONANCE SPECTRUM The aim of this section is to show how we compute numerically the EPR spectra. For this we need to extend theresults of [68] to the nonsecular case in which the stationary density matrix ρ stat possesses non-zero coherences inthe eigenbasis of the system Hamiltonian. This provides the formula (S63) for a smoothed EPR spectrum which isappropriate for numerical evaluation.We recall that after the π/ τ = 0) the densitymatrix gets transformed into ρ π/ = e i π S xi ρ stat e − i π S xi . The polarization of spin i is given at later times by P µi ( τ ) = 2Tr( S µi ( τ ) ρ π/ ) , S µi ( τ ) = e iH S τ S µi e − iH S τ (S55)and we recall the quantity defined in the main text: g i ( τ ) = ( P yi − iP xi )( τ ) = − i Tr (cid:104) S + i ( τ ) ρ π/ (cid:105) (S56)Note that in our formalism the steady-state density matrix ρ that we get numerically is expressed in the ( x, y )plane in the rotating frame at frequency ω MW , whereas ρ stat is in the fixed frame. Their relationship is ρ stat = U MW ( τ ) ρU MW ( τ ) † .The polarization on the z axis (which is an invariant axis) in the steady-state can be computed with eitherdensity matrix. Indeed,Tr ( S zi ρ stat ) = (cid:88) n,m e iω MW ( s zm − s zn ) τ ρ nm (cid:104) m | S zi | n (cid:105) = (cid:88) n,ms zn = s zm ρ nm (cid:104) m | S zi | n (cid:105) = Tr ( S zi ρ ) (S57)as (cid:104) m | S zi | n (cid:105) vanishes outside the blocks of constant S z (defined by s zm = s zn ).In the following we focus on the EPR spectrum instead. First, to express the post-pulse polarization with ρ instead of ρ π/ , we write e i π S xi = cos π iS xi sin π √ + 2 iS xi ) (S58)5Starting from the definition Eq. (S56) : g i ( τ ) = − i (cid:88) k,l,m,n e iτ [ ε k − ε l + ω MW ( s zm − s zn )] (cid:104) m | − iS xi | k (cid:105)(cid:104) k | S + i | l (cid:105)(cid:104) l | iS xi | n (cid:105) ρ nm (S59)We thus expect the frequency spectrum to be a sum of exponentially many Dirac delta peaks. In order to smooththe spectrum we shall average over a small frequency window δω , i.e. we compute the EPR spectrum as˜ f ( ω ) = 1 N N (cid:88) i =1 δω ˆ ω + δωω f i , f i ( ω ) = Re (cid:20) ˆ ∞ d τπ g i ( τ ) e − iωτ − ητ (cid:21) (S60)where we introduce a small cutoff η > 0. In essence we have to regularize integrals using a Cauchy principal valueover finite intervals of length δω . We indeed note that ˆ ∞ d τπ e − i ( ω − ω ) τ − ητ = 1 π i ( ω − ω ) + η = η/π ( ω − ω ) + η − iπ ω − ω ( ω − ω ) + η (S61)The real part is a Lorentzian and will act as δ ( ω − ω ) for η → + while the other (principal value) imaginary partremains integrable around ω = ω ∀ η > ∝ η − ).Performing the Fourier transforms yields˜ f ( ω ) = 1 N δω (cid:88) k,l,m,n Re (cid:34) ( ω klmn ∈ I ω ) − iπ ln (cid:12)(cid:12)(cid:12)(cid:12) ω + δω − ω klmn ω − ω klmn (cid:12)(cid:12)(cid:12)(cid:12)(cid:35) × ρ nm (cid:88) i (cid:104) m | − iS xi | k (cid:105)(cid:104) k | S + i | l (cid:105)(cid:104) l | iS xi | n (cid:105) (S62)where I ω = [ ω, ω + δω ] and ω klmn = ε k − ε l + ω MW ( s zm − s zn ). The index function is 1 if true and 0 else.The interval of EPR frequencies experimentally probed is centered around ω e . So, as the energy scales betweendifferent polarization sectors are well separated, the terms between brackets in the first line are negligible unless s zk − s zl + s zm − s zn = 1. Yet the factor (cid:104) k | S + i | l (cid:105) implies s zk − s zl = 1, thus considering only matrix elements suchthat s zm = s zn is enough to compute the EPR spectrum.This means in the second line of Eq. (S62) the polarization from | n (cid:105) to | m (cid:105) must not change. As S + i performspolarization jumps of +1, we can get rid of the non-polarization-conserving terms and get, relabeling indices,˜ f ( ω ) = 1 N δω (cid:88) n,k Re (cid:34) ( ω nk ∈ I ω ) − iπ ln (cid:12)(cid:12)(cid:12)(cid:12) ω + δω − ω nk ω − ω nk (cid:12)(cid:12)(cid:12)(cid:12)(cid:35) × (cid:88) ms zm = s zn ρ mn (cid:88) i (cid:104) n | S + i | k (cid:105)(cid:104) k | S − i | m (cid:105) − (cid:88) ms zm = s zk ρ km (cid:88) i (cid:104) m | S − i | n (cid:105)(cid:104) n | S + i | k (cid:105) (S63)Note that if the density matrix is diagonal in the eigenbasis, this implies˜ f ( ω ) = 1 N δω (cid:88) n,k ( ω nk ∈ I ω )( ρ nn − ρ kk ) (cid:88) i (cid:12)(cid:12)(cid:12) (cid:104) n | S + i | k (cid:105) (cid:12)(cid:12)(cid:12) (S64)which agrees with the result of [68, Eq.(23)]. In this case we have, noting f ( ω ) = (cid:80) i f i ( ω ) /N (for η → + ): f ( ω ) = 1 N (cid:88) i (cid:88) n,m δ ( ω − ω nm )( ρ nn − ρ mm ) (cid:12)(cid:12)(cid:12) (cid:104) n | S + i | m (cid:105) (cid:12)(cid:12)(cid:12) (S65)and only the gaps ω nm ≈ ω e with s zn = s zm + 1 actually contribute to the sum. So writing this set as anapproximate interval (for large N ) I = { ω nm with s zn = s zm + 1 } where the EPR spectrum is non zero, one gets,using (cid:104) S + i , S − i (cid:105) = 2 S zi , ˆ I d ω f ( ω ) = 1 N (cid:88) i (cid:34)(cid:88) n ρ nn (cid:104) n | S + i S − i | n (cid:105) − (cid:88) m ρ mm (cid:104) m | S − i S + i | m (cid:105) (cid:35) = 2 N (cid:88) i Tr( S zi ρ ) = 2 (cid:104) S z (cid:105) N = P z (S66) i.e. the integral of the EPR spectrum gives the total polarization along z , P z .Finally, for the more trivial case of non-interacting spins with Hamiltonian H S = (cid:80) i ( ω e + ∆ i ) S zi = (cid:80) i ω i S zi ,the previous integral holds as it is a particular case, but more precisely in Eq. (S65) at fixed spin i the transition Note however that there the gaps are largely degenerated and one must take this into account in our dynamical equations. (cid:104) n | S + i | m (cid:105) is non-zero iff ω nm = ω i , so that we can factor out the Dirac delta function and get f ( ω ) = 1 N (cid:88) i δ ( ω − ω i ) (cid:88) n,m ( ρ nn − ρ mm ) (cid:12)(cid:12)(cid:12) (cid:104) n | S + i | m (cid:105) (cid:12)(cid:12)(cid:12) = 1 N (cid:88) i P zi δ ( ω − ω i ) (S67)where P zi = 2Tr( S zi ρ ). This last equation exhibits a connection between EPR and polarization profiles for non-interacting spins. III. BATH-INDUCED LOCALIZATION BY TEMPERATURE VARIATION: ADDITIONALNUMERICAL RESULTSA. Thermal mixing and spin temperature behavior It has been shown in Refs. [41–43] within the secular GKLS equation and the Hilbert approximation that thereexists a remarkable situation of thermal mixing where the driven system behaves as an effective equilibrium steadystate. It can be thought as an effective Gibbs ensemble with two parameters conjugated to the two conservedquantities of the isolated (microcanonical) system with Hamiltonian H S , which are the energy and the total spinalong z . We note ρ stat the steady-state distribution, which is then diagonal in the eigenstate basis: ρ stat nn ∝ e − β S ( ε n − hs zn ) (S68)where β S is the so-called spin temperature conjugated to the energy and h is the effective magnetic field conjugatedto the spin along z .In the Hilbert approximation Eq. (S52), the transition rates define several timescales. In the following forsimplicity we assume that the two S zi transitions timescales are equal, and we note their common value T z = T (0) = T ∗ . We shall comment later what happens when they differ. T = 2 T ( ω e ) is the relaxation time of thesystem , and therefore can be considered as the longest timescale. If we assume T z (cid:28) T then the main timescalesare: • T the bath relaxation timescale induced by S xi or S yi flipping of the spins with associated energy ω ≈ ± ω e .At T = 1 . h ( ω (cid:39) ω e ) (cid:39) . h ( ω (cid:39) − ω e ) ≈ . • T z the bath relaxation timescale due to S zi flips with energy | ω nm | (cid:28) ω e . Such flips do not particularly favorany sign of the spin i , as h ( ω nm ) (cid:39) / ω nm . • The microwave time 1 W MW nm = 1 + ( T nm ∆ ω nm ) ω T nm (cid:12)(cid:12) (cid:104) n | S x | m (cid:105) (cid:12)(cid:12) (S69)defined in Eq. (S53). These timescales depend on n and m but the shortest ones ( i.e. most efficient transitions)are such that (cid:12)(cid:12) (cid:104) n | S x | m (cid:105) (cid:12)(cid:12) is not too low and ∆ ω nm (cid:39) T nm ∼ T z , this timemay be assessed through a simple Lorentzian rate1 T MW (∆ ω ) = ω T z ωT z ) (S70)of width ∆ ω ∼ /T z and minimal typical value T MWmin = ( ω T z ) − . • The Thouless time T D is the typical dephasing timescale 1 / ∆ ω nm in the coherence dynamics Eq. (S48), whichis roughly given by the interaction and disorder timescales T D = min(1 /U, / ∆ ω e ).This phase is characterized by an inhomogeneity of the polarization profile between the spins when the spintemperature is low. In a nutshell this occurs when the microwaves are effective in a short frequency range, whichthus resonate with a small number of transitions ω nm . Consequently they affect a small number of spins. But asshown in Refs. [41–43], thermal mixing happens when the interactions strong enough, typically where the system The factor 2 is conventionally defined to recover the Bloch so-lution [52] by solving exactly the single-spin case N = 1 with H S = ω e S z , Tr (cid:16) S z ρ stat (cid:17) = M T T MW , T MW = 1 + ( ω e − ω MW ) T T ω , T = 14 T (0) + 14 T ( ω e ) ,M = − 12 tanh (cid:18) βω e (cid:19) = Tr (cid:16) S z e − βω e S z (cid:17) Tr e − βω e S z . P ( ω i ) with ω i = ω e + ∆ i the typical frequency of spin i .A difficulty in observing the spin-temperature behavior is that the same timescale T z controls several competingphysical processes. It is the main timescale in the nonsecular terms for T z (cid:28) T . In the secular terms (S48), it (i) enters in the microwave timescales as is clearly visible in the Hilbert approximation (S52) (ii) it enters theevolution of the diagonal elements, inducing the homogenization of the polarization profile as these S z transitionsflip equally spins + or − (for this reason it is called spectral diffusion ) (iii) participates out of the diagonal in thedecoherence process through T nm .Within the nonsecular GKSL equation (S46), we are able to recover a marked spin-temperature behavior underthe following necessary conditions:1. ∆ ω e (cid:46) U (S71)meaning that the spin-temperature phase needs enough interactions for ETH to hold. Note that especially if N is small the inequality is not sharp.2. T z (cid:29) T D (Thouless time) (cid:39) min (cid:18) U , ω e (cid:19) (S72)In the coherence dynamics, the dephasing dominate the nonsecular terms and the Hilbert approximation isvalid.3. T T z ω (cid:38) ⇔ T MWmin (cid:46) T (S73)The microwaves are effective at least for the resonant gaps (otherwise the system would thermalize withBoltzmann distribution).4. ω (cid:114) T T z (cid:46) max( U, ∆ ω e ) (S74) ω (cid:112) T /T z is the maximal frequency range for which T MW (∆ ω ) < T under the previous condition Eq. (S73), i.e. the frequency window over which the microwaves are effective. This range needs to be narrower thanthe typical frequency range max( U, ∆ ω e ) which is probed in the EPR experiment ( i.e. frequencies around ω e which correspond to gaps of S + i ). Otherwise, most transitions are irradiated and all spins feel the microwaves,which has a tendency to thermalize them at infinite temperature ( i.e. Tr( S zi ρ stat ) = 0). This is needed tocreate some inhomogeneity within the spins.5. T z (cid:46) T (S75)This emphasizes that if T is too large then spectral diffusion will dominate and the profile will be homogeneous(high spin temperature). A way to keep spectral diffusion low is to make the two timescale differ such that T (0) (cid:28) T ∗ (S76)Indeed T (0) enters only in T nm (S19), i.e. in the decoherence process off the diagonal of Eq. (S48), while T ∗ appears in the bath transitions (S19) on the diagonal which determine directly the steady-state values.6. ω T z (cid:38) ⇔ T MWmin (cid:46) T z (S77)This is again to prevent spectral diffusion to destroy the spin-temperature behavior.The spin-temperature shapes of the EPR spectrum are displayed in Fig. 3 of the main text. In Fig. S6 we showthe corresponding polarization profiles in the ETH phase.8 ω e + ∆ i (2 π GHz) ω e + ∆ i (2 π GHz)FIG. S6. Polarization profiles (red : Hilbert, black : nonsecular) for N = 10 spins with ∆ ω e = 5 · π MHz and U = 0 . · π MHz (ETH eigenstates). The bath and microwave parameters are given in Table 1 of the main text ( T = 10 (2 π GHz) − , T ∗ = T (0) = 10 (2 π GHz) − ). The vertical dashed line points out the microwave frequency. The horizontal blue line is theBoltzmann prediction. Full dots are the numerical results ; hollow circles are obtained by the spin-temperature ansatz (S78).(left) β − = 1 . β − = 12 K. Here the bath timescales are short and the full (nonsecular) dynamics gets localized : only thenear-resonant spins feel the microwaves, as if they were non interacting, causing the hole burning. Each profile has beenaveraged over 1000 realizations of disorder and interactions for the nonsecular case and 3000 realizations for the Hilbertcase. B. Estimating the spin-temperature ansatz parameters from simulation data Here we give details about the numerical computation of the spin-temperature ansatz parameters. The spin-temperature ansatz is ρ ans nn ( β s , h ) ∝ e − β s ( ε n − hs zn ) (S78)We solve numerically for ( β s , h ) to match the average energy and polarization in the stationary state [42] Tr( H S ρ stat ) = N (cid:88) n =1 ε n ρ ans nn ( β s , h )Tr( S z ρ stat ) = N (cid:88) n =1 s zn ρ ans nn ( β s , h ) (S79)through an iterative procedure, starting from the initial guess β s = 2 N ∆ ω e argth P statn h = ω MW (S80)where the nucleus polarization is [69] P statn = ´ d ω d ( ω ) d ( ω + ω n )[ P z stat ( ω ) − P z stat ( ω + ω n )] ´ d ω d ( ω ) d ( ω + ω n )[1 − P z stat ( ω ) P z stat ( ω + ω n )] (cid:39) (cid:80) N − i =1 (cid:2) Tr( S zi ρ stat ) − Tr( S zi +1 ρ stat ) (cid:3)(cid:80) N − i =1 (cid:104) − S zi ρ stat )Tr( S zi +1 ρ stat ) (cid:105) (S81) ω n is the nucleus Zeeman gap. d ( ω ) is the disorder distribution at energy ω ; the second equality comes fromconsidering the disorder as uniform and the typical frequency between two spins ω i − ω i +1 ≈ ω n , where ω i = ω e +∆ i is the frequency of spin i in a non-interacting picture. We finally check that for all spins Tr( S zi ρ ans ) = Tr( S zi ρ stat )in order to control that the ansatz reproduces well the data.Here the polarization of a coupled nucleus enters explicitly. Note that in this work we focused on the electronicspins only, as their out-of-equilibrium steady state is the crucial feature concerning DNP in the thermal mixingregime. The electronic polarization is indeed transferred to the nuclei, the polarization of interest in DNP in fine .One can actually add a nucleus to the spin system to check the polarization transfer; this has been studied inRef. [42]. C. Hole burning A very different situation where spins do not have collective behavior is for example when the disorder isstrong [41–43], in the phase where the eigenstates are many-body localized. Then only spins resonant with the9microwave feel them and acquire a vanishing polarization, while the others are unaffected by the driving, as if theywere independent spins. The polarization profile results in a hole burning shape. We showed in this article that ifthe coupling to the bath is strong enough, the nonsecular terms may induce this hole burning shape as well in theETH phase, as if spins resonant with the microwaves were again localized with respect to the others.The necessary conditions for this bath-induced localization to take place are 1, 3 and 4 of Sec. III A but nowthe nonsecular terms need to be as important as possible to dominate over dephasing. This means we take T ∗ = T (0) = T z with T z (cid:46) T D (Thouless time) (cid:39) min (cid:18) U , ω e (cid:19) (S82)Note that the rate of the 1 /T nm decoherence term is as well controlled by T z .The hole burning at high temperature is visible in the ETH phase in Fig. S6(right) for the polarization and inFig. 3(bottom) of the main text for the EPR spectrum. In Fig. S7 we show that if the Hamiltonian eigenstatesare MBL, one gets, independently of the chosen dynamics and temperatures, a hole burning behavior. (a) (c) ω (2 π GHz) ω e + ∆ i (2 π GHz)(b) (d)FIG. S7. Numerical profiles (red : Hilbert, black : nonsecular) for N = 10 spins with ∆ ω e = 5 · π MHz and U = 0 . · π MHz(MBL eigenstates). The bath and microwave parameters for each temperature are given in the main text. The verticaldashed line points out the microwave frequency. Average is taken over 1000 realizations of disorder and interactions for thenonsecular case and 10000 realizations in the Hilbert case. β − = 1 . β − = 12 K in the bottomfigures. (a)(b) EPR spectra. Frequency bins are ten times larger in the nonsecular case than in the Hilbert case. (c)(d)Polarization profiles. The horizontal blue line is the Boltzmann prediction.Here we note a hole burning shape for both temperatures. There is no qualitative effect of nonsecularity : nonsecular termsare small for MBL eigenstates as mentioned in Sec. I F. In the data exhibited in the main text we have chosen T to be the longest timescale (the other timescale of thenonsecular terms being T ∗ ). In Fig. S8 we display the EPR and polarization profiles obtained for the same systemas in Fig. 3(bottom) of the main text (see Fig. S6) but with T (cid:28) T ∗ . We obtain again a breakdown of thermalmixing.In Fig. S9 we show that one can get an even more impressive localization effect from the nonsecular GKSLequation without disorder. This is achieved for N = 12 spins by pushing the bath relaxation times to less realisticvalues and irradiating only a single spin (namely, the third one). In the Hilbert approximation, the 12 spinsshare the same polarization owing to the strong interaction. All spins feel the microwaves. But in the nonseculardynamics, despite the strong interaction, the third spin is put to infinite temperature (Tr( S z ρ stat ) = 0) whereasall other spins are at Boltzmann equilibrium.0 ω (2 π GHz) ω e + ∆ i (2 π GHz)FIG. S8. Polarization profiles (red : Hilbert, black : nonsecular) for N = 10 spins with ∆ ω e = 5 · π MHz and U =0 . · π MHz (ETH eigenstates) β − = 12 K. The bath and microwave parameters are T = 10 (2 π GHz) − = 1 . µ s, T ∗ = 5 · (2 π GHz) − = 80 µ s, T (0) = 2 · (2 π GHz) − = 3 . µ s. (left) EPR spectrum. (right) Polarization profiles. Thevertical dashed line points out the microwave frequency. The horizontal blue line is the Boltzmann prediction. Each profilehas been averaged over 1000 realizations of disorder and interactions for the nonsecular case and 3000 realizations for theHilbert case. ω e + ∆ i (2 π GHz)FIG. S9. Polarization profile for a setup exhibiting a very strong localization, for N = 12 spins at β − = 12 K (red : Hilbert,black : nonsecular). The data corresponds to a single typical sample where disorder is negligible (∆ ω e = 10 − · π MHz,essentially to ensure that the gaps are non degenerate) and the interaction is U = 2 π MHz. The microwaves irradiate only the third spin, with frequency ω MW = ω e and strong amplitude ω = 0 . · π GHz. The nonsecular terms are very effectivedue to the very low value of T (0) = T ∗ = 0 . · (2 π GHz) − = 10 − / π s ; T = 10 − / π s. D. Spectral properties of the stationary state In this section we provide additional data concerning the spectral properties of the steady-state density matrix ρ stat in the Hilbert or nonsecular cases. We emphasize that, for all the N = 10 data, the sample-to-samplefluctuations are quite well concentrated on the averages displayed.In Fig. S10 we provide the eigenvalue distribution of the stationary density matrix, for parameters where thenonsecular dynamics leads to a hole burning while the Hilbert one yields a spin-temperature shape. This is toshow that the distribution is not dominated by a small number of eigenstates ; all of them are relevant. Moreover,the distribution in both dynamics are similar : the qualitative difference in steady-state observables does not comefrom this population distribution but finds an explanation in the eigenstates’ statistics, as explained in the maintext.It is pointed out in the main text that the when when the bath induces localization in the system, the eigenstatesreached by the density matrix in the long-time limit are different from the Hamiltonian ones. They do not satisfyETH : are less entangled and more akin to MBL eigenstates. Now we focus on the more extreme case introducedat the end of Sec. III C: 12 homogeneous spins, only spin 3 irradiated with strong nonsecular effects. Here thelocalization is “complete” : spin 3 decouples from the rest of the system, in spite of the strong interactions. InFig. S11 we demonstrate that the eigenstates of the density matrix is not anymore the ETH Hamiltonian eigenstatesbut are pure states factorized on each individual spin basis, i.e. all S zi are good quantum numbers to describethis eigenbasis ; it corresponds to the Hamiltonian basis for non-interacting spins U = 0. Here the measurementsperformed by the bath modes wipe out all interaction in the many-body system. The entanglement entropy ofeach eigenstate is thus zero. The eigenvalue distribution of the matrix is much closer to the Boltzmann one thanthe Hilbert one.1 λ FIG. S10. Eigenvalues (populations) (cid:104) λ | ρ stat | λ (cid:105) of the stationary density matrix ( | λ (cid:105) are eigenvectors), classified on thehorizontal axis by increasing energy (cid:104) λ | H S | λ (cid:105) . Parameters correspond to those of Fig. S6 for N = 10 spins. Vertical gray linesdelimit sectors of constant polarization (for Hamiltonian eigenstates). Blue : Boltzmann equilibrium ρ stat ∝ exp( − βH S ).Red : Hilbert solution. For the latter cases, the eigenvectors | λ (cid:105) = | n (cid:105) are those of the Hamiltonian. Black : nonsecularsolution. (a) (c) λ λ (b) (d)FIG. S11. The setting corresponds to the fully localized one with 12 spins introduced at the end of Sec. III C. Index λ isclassified on the horizontal axis by increasing energy (cid:104) λ | H S | λ (cid:105) ( | λ (cid:105) are eigenvectors of ρ stat ). Blue : Boltzmann equilibrium ρ stat ∝ exp( − βH S ). Red : Hilbert solution. For the latter cases, the eigenvectors | λ (cid:105) = | n (cid:105) are those of the Hamiltonian.Black : nonsecular solution. Vertical gray lines delimit sectors of constant polarization (for Hamiltonian eigenstates). (a)Eigenvalues (populations) (cid:104) λ | ρ stat | λ (cid:105) of ρ stat . (b) Entanglement entropy of each eigenvector. Partial trace is taken overspins 6 to 12. (c) Expectation value of the observable S z . (d) Same for S z , qualitatively close to any other spin. IV. A TOY MODEL FOR ZENO LOCALISATION In this section, we analyze a toy model introduced in Fig. 2 of the main text. It presents similar features to themany-body DNP model with the advantage of possessing an analytical solution. We put hats on some operatorsto avoid possible confusion with scalar quantities.Let us consider a three-level system (states labeled 0, 1, 2) subject to a monochromatic microwave perturbationthat only couples the ground state | (cid:105) and the state | (cid:105) . The system is coupled to two uncorrelated harmonic bathsthat allow exchange of energy between the ground state and the excited states. The Hamiltonian reads H = H S + H MW ( t ) + H int + H B == V ∗ e iω MW t + ˆ B † ˆ B † V e − iω MW t + ˆ B ω e − ∆ J ˆ B J ω e + ∆ + (cid:88) α =1 , (cid:88) k Ω k ˆ a † αk ˆ a αk (S83)ˆ B α = (cid:88) k g αk ˆ a αk (S84)The Born-Markov master equation (S7) (no secular approximation is involved) for the system density matrix ρ inthe Schr¨odinger representation has the form dρdt = − i [ H S + H MW , ρ ] − ∞ ˆ d τ Tr B (cid:40)(cid:20) H int , (cid:104) e − i ( H S + H B ) τ H int e i ( H S + H B ) τ , ρ ⊗ ρ B (cid:105)(cid:21)(cid:41) . (S85)It does not depend on the choice of basis in the system subspace and thus does not privilege the eigenbasis of H S .Assuming J, ∆ (cid:28) ω e , we approximate e − i ( H S + H B ) τ H int e i ( H S + H B ) τ ≈ (cid:88) α =1 , | α (cid:105)(cid:104) | (cid:88) k g αk ˆ a αk e i (Ω k − ω e ) τ + H . c . ≡ (cid:88) α =1 , | α (cid:105)(cid:104) | ˆ B α ( − τ ) e − iω e τ + H . c . (S86)with ˆ B α ( τ ) = (cid:80) k g αk e − i Ω k τ ˆ a αk the dynamical evolution of the bath modes. When performing the trace over thebath modes, assuming the bath has no fine structure on the scale J, ∆, we must only evaluate the following bathcorrelators at frequency ω e (and their complex conjugates): ∞ ˆ d τ e iω e τ (cid:68) ˆ B † α (0) ˆ B β ( τ ) (cid:69) B = δ αβ (cid:88) k i | g αk | n αk ω e − Ω k + i + ≡ δ αβ [ γ α ¯ n α + iδ α ] , (S87) ∞ ˆ d τ e iω e τ (cid:68) ˆ B α ( τ ) ˆ B † β (0) (cid:69) B = δ αβ (cid:88) k i | g αk | ( n αk + 1) ω e − Ω k + i + ≡ δ αβ [ γ α (¯ n α + 1) + iδ (cid:48) α ] , (S88) n αk = ( e β Ω k − − and ¯ n α = ( e βω e − − are Bose-Einstein distributions, γ α ∝ | g αk | k ↔ ω e and the Lamb-shifts δ α , δ (cid:48) α can be expressed through Cauchy principal values. The resulting master equation has the GKSL form˙ ρ = − i [ H S + H LS + H MW , ρ ] + D ( ρ ) + D ( ρ ) H LS ≡ (cid:88) α =1 , δ (cid:48) α | α (cid:105)(cid:104) α | − δ α | (cid:105)(cid:104) | D α ( ρ ) ≡ L α ρL † α + (cid:101) L α ρ (cid:101) L † α − (cid:110) L † α L α + (cid:101) L † α (cid:101) L α , ρ (cid:111) L α = (cid:112) γ α ¯ n α | α (cid:105)(cid:104) | (cid:101) L α = (cid:112) γ α (¯ n α + 1) | (cid:105)(cid:104) α | (S89)Let us make here a few comments on this GKSL equation to make contact with our DNP model:1. This equation is identical to the one obtained through Ref. [58]’s approach employed in the DNP model.Indeed one can project the usual Born-Markov equation (S7) on eigenstates and get the same jump operatorsfor ∆ , J (cid:28) ω e . Alternatively, one can directly compute Ref. [58]’s nonsecular jump operators (defined inEq.(6) of the main text) by redefining slightly the system in order to satisfy the hypothesis that H int is writtenas a sum of Hermitian operators : H int = (cid:80) α =1 , ˆ S α ⊗ ˆ B α with ˆ S α = | α (cid:105)(cid:104) | + H . c . , ˆ B α = (cid:80) k g αk ˆ b αk + H . c . ,yielding the same result.2. The toy model corresponds to a variant of our DNP model for N = 2 spins truncated to the three lowest levels, i.e. taking | ++ (cid:105) = 0 or density matrix elements (cid:104) λ | ρ | ++ (cid:105) for any vector | λ (cid:105) . More precisely, the relationshipis U ≡ J with fixed local disorder ∆ i ≡ ± 2∆ and energy levels shifted so that H S |−−(cid:105) = 0, only the spin 1is irradiated at frequency ω MW (with V ≡ ω / √ 2) and the bath do not couple to S zi ( i.e. T ∗ , T (0) → ∞ ); forexample one can consider H int = (cid:80) i S xi ⊗ B xi which defines the parameters γ , ≡ γ ( ω e ) / 8. In this setting,the population of levels 1 and 2 in the 3-level model plotted in Fig. 2 of the main text gives directly thesteady-state polarization of each spin: Tr( S zi ρ ) ≡ ρ ii − .3In the following we neglect the Lamb-shifts δ . Explicitly in components,˙ ρ = iV ∗ e iω MW t ρ ∗ − iV e − iω MW t ρ + iJ ( ρ − ρ ∗ ) + 2 γ ¯ n ρ − γ (¯ n + 1) ρ , (S90)˙ ρ = iJ ( ρ ∗ − ρ ) + 2 γ ¯ n ρ − γ (¯ n + 1) ρ , (S91)˙ ρ = iω e ρ + iV ∗ e iω MW t ( ρ − ρ ) − i ∆ ρ + iJρ − [ γ (2¯ n + 1) + γ ¯ n ] ρ , (S92)˙ ρ = iω e ρ − iV ∗ e iω MW t ρ + iJρ + i ∆ ρ − [ γ (2¯ n + 1) + γ ¯ n ] ρ , (S93)˙ ρ = 2 i ∆ ρ − iV e − iω MW t ρ + iJ ( ρ − ρ ) − [ γ (¯ n + 1) + γ (¯ n + 1)] ρ . (S94)Looking for the solution with time-independent ρ , ρ , ρ and ρ , ρ ∝ e iω MW t (we rename ρ α → ρ α e − iω MW t ),we arrive at a time-independent linear system. Let us for simplicity take V = V ∗ , γ = γ ≡ γ , and set ¯ n α = 0( i.e. zero temperature). Then the system has the form0 = V ( ρ − ρ ∗ ) − J ( ρ − ρ ∗ ) − iγρ , (S95)0 = J ( ρ − ρ ∗ ) − iγρ , (S96)0 = ( ω MW − ω e + ∆ − iγ ) ρ − Jρ + V (2 ρ + ρ − , (S97)0 = − Jρ + ( ω MW − ω e − ∆ − iγ ) ρ + V ρ , (S98)0 = V ρ − iγ ) ρ + J ( ρ − ρ ) . (S99)Let us solve the last three equations (we define ∆ ω = ω e − ω MW for compactness): ρ D = − J V ( ρ − ρ ) + [ V + 2 V (∆ + iγ )(∆ ω + ∆ + iγ )](1 − ρ − ρ ) ,ρ D = − JV (∆ − ∆ ω − iγ )( ρ − ρ ) − JV (∆ + iγ )(1 − ρ − ρ ) ,ρ D = [ J (∆ ω + ∆ + iγ )(∆ ω − ∆ + iγ ) − J ]( ρ − ρ ) − JV (1 − ρ − ρ ) ,D = 2 J (∆ + iγ ) + V (∆ ω − ∆ + iγ ) − iγ )(∆ ω + ∆ + iγ )(∆ ω − ∆ + iγ ) . (S100)In the linear response regime, we seek ρ , ρ = O ( V ), ρ , ρ , ρ = O ( V ) and neglect V in the denominator,which gives ρ = V [(∆ ω + ∆) + γ ](∆ ω − ∆ − J − γ ) + 4 γ ∆ ω , ρ = V J (∆ ω − ∆ − J − γ ) + 4 γ ∆ ω . (S101)This is the expression plotted in the Fig. 2 in the main text. We have good mixing at γ, ∆ (cid:28) J , Andersonlocalisation for γ (cid:28) J (cid:28) ∆, Zeno localisation at ∆ < J (cid:28) γ , and weak incoherent coupling for J (cid:28) γ, ∆. V. MICROSCOPIC MODEL OF THE ELECTRONIC SPINS AND BATH In this section we comment on the choice of the microscopic Hamiltonian of the whole system. However thepresent treatment of the dissipation modes is crude and does not allow to recover realistic orders of magnitudesfor relaxation timescales ; nevertheless it shows a simple microscopic example with bath rates that are power-lawincreasing with temperature, as expected from NMR experiments. A. Coupling between spins and the magnetic field In absence of lattice vibrations, the electron spins are frozen in an amorphous matrix. We denote the positionof a given spin by R . The spin-orbit interaction between this spin and the external magnetic field depends onthe orientation of the radical with respect to the field, and in general is thus written as a tensorial coupling − g zµ µ B BS µi [49, Chap. VI]-[70] where g zµ = g e δ zµ + g (0) zµ (S102)The first term in Eq. (S102) is the dominant isotropic contribution, resulting in the Zeeman part of the spinHamiltonian H S , i.e. (cid:80) i ω e S zi = ω e S z with ω e = − g e µ B B the Zeeman frequency, related to the electron Land´e g -factor g e ( (cid:39) − 2) through µ B = e m e . e is the unit charge, m e the electron mass and µ B the Bohr magneton.In practice we take ω e = 93 . · π GHz, meaning B (cid:39) . 35 Tesla, a standard value for DNP. The second term inEq. (S102) is the so-called g -factor anisotropy of the disordered sample, which depends on the radical orientationand therefore appears as a random quantity. This term contributes to H S as (cid:80) i ∆ i S zi with ∆ = − g (0) zz µ B B . Notethat we discarded directions µ (cid:54) = z : when the magnetic field is large one can resort to the secular approximationof the Hamiltonian, which consists in keeping only the terms conserving the total polarization along z . Thereason in that hybridization between different polarization sectors is very weak in perturbation theory [49, Chap.IV.II.A]-[42, 51]. Repeated Greek indices are summed over in the whole section.Note that the magnetic field is along z , hence the z index. For notational simplicity we omit explicit reference to the spin i , implicitly born by g (0) zz . B. Vibrational modes of the embedding material The position of an electron spin in the amorphous matrix is r = R + u ( R ) with R an equilibrium position,and u ( R ) describes a small vibrational motion around it. The latter motion affects all space-dependent quantities,such as the Zeeman interaction which is the strongest term in the Hamiltonian. The field u is expected to varyslowly on the scale of the electron distances for extended vibrational modes. The distance between two particles 1and 2 in the matrix is indeed r = r − r = R − R + ( ∂ µ u )( R − R ) µ at first order in the u derivatives. Thetensorial coupling in Eq. (S102) is now modified by the vibrations as [71, Chap. 22]-[72]. g zµ = g e δ zµ + g (0) zµ + g (1) zµγδ ∂ γ u δ + g (2) zµγδγ (cid:48) δ (cid:48) ∂ γ u δ ∂ γ (cid:48) u δ (cid:48) + . . . (S103)Vibrational modes in the glass originate from several processes [73, 74], but in the following we shall restrictourselves for simplicity to model them as low-energy excitations arising from acoustic phonon modes . In theglassy sample, the arrangement of the different atoms is not periodic. This implies a continuous set of wavevectors;we nevertheless use the standard theory of phonons of a periodic lattice for convenience, as it should not affectmuch the results. One can write the quantized displacement field [71] u ( R ) = 1 √ N (cid:88) k ,s e k /k,s (cid:112) m Ω k ,s a k ,s e i k · R + H.c. (S104) s = 1 , , k are the wavevector (quantized in the first Brillouin zone due to the assumedperiodicity of the lattice), m is the mass of the glassy molecule Ω k ,s are the phonon frequencies, e k /k,s arepolarization unit eigenvectors, and a k ,s , a † k ,s are phonon annihilation and creation operators. The bath Hamiltonianis thus a collection of harmonic oscillators H B = (cid:80) k ,s Ω k ,s a † k ,s a k ,s .In the following we consider separately the one- and two-phonon processes as they are incoherent owing to Wick’stheorem [67], i.e. their contribution to the correlation function (cid:10) B µi ( t ) B µi (cid:11) B is additive. 1. Direct process The direct process concerns the exchange of a single phonon between the bath and the system. It is due to thefirst-order interaction between the spin and the bath modes in Eq. (S103) involving the tensor g (1) , substituting ∂ γ u δ via Eq. (S104). The interaction Hamiltonian is thus of the form (S2) where B µi a linear combination of annihilationand creation operators. The equilibrium bath correlation function (S15) then involves only quadratic correlators inthe a k ,s . For simplicity we drop tensor indices and the polarization vectors, as these factors only contribute O (1)proportionality constants. The calculation is standard (see (S87)-(S88)) using the dispersion relation Ω k = kv with v the sound velocity. We get in the continuous limit N → ∞ T ( | ω | ) = 12 π (cid:32) g (1) g e (cid:33) ω ρ v ω coth (cid:18) βω (cid:19) (S105)where ρ = N m/L is the mass density. 2. Two-phonon processes The next term in the system-bath interaction is a two-phonon process caused by the second-order interactionbetween the spin and the bath modes in Eq. (S103) involving the tensor g (2) . It takes the form (S2) where B µi is quadratic in the annihilation and creation operators. We thus have to deal with 4-point canonical averagesapplying Wick’s theorem [67]. In the continuous limit for wavevectors and at low temperature , dropping againindices and polarization vectors, we obtain1 T ( | ω | ) = 16 π (cid:32) g (2) g e (cid:33) ω ( ρ v ) T cosh (cid:18) βω (cid:19) ˆ β | ω | / d y y (cid:16) β | ω | − y (cid:17) sinh y sinh (cid:16) β | ω | − y (cid:17)(cid:124) (cid:123)(cid:122) (cid:125) absorption of two phonons + 2 ˆ ∞ d y y (cid:16) β | ω | + y (cid:17) sinh y sinh (cid:16) β | ω | + y (cid:17)(cid:124) (cid:123)(cid:122) (cid:125) Raman process: absorption and emission (S106) One may expect this physical picture to dominate the system-bath interaction at low temperature for energies below 100GHz, based on Fermi’s golden rule which involves vibrational density of states estimates from Refs. [73, 74]. i.e. approximating pyruvic acid as a monoatomic substance. The Debye frequency Ω D must be large, i.e. β Ω D (cid:29) 3. Discussion of the one- and two-phonon bath timescales For our range of frequencies and temperature, the direct rate (S105) is ∝ ω T . The two-phonon contribu-tion (S106) is dominated by the Raman contribution roughly independent of the frequency and ∝ T . In experi-ments the dependence in frequency is not known while the dependence in temperature of T = 2 T ( ω e ) is roughly T [75, 76]. As we do not know realistic values of the ratios g ( i ) /g e , we cannot assess the order of magnitude ofthe predicted timescales, as well as the relative weight of the Raman rate with respect to the direct rate whichdetermines the temperature dependence of the bath correlation function γ ( ω ). Yet it is known experimentally thatthe relaxation time T of the system is ∼ β − ∼ g ( i ) /g e . Finally let us look at the decorrelation assumption of the bath degrees of freedom.The fact that (cid:68) B µi ( t ) B νj (cid:69) B ∝ δ µν owes to the isotropy of the material. If i (cid:54) = j the only difference is that allintegrals over wavevectors (say k ) get an additional phase factor e i k · ( R i − R j ) . Therefore decorrelation happens ifthe wavelength of the phonon is much smaller than the inter-electron distance | R i − R j | . In other words, thecriterion for decorrelation is that the phonon frequencies ω involved are such that ω (cid:29) v/ | R i − R j | i.e. ω (cid:29) vρ / where ρ e is the electron density in the material. In practice ( ρ e ∼ 10 mmol/L, v ∼ m / s) this threshold isof the same order of magnitude as ω e . We conclude that this model of ballistic phonons does not give a realisticdescription but represents a simple example with a power-law decrease of the bath relaxation timescales neededfor Zeno localization.6 [1] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol,Advances in Physics , 239 (2016).[2] E. M. Lifshitz and L. P. Pitaevskii, Statistical Physics,Part 2 , vol. 9 of L. D. Landau and E. M. Lifshitz Courseof Theoretical Physics (Pergamon Press, 1980), sec. 41.[3] J. M. Deutsch, Physical Review A , 2046 (1991).[4] M. Srednicki, Physical Review E , 888 (1994).[5] R. Nandkishore and D. A. Huse, Annual Review of Con-densed Matter Physics , 15 (2015).[6] F. Alet and N. Laflorencie, Comptes Rendus Physique , 498 (2018), ISSN 1631-0705, quantum simulation/ Simulation quantique.[7] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn,Reviews of Modern Physics (2019).[8] B. L. Altshuler, Y. Gefen, A. Kamenev, and L. S. Lev-itov, Phys. Rev. Lett. , 2803 (1997).[9] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Annalsof Physics , 1126 (2006), ISSN 0003-4916.[10] M. Pollak and B. Shklovskii, eds., Hopping Transport inSolids , vol. 28 of Modern Problems in Condensed MatterSciences (Elsevier, 1991).[11] A. Beskow and J. Nilsson, Arkiv f¨ur Fysik , 561(1967).[12] L. Khalfin, Jetp Lett , 65 (1968).[13] B. Misra and E. C. G. Sudarshan, Journal of Mathe-matical Physics , 756 (1977).[14] W. M. Itano, D. J. Heinzen, J. J. Bollinger, and D. J.Wineland, Phys. Rev. A , 2295 (1990).[15] P. Facchi and S. Pascazio, in Progress in Optics (Else-vier, 2001), pp. 147–217.[16] H.-P. Breuer and F. Petruccione, The theory of openquantum systems (Oxford University Press, 2002).[17] D. A. Huse, R. Nandkishore, F. Pietracaprina, V. Ros,and A. Scardicchio, Physical Review B , 014203(2015).[18] V. S. Shchesnovich and V. V. Konotop, Phys. Rev. A , 053611 (2010).[19] P. Barmettler and C. Kollath, Phys. Rev. A , 041606(2011).[20] D. A. Zezyulin, V. V. Konotop, G. Barontini, andH. Ott, Phys. Rev. Lett. , 020405 (2012).[21] H. Fr¨oml, A. Chiocchetta, C. Kollath, and S. Diehl,Phys. Rev. Lett. , 040402 (2019).[22] P. E. Dolgirev, J. Marino, D. Sels, and E. Demler, Phys.Rev. B , 100301 (2020).[23] A. J. Bray and M. A. Moore, Phys. Rev. Lett. , 1545(1982).[24] S. Chakravarty, Phys. Rev. Lett. , 681 (1982).[25] A. Schmid, Physical Review Letters , 1506 (1983).[26] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A.Fisher, A. Garg, and W. Zwerger, Reviews of ModernPhysics , 1 (1987).[27] L. F. Cugliandolo, D. R. Grempel, G. Lozano, H. Lozza,and C. A. da Silva Santos, Physical Review B ,014444 (2002).[28] Y. Li, X. Chen, and M. P. A. Fisher, Phys. Rev. B ,205136 (2018).[29] A. Chan, R. M. Nandkishore, M. Pretko, and G. Smith,Phys. Rev. B , 224307 (2019).[30] B. Skinner, J. Ruhman, and A. Nahum, Phys. Rev. X , 031009 (2019).[31] M. Szyniszewski, A. Romito, and H. Schomerus, Phys.Rev. B , 064204 (2019).[32] X. Turkeshi, R. Fazio, and M. Dalmonte, Phys. Rev. B , 014315 (2020).[33] C.-M. Jian, Y.-Z. You, R. Vasseur, and A. W. W. Lud-wig, Phys. Rev. B , 104302 (2020).[34] A. Zabalo, M. J. Gullans, J. H. Wilson, S. Gopalakrish-nan, D. A. Huse, and J. H. Pixley, Phys. Rev. B ,060301 (2020). [35] Y. Bao, S. Choi, and E. Altman, Phys. Rev. B ,104301 (2020).[36] L. Zhang, J. A. Reyes, S. Kourtis, C. Chamon, E. R.Mucciolo, and A. E. Ruckenstein, Phys. Rev. B ,235104 (2020).[37] X. Cao, A. Tilloy, and A. D. Luca, SciPost Phys. , 24(2019).[38] Q. Tang and W. Zhu, Phys. Rev. Research , 013022(2020).[39] S. Goto and I. Danshita, Phys. Rev. A , 033316(2020).[40] Y. Fuji and Y. Ashida, Phys. Rev. B , 054302(2020).[41] A. De Luca and A. Rosso, Physical review letters ,080401 (2015).[42] A. D. Luca, I. Rodr´ıguez-Arias, M. M¨uller, andA. Rosso, Physical Review B , 014203 (2016).[43] I. Rodr´ıguez-Arias, M. M¨uller, A. Rosso, andA. De Luca, Phys. Rev. B , 224202 (2018).[44] Z. Lenarˇciˇc, E. Altman, and A. Rosch, Physical ReviewLetters (2018).[45] Z. Lenarˇciˇc, O. Alberton, A. Rosch, and E. Altman,Phys. Rev. Lett. , 116601 (2020).[46] F. Lange, Z. Lenarˇciˇc, and A. Rosch, Nature Commu-nications (2017).[47] D. Guarin, S. Marhabaie, A. Rosso, D. Abergel, G. Bo-denhausen, K. L. Ivanov, and D. Kurzbach, The Journalof Physical Chemistry Letters , 5531 (2017).[48] V. Atsarkin, Soviet Phys.-JETP , 1012 (1970).[49] A. Abragam, The principles of nuclear magnetism , 32(Oxford university press, 1961).[50] A. Abragam and M. Goldman, Reports on Progress inPhysics , 395 (1978).[51] S. A. Smith, W. E. Palke, and J. T. Gerig, Concepts inMagnetic Resonance , 107 (1992).[52] F. Bloch, Physical Review , 460 (1946).[53] P. Schosseler, T. Wacker, and A. Schweiger, ChemicalPhysics Letters , 319 (1994).[54] J. Granwehr and W. K¨ockenberger, Applied MagneticResonance , 355 (2008).[55] Y. Hovav, I. Kaminker, D. Shimon, A. Feintuch,D. Goldfarb, and S. Vega, Physical Chemistry ChemicalPhysics , 226 (2015).[56] R. Alicki, Physical Review A , 4077 (1989).[57] R. Alicki, D. Gelbwaser-Klimovsky, and G. Kurizki,arXiv:1205.4552 (2012).[58] F. Nathan and M. S. Rudner, Phys. Rev. B , 115109(2020).[59] G. Kirˇsanskas, M. Francki´e, and A. Wacker, PhysicalReview B (2018).[60] E. Kleinherbers, N. Szpak, J. K¨onig, and R. Sch¨utzhold,Phys. Rev. B , 125131 (2020).[61] Y. Hovav, A. Feintuch, and S. Vega, Phys. Chem.Chem. Phys. , 188 (2013).[62] J. Bezanson, A. Edelman, S. Karpinski, andV. B. Shah, SIAM Review , 65 (2017),https://doi.org/10.1137/141000671.[63] Y. Saad, Iterative Methods for Sparse Linear Sys-tems (Society for Industrial and Applied Mathematics,2003).[64] N. Pottier, Nonequilibrium statistical physics: linear ir-reversible processes (Oxford University Press, 2010).[65] M. Le Bellac, Physique quantique (EDP sciences, 2012).[66] I. Rodr´ıguez-Arias, A. Rosso, and A. D. Luca, MagneticResonance in Chemistry , 689 (2018).[67] P. Coleman, Introduction to Many-Body Physics (Cam-bridge University Press, 2015).[68] F. Caracciolo, M. Filibian, P. Carretta, A. Rosso, andA. D. Luca, Physical Chemistry Chemical Physics ,25655 (2016). [69] S. C. Serra, M. Filibian, P. Carretta, A. Rosso, andF. Tedoldi, Phys. Chem. Chem. Phys. , 753 (2014).[70] L. F. Chibotaru, A. Ceulemans, and H. Bolvin, Phys.Rev. Lett. , 033003 (2008).[71] N. W. Ashcroft and N. D. Mermin, Solid state physics (Holt, Rinehart and Winston, 1976).[72] L. D. Landau and E. M. Lifshitz, Theory of Elasticity (Elsevier, 1986).[73] H. Mizuno, H. Shiba, and A. Ikeda, Proceedings of theNational Academy of Sciences , E9767 (2017). [74] L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel,and E. Flenner, Nature Communications (2019).[75] M. Filibian, S. C. Serra, M. Moscardini, A. Rosso,F. Tedoldi, and P. Carretta, Phys. Chem. Chem. Phys. , 27025 (2014).[76] M. Filibian, E. Elisei, S. C. Serra, A. Rosso, F. Tedoldi,A. Ces`aro, and P. Carretta, Physical Chemistry Chem-ical Physics18