Dynamical properties of a driven dissipative dimerized S=1/2 chain
M. Yarmohammadi, C. Meyer, B. Fauseweh, B. Normand, G. S. Uhrig
DDynamical properties of a driven dissipative dimerized S = 1 / chain M. Yarmohammadi, C. Meyer, B. Fauseweh, B. Normand,
4, 1, 5 and G. S. Uhrig Lehrstuhl f¨ur Theoretische Physik I, Technische Universit¨at Dortmund,Otto-Hahn-Strasse 4, 44221 Dortmund, Germany Institut f¨ur Theoretische Physik, Georg-August-Universit¨at G¨ottingen,Friederich-Hund-Platz 1, 37077 G¨ottingen, Germany Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA Paul Scherrer Institute, CH-5232 Villigen PSI, Switzerland Institute of Physics, Ecole Polytechnique F´ed´erale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
We consider the dynamical properties of a gapped quantum spin system coupled to the electricfield of a laser, which drives the resonant excitation of specific phonon modes that modulate themagnetic interactions. We deduce the quantum master equations governing the time-evolution ofboth the lattice and spin sectors, by developing a Lindblad formalism with bath operators providingan explicit description of their respective phonon-mediated damping terms. We investigate the non-equilibrium steady states (NESS) of the spin system established by a continuous driving, delineatingparameter regimes in driving frequency, damping, and spin-phonon coupling for the establishmentof physically meaningful NESS and their related non-trivial properties. Focusing on the regime ofgeneric weak spin-phonon coupling, we characterize the NESS by their frequency and wave-vectorcontent, explore their transient and relaxation behavior, and discuss the energy flow, the systemtemperature, and the critical role of the type of bath adopted. Our study lays a foundation for thequantitative modelling of experiments currently being designed to control coherent many-body spinstates in quantum magnetic materials.
I. INTRODUCTION
Both the advent of powerful new laser sources and theincreasing demand for next-generation magnetic devices,required to power the information revolution, are focus-ing intensive research efforts on time-dependent phenom-ena in condensed matter. On the laser side, x-ray free-electron laser sources in the US and Europe now allowthe “ultrafast” probing of materials on the femtosecondtimescales of their fundamental electronic and magneticprocesses. On the device side, the immediate target isdesigner materials for antiferromagnetic (AF) spintronics[1, 2], to enable the writing, storage, and reading of large-scale classical magnetic information with factor-1000 im-provements over the current levels of speed and powerconsumption. Already on the horizon, however, is thedevelopment of magnetic materials as a route to encod-ing and manipulating quantum information, and indeedquantum information processing in systems with stronginteraction energies would ensure very high-frequency op-eration at the lowest possible dissipation.The concept of laser driving generalizes the pump-probe paradigm from simple pulse-delay schemes to theimprinting of arbitrary dynamics (within the limits offield control). The laser excitation of quantum systemshas generated theoretical proposals for uniquely out-of-equilibrium states of matter, including non-equilibriumsteady states (NESS) [3, 4], non-equilibrium topologicalstates [5], and many-body localization (MBL) [6, 7]. Todate these ideas have been tested largely on systems of ul-tracold atoms [8–11], where the laser controls the “opticallattice” on which the atoms reside [12]. The undeniablebeauty of both the physical concepts and the technolog-ical achievements aside, these systems are neither very large nor very readily miniaturized.Laser facilities operating on the energy and ultrafasttime scales of condensed-matter systems have been de-ployed recently to observe a wide array of novel phe-nomena in graphene [13], superconductors [14], charge-density-wave materials [15, 16], and correlated insulatorsnear their metallic transition [17, 18]. Beyond inducing,enhancing, or destroying a symmetry-broken state, a keyfocus of these experiments has been the high-frequencyFloquet regime, where steady laser driving can inducenew topological states [19, 20], the “time crystal” [21],or more generally allow the “Floquet engineering” of theelectronic bands [22, 23].While any material can be laser-driven, the key ques-tion is whether this driving creates a coherent quantumstate [4]. Some of these new phenomena, notably photo-enhanced superconductivity [14], occur because the laserdrives particular phonon excitations of the lattice hostingthe electrons. Because strong laser driving can lead tovery high populations of any targeted mode, exploitingthe anharmonic part of the lattice restoring force leads tothe concept of nonlinear phononics [24–26]. However, thephonon ensemble determines the temperature, and henceheating of the system is a fundamental issue in determin-ing whether any of these novel laser-driven phenomena,and particularly their quantum nature, can survive be-yond the initial ultrafast laser pulses.Among the extensive body of theoretical studies ofnon-equilibrium quantum systems are analyses of short-time transient behavior caused by quenches [27–29], in-cluding those due to laser pulses [30, 31], and of long-time thermalization behavior [28, 32]. Ideas from (near-)integrable systems include MBL, which is known atleast in one dimension (1D) [7], and pre-thermalization a r X i v : . [ c ond - m a t . s t r- e l ] S e p [33], while numerous studies have explored the Floquetregime [34]. Of the many numerical approaches toquenched or driven models, one of the most successfulis non-equilibrium dynamical mean-field theory (DMFT)[35, 36], which has been applied to many problems incold atoms [37] and condensed matter [38–40]. However,these studies are largely restricted to fermionic systemsand focus mostly on leading qualitative effects due to in-trinsic system dynamics, rather than on the dynamics inthe presence of dissipation. By contrast, a realistic NESSrequires a path for outflow of the injected energy [4]. Oneapproach is to use a phenomenological Gilbert damp-ing [41], while a more sophisticated treatment invokesa damping bath that maintains a constant temperature[42, 43]. In general, a dissipative (open) system is treatedby the Lindblad formalism [44], in which damping is pro-vided by bath operators whose Hamiltonian dynamics arenot required to formulate the equations of motion gov-erning the time-evolution of physical observables in theHeisenberg representation [45, 46]. Providing a quanti-tative treatment of a driven quantum system subject toLindblad dissipation processes is the aim of the currentwork.For this purpose we will focus on quantum magneticsystems, which historically have provided a clean, read-ily realized, low-dissipation test bed for many conceptsin condensed-matter and statistical physics. The smallnumber and unique behavior of the spin degrees of free-dom lead to exact solutions including the Heisenberg spinchain, the transverse-field Ising model, and the Kitaevmodel. In non-equilibrium physics, idealized (and oftenintegrable) spin-chain models as the Hamiltonian part ofa Lindblad system have provided the framework for il-lustrating NESS [47], MBL [48, 49], Floquet prethermal-ization [50], and dynamical quantum phase transitions[51], as well as lending themselves very well to numericalinvestigation. With a view to future device application,single spins have long been considered as excellent candi-date qubits and the application of suitable laser controlschemes [4] has been attempted in ensembles of quan-tum dots [52, 53]. The entangled quantum many-bodystates available in magnetic materials present not only arich variety of options for encoding (protected) quantuminformation, keywords including (topological) magnonics[2, 54–57], quantum spin liquids (QSLs) [58], and mag-netic textures such as vortices [59] and skyrmions [60],but also many routes for exploiting intrinsic interactions[61–63] or extrinsic materials-design flexibility [64] to ob-tain “handles” for manipulating magnetism using laserlight [65–67].The reason why insulating quantum magnets are arelative late-comer to the game of laser excitation andpump-probe physics is the weak direct coupling of lightto spin. In general one may consider four routes for thecreation of magnetic excitations by incident light. (1)The response of metallic magnetic systems is usually de-scribed in terms of the inverse Faraday effect, and thismechanism remains (in the form of virtual electronic pro- cesses) in insulators. It was exploited recently [68] tostudy the coherent transport of GHz precession modes ofthe magnetization over 100 µ m distances in ferromag-netic iron garnet films. (2) At the intrinsic frequen-cies of magnetic modes in condensed matter, which areof order 1 THz, processes by which one photon createsone magnon depend on anisotropies in the spin Hamil-tonian and thus are rather weak in the most familiarquantum magnetic materials, whose magnetic ions are3 d transition metals. However, they are present in type-II multiferroics and other systems with finite magneto-electric [69] and thermomagnetic coupling [70], and onesuch anisotropy was exploited in a recent discussion of alaser-pumped spin chain as a test case for a GeneralizedGibbs Ensemble approach to near-integrable dissipativesystems [42].(3) The mechanism invoked most commonly in con-densed matter is the Peierls coupling of the magneticfield to the electron phase. In insulating magnets, thisleads to Raman-type processes where one photon excitestwo magnon modes [71, 72] and thus spin is conserved.For this type of process, the light frequency, ω , shouldexceed the band gap, meaning that it should be a sig-nificant fraction of the on-site Coulomb repulsion, U , ofthe electrons being excited virtually; because U is of or-der 5 eV, the incident light should be around the visiblerange. At lowest order, incident photons with frequency ω modify U to U − ω or U + ω , thereby affecting themagnetic (super)exchange interaction. All of routes (1-3) are nonresonant, and none make direct reference tothe fact that spins are hosted on a lattice, although thelattice may create the anisotropic magnetic environmentrequired for route (2). However (4), because the latticegeometry is fundamental to the magnetic interactions,selective excitation of specific phonon modes would pro-vide direct control of magnetic interactions through amechanism resonant both between laser and phonon andbetween the selected phonon and the spin sector. Whileanharmonic phonons have also been invoked for the cre-ation of effective static magnetic fields [73], our focus willbe on the coherent driving of harmonic phonon modes,which to date has been considered only theoretically, forclassical magnets with phenomenological damping [74].We consider the dynamics of the driven dissipativequantum spin system using the example of the alternat-ing spin chain shown in Fig. 1. The driving is effected bylaser excitation of an Einstein phonon that couples to oneof the magnetic interactions in the spin chain and the dis-sipation is modelled in the Lindblad formalism by bathoperators that damp both the lattice and spin sectorsdirectly. We establish the equations of motion govern-ing the basic physics of quantum NESS in this system,in terms of the driving frequency, the system parame-ters, and the response of the separate lattice and spinsectors. These equations enable us to discuss the dif-ferent regimes of weak and strong spin-phonon coupling,of weak and strong damping, and all the timescales as-sociated with driving, NESS formation, and relaxation. FIG. 1. Schematic representation of a lattice spin system, herea structurally dimerized chain with antiferromagnetic inter-action parameters
J > J (cid:48) , driven by the selective excitationof one specific phonon mode of the lattice. Both the drivenphonon and the spin system are damped by the ensemble oflattice phonons. We envisage an experimental geometry withthe sample attached to an efficient heat sink for thermal reg-ulation.
Thus our study establishes a foundation for many typesof extension, specifically to different types of spin system,to different types of bath (characterized by whether theyconserve spin and momentum), to finite system tempera-tures and thus to driving protocols for the management ofheat and of coherence, and to quantitative studies of ma-terials and device geometries for practical experiments.The structure of this article is as follows. In Sec. IIwe present our model for the quantum spin chain, forthe laser-coupled phonon mode that drives it, and forthe Lindblad bath operators that damp it. We derivethe equations of motion for the coupled lattice and mag-netic sectors and comment on their structure. SectionIII contains a preliminary analysis of the content of thesemaster equations, with specific attention paid to NESS.We demonstrate numerically that NESS can indeed beestablished, and illustrate how their basic properties aregoverned by the primary system parameters, namely thedriving power, the driving frequency, the lattice and spindamping coefficients, and the spin-lattice coupling. Withthis basis, in Sec. IV we concentrate on the regime of lowspin-phonon coupling to perform a complete investiga-tion of the dynamical properties of the NESS in the spinsector, characterizing their response by frequency, wave-vector components, and spin damping.In Sec. V we turn to a different but essential aspect ofNESS, namely the transient processes occurring as theyare established, from the moment the laser driving isswitched on, and the relaxation processes by which equi-librium is restored when the drive is removed. At highernet occupancies of lattice and spin excitations we find anomalously slow convergence to NESS, and in Sec. VAwe apply analytical arguments to discuss the underlyingphysics. Section VB extends this analysis to the questionof limits in parameter space for the existence of NESSwithin our model framework and Sec. VC provides a briefdiscussion of relaxation and temperature. In Sec. VI weanalyze the energy flow in the NESS, considering both itsuptake by the spin system as a function of laser powerand frequency and its dissipation by the Lindblad terms.This allows us to provide experimentally oriented esti-mates for the rate of temperature increase in the drivensystem, for its control by the heat sink shown in Fig. 1,and the resulting timescales for read-out and control pro-cesses. In Sec. VII we discuss the context of our resultsfrom a number of angles, including methodology, the in-fluence of the bath model, timescales and heating effects,and laser experiments on real materials. Section VIIIconsists of a brief summary and perspectives for futureextensions of the framework established in this study.
II. MODEL AND METHODS
We begin by representing the Hamiltonian of the cou-pled system shown in Fig. 1 as H = H s + H sp + H p + H l , (1)where the four terms describe respectively the spin sys-tem, the spin-phonon coupling, the Einstein phonon, andthe effect of the laser electric field on this phonon. Thebath operators damping the spin and lattice systems donot enter Eq. (1) explicitly, but are introduced at the levelof the Lindblad formalism. While some authors have in-vestigated the strong and controllable effects obtained byconsidering the quantum nature of the light field, gener-ally referred to as “cavity QED” [75, 76], for the purposeof driving phonon modes we treat the laser light field asclassical. A. Spin system
We express the Hamiltonian for the structurally dimer-ized, antiferromagnetic spin chain as H s = (cid:88) i J (cid:126)S ,i · (cid:126)S ,i + J (cid:48) (cid:126)S ,i · (cid:126)S ,i +1 , (2)with J > J (cid:48) >
0. For simplicity we consider onlyHeisenberg interactions between the spins and neglectany anisotropy terms; in real materials these could be ofexchange, Dzyaloshinskii-Moriya, single-ion, g -tensor, orother origin, and as noted in Sec. I are generally weak in3 d transition-metal compounds. A representation par-ticularly useful for dimerized spin systems is the bond-operator description [77, 78], in which the Hamiltonianis transformed by expressing the two spin operators oneach dimer using the identity S α , = ± ( s † t α + t † α s ) − i (cid:88) βζ (cid:15) αβζ t † β t ζ , (3)where s and t α ( α = x, y, z ) are operators for the sin-glet and triplet states of each dimer ( J ) bond. Theseoperators have bosonic statistics, required to reproducethe spin algebra of S α , ; however, because each dimermay only be in a singlet state or one of the three triplets(equivalent to the four possible states of two spin-1/2 en-tities), the bond operators must also obey a local hard-core constraint, s † i s i + (cid:88) α t † i,α t i,α = 1 , (4)on each dimer i , and hence are hard-core bosons.For an inversion-symmetric system, the minimalHamiltonian of Eq. (2) takes the form H s = H + H + H ,where [78, 79] H = (cid:88) i − J ( s † i s i − t † i,α t i,α ) (5) − µ i ( s † i s i + t † i,α t i,α − , with summation over the repeated index α , H = − J (cid:48) (cid:88) i,α (cid:0) t † i,α t i +1 ,α s † i +1 s i (6)+ t † i,α t † i +1 ,α s i s i +1 + H . c . (cid:1) , and H = J (cid:48) (cid:88) i,α (cid:54) = β (cid:0) t † i,α t † i +1 ,β t i +1 ,α t i,β (7) − t † i,α t † i +1 ,α t i +1 ,β t i,β + H . c . (cid:1) . The second term in H enforces the constraint [Eq. (4)]using the Lagrange multipliers µ i . At zero applied mag-netic field, the term quadratic in the singlet operators isnegative, which ensures a singlet condensation and jus-tifies their replacement by a constant, s i = (cid:104) s i (cid:105) , on eachdimer. The ground state of the system is then a conden-sate of singlets with a spin gap to all triplet excitations,whose dispersion is specified by the quadratic terms in H . Here we will not consider any spatial gradients (forexample in temperature, magnetic field, or laser flux) andhence (cid:104) s i (cid:105) = s and µ i = µ ; the latter condition enforcesthe hard-core constraint at a global level, but not locally.For the purposes of the present analysis we will not con-sider triplet-triplet interactions, and thus we neglect H .We transform the quadratic triplet Hamiltonian H + H to reciprocal space using t i,α = 1 √ N (cid:88) k t k,α e − ikr i , (8) where N is the number of dimers, and express the resultin the form H mf = E + (cid:88) k,α (cid:2)(cid:0) J − µ (cid:1) t † k,α t k,α − J (cid:48) s cos k (cid:0) t † k,α t k,α + t − k,α t †− k,α + t † k,α t †− k,α + t − k,α t k,α (cid:1)(cid:3) , (9)where E = N (cid:2) ( − J − µ ) s + µ − J (cid:3) . (10)We note that the only terms generated are those couplingoperators at wave vectors k and − k , and there is no mix-ing of the triplet indices, α . The conventional approach[78, 79] is to symmetrize the Hamiltonian matrix, diago-nalize it to obtain a new bosonic quasiparticle, known asthe triplon [80], form two mean-field equations, and solvethese for µ and s . By the use of effective quasiparticlestatistics, this procedure may also be followed at finitetemperatures [81, 82].Here we adopt one further simplification with a view toapplying equation-of-motion methods. In the “Holstein-Primakoff” approximation [79], the singlet occupationis replaced directly by invoking the local constraint[Eq. (4)], giving s = 1 − N (cid:88) k,α t † k,α t k,α . (11)At quadratic order, this substitution reduces to the ap-proximation s = 1, µ = − J , which is clearly valid inthe limit of a strongly dimerized chain. From extensivestudies of the spin ladder [82], it is generally recognizedthat the bond-operator description retains semiquanti-tative validity for interaction ratios J (cid:48) /J (cid:46) / H s = (cid:88) k,α (cid:2) Jt † k,α t k,α − J (cid:48) cos k (cid:0) t † k,α t k,α + t † k,α t †− k,α + t k,α t − k,α (cid:1)(cid:3) , (12)we diagonalize it by applying the Bogoliubov transfor-mation t k,α = ˜ t k,α cosh θ k + ˜ t †− k,α sinh θ k , (13a) t † k,α = ˜ t † k,α cosh θ k + t − k,α sinh θ k , (13b)where tanh 2 θ k = λ cos k − λ cos k (14)with λ = J (cid:48) /J , to obtain H s = (cid:88) k,α ω k ˜ t † k,α ˜ t k,α . (15)The operators ˜ t † k,α and ˜ t k,α create and destroy the triplonmodes of the dimerized chain and have dispersion relation ω k = J √ − λ cos k. (16) B. Phonon system and spin coupling
As shown in Fig. 1, we consider a situation in whichthe interaction J is modulated by the oscillations of onespecific phonon mode on every bond. We take this to bean Einstein phonon with wavevector q = 0 and a finiteenergy, ω . This optical phonon is driven by the lightfield of a laser that illuminates the entire sample, mean-ing that we take the driving to be a bulk effect. In a realmaterial, many different phonon modes are present in ad-dition to the driven phonon, and all of them, in particularthe acoustic phonons, are responsible for the dissipationof energy from both the lattice and spin sectors.The Hamiltonian terms involving the driven phononare H p + H sp + H l = (cid:88) j (cid:2) ω b † j b j + g ( b j + b † j ) (cid:126)S ,j · (cid:126)S ,j + E ( t )( b j + b † j ) (cid:3) , (17)where g is the spin-phonon coupling constant and E ( t ) = a cos( ωt ) is the oscillating electric field of the laser, whichwe assume to contain a single driving frequency, ω ; asnoted above, for the amplitudes a we consider, E ( t ) maysafely be treated as classical field. For our purposes, E ( t )is an internal field, meaning it is the fraction of the in-cident laser light transmitted into the sample, and wedo not concern ourselves with the reflected component.The dissipative terms do not enter Eq. (17), but will beincluded using the Lindblad formalism in Sec. IIC.The transformation of Eq. (17) includes single- andtriple-operator terms. For pedagogical accuracy we takea conventional definition of the Fourier transform, b j = 1 √ N (cid:88) q b q e − iqr j , (18)under which the electric-field term becomes E ( t ) √ N (cid:88) j,q ( b q e − iqr j + b † q e iqr j ) = E ( t ) √ N (cid:88) q ( b + b † ) , (19a)= N E ( t ) d, (19b)where only the q = 0 mode is selected, but we expressit as an intensive quantity summed over q , with effec-tive displacement operator d = √ N ( b + b † ). Even moresimply, the phonon term becomes (cid:80) q ω b † q b q .Finally, the spin-phonon coupling term becomes1 N √ N (cid:88) j,q,k,k (cid:48) ,α (cid:0) b q t † k,α t k (cid:48) ,α e i ( k − k (cid:48) − q ) r j + H.c. (cid:1) = 1 √ N (cid:88) q,k,α (cid:0) b q t † k,α t k − q,α + b † q t † k,α t k + q,α (cid:1) , (20)with q = 0 as the only relevant phonon mode. At the mean-field level one obtains the decoupled terms H sp = H sp,s + H sp,p , (21a) H sp,s = (cid:104) d (cid:105) (cid:104) (cid:88) k,α t † k,α t k,α − (cid:104) t † k,α t k,α (cid:105) eq (cid:105) , (21b) H sp,p = (cid:68) (cid:88) k,α t † k,α t k,α − (cid:104) t † k,α t k,α (cid:105) eq (cid:69) d, (21c)where we have omitted the product of the two expecta-tion values in Eq. (21a) because it has no influence atall on the dynamics of the system. Here H sp,s containsthe operator part acting on the spin degrees of freedom,expressed in triplon operators, while H sp,p contains theoperator part acting on the driven phonon. In both termsthe spin-phonon interaction is expressed by deducting theequilibrium value of the triplon occupation, such that ithas no effect when the system is not driven. While thismean-field decoupling is an approximation, we will showin Sec. VI that its quantitative limitations are minor.To transform H sp into the diagonal (triplon) basis ofthe spin sector, we use the identity t † k,α t k,α = y k (cid:0) ˜ t † k,α ˜ t k,α + (cid:1) − + y (cid:48) k (cid:0) ˜ t † k,α ˜ t †− k,α + ˜ t k,α ˜ t − k,α (cid:1) , (22)in which y k = 1 − λ cos k √ − λ cos k = J ω k /J ω k and (23a) y (cid:48) k = λ cos k √ − λ cos k = J − ω k /J ω k , (23b)to obtain the expression H sp,s = g (cid:104) d (cid:105) (cid:88) k,α (cid:2) y k [˜ t † k,α ˜ t k,α − n ( ω k )]+ y (cid:48) k (cid:0) ˜ t † k,α ˜ t †− k,α + ˜ t k,α ˜ t − k,α (cid:1)(cid:3) . (24)Here the bosonic occupation function, n ( ω k ) =[exp( (cid:126) ω k /k B T ) − − , provides an accurate value forthe equilibrium occupancy of the triplon mode with fre-quency ω k [Eq. (16)] at the low temperatures we consider,despite the hard-core nature of these modes [83].We define the operators u k = (cid:88) α ˜ t † k,α ˜ t k,α and (25a)˜ v k = (cid:88) α ˜ t † k,α ˜ t †− k,α (25b)in the triplon sector and denote their expectation valuesat any given time, t , by u k ( t ) = (cid:104) u k (cid:105) ( t ) , (26a)˜ v k ( t ) = (cid:104) ˜ v (cid:105) ( t ); (26b)the expectation value of the product of two annihilationoperators is manifestly the complex conjugate of ˜ v k ,˜ v ∗ k ( t ) = (cid:88) α (cid:104) ˜ t k,α ˜ t − k,α (cid:105) ( t ) . (27)We comment that the triplon branch, α , is summed overhere and, because we do not consider an applied magneticfield or any anisotropy in the spin Hamiltonian, will notenter our considerations again. u k ( t ) is clearly a real variable due to the hermiticity ofthe operator on the right-hand side of Eq. (25a), whilethe complex variable ˜ v k ( t ) is conveniently separated intoits real and imaginary parts, v k ( t ) = Re ˜ v k ( t ) (28a) w k ( t ) = Im ˜ v k ( t ) . (28b)In the equations of motion to be derived in Sec. IIC,the spin-phonon coupling introduces two quantities com-posed of the above expectation values, which we includeby defining U ( t ) = 1 N (cid:88) k y k [ u k ( t ) − n ( ω k )] , (29a) V ( t ) = 1 N (cid:88) k y (cid:48) k v k ( t ) , (29b)both of which are real by construction. For the descrip-tion of the spin sector in the driven system, we define thenumber, n x , of elementary (Bogoliubov, or “dressed”)triplons per site, n x ( t ) = 1 N (cid:88) k u k ( t ) , (30)and it will be helpful to compare this with the number oforiginal (or “bare”) triplons per site in the starting basisof Eq. (1), n b ( t ) = 1 N (cid:88) k,α (cid:104) t † k,α t k,α (cid:105) ( t ) . (31)Using Eq. (22), this last definition is equivalent to n b ( t ) = U ( t ) + V ( t ) + 1 N (cid:88) k,α (cid:104) t † k,α t k,α (cid:105) eq , (32)in which the last term is given by1 N (cid:88) k,α (cid:104) t † k,α t k,α (cid:105) eq = 32 N (cid:88) k [(2 n ( ω k ) + 1) y k − . (33)At zero temperature and for λ = 1 /
2, which will beour test case in what follows, the equilibrium expecta-tion value is n b0 = 0 . C. Equations of motion
The time evolution of an open quantum system is spec-ified by adjoint quantum master equations [45] of the form ddt A H ( t ) = i [ H, A H ( t )] + (34a) (cid:88) l ˜ γ l (cid:2) A † l A H ( t ) A l − A H ( t ) A † l A l − A † l A l A H ( t ) (cid:3) (34b)for any operator A H ( t ) describing a physical observable.In these Heisenberg equations of motion, H is the Hamil-tonian of the “reduced” system under consideration, bywhich is meant the quantum system with no environ-ment. The “Lindblad” operators, { A l } , are formed fromthe Liouville space of the reduced system to describe itsinteraction with the environment (the “bath”), which isexcluded from explicit consideration. It was proven byLindblad [44] that Eq. (34) is the most general form of thedissipation term for a separable (system-bath) Hilbertspace when l describes a bounded set of operators. Thecoefficients ˜ γ l play the role of damping parameters.The driven phonon exemplifies the textbook case[45, 46] of the Lindblad equations, namely those of thedamped harmonic oscillator. The Lindblad operators inthis case are A = b † and A = b , with damping rates γ and γ . For the system to relax back to its equilib-rium state in the absence of driving, it is known [45] thatthe ratio of the two rates must be given by the ratio n ( ω ) / (1 + n ( ω )), and hence the conventional parame-terization is ˜ γ = γ n ( ω ) (35a)˜ γ = γ (1 + n ( ω )) , (35b)leaving only one damping parameter, γ . For physicaltransparency we separate the Lindblad operators intothose that excite the system by an energy ω l , which wedenote by B l , and those that de-excite, which are givenby the Hermitian conjugates, B † l . The dissipative part ofEq. (34), which is the second line, may then be separatedinto the two contributions T = (cid:88) l γ l n ( ω l ) (cid:8)(cid:2) B l , [ A H ( t ) , B † l ] (cid:3) + (cid:2) B † l , [ A H ( t ) , B l ] (cid:3)(cid:9) , (36a) T = (cid:88) l γ l (cid:8) [ B l , A H ( t )] B † l + B l [ A H ( t ) , B † l ] (cid:9) . (36b)The commutators in these expressions facilitate theirrapid evaluation in comparison with the expression inEq. (34); if the observable and the Lindblad operatorsare linear bosonic operators, as for the damped phonon,it can be seen without explicit calculation that the term T vanishes and hence no dependence on the bosonic oc-cupation, n ( ω l ), arises.To describe the driven phonon we consider the realvariables q ( t ) = (cid:10) √ N ( b + b † ) (cid:11) ( t ) , (37a) p ( t ) = (cid:10) i √ N ( b † − b ) (cid:11) ( t ) , (37b) n ph ( t ) = (cid:10) N b † b (cid:11) ( t ) , (37c)describing respectively the displacement of the Einsteinphonon ( d in Sec. IIB), the conjugate phonon momen-tum, and the number operator. We recall that, despitethe presence of all phonon modes, b † q , in H p , only theoperators b and b † appear elsewhere in the Hamilto-nian of the reduced system, and hence are candidates forthe formation of Lindblad operators. The other phononsform the environment and their presence gives rise to thedamping, which is contained in the single parameter γ .By evaluating Eqs. (34) or (36) with the expectationvalues from the spin sector [Eqs. (29)] in the spin-phononcoupling term, one obtains the closed set of equations ofmotion [45] ddt q ( t ) = ω p ( t ) − γq ( t ) , (38a) ddt p ( t ) = − ω q ( t ) − γp ( t ) − E ( t ) + g ( U ( t ) + V ( t ))] , (38b) ddt n ph ( t ) = − [ E ( t ) + g ( U ( t ) + V ( t ))] p ( t ) − γ [ n ph ( t ) − n ( ω )] . (38c)One observes the characteristic structure in Eqs. (38a)and (38b) of the displacement and momentum servingas conjugate time derivatives, but damping themselvesthrough the γ/ n ( ω ) is theoccupation of phonon mode ω at thermal equilibrium.Here we do not extend these considerations to bilinearLindblad operators in either the phonon or the spin sec-tor.Turning to the spin degrees of freedom, we considerthe real expectation values u k ( t ), v k ( t ), and w k ( t ) intro-duced in Eqs. (26a) and (28) to describe the dynamicalspin processes diagonal and off-diagonal in the triplonnumber basis. To determine the equations of motion, itis expected that the spin sector will be subject to a di-rect damping due both to weak spin-anisotropic termsand to phononic processes arising from the many acous-tic and optical phonon modes in the Hamiltonian of anyreal material. The specific nature of these damping pro-cesses will be the subject of more extended discussionin Secs. IV and VII, but the available Lindblad opera-tors will in general be linear and bilinear combinationsof ˜ t k and ˜ t † k . In the present analysis we focus on lin-ear operators, in order to present the primary phenom-ena associated with the driven dissipative quantum spinchain. The equations of motion we will deduce have theanalytical advantage of maintaining a simple form withtransparent physical consequences. However, it is truethat such one-triplon Lindblad operators are spin non-conserving, meaning that this type of bath is appropri-ate for materials with the non-negligible spin anisotropiesmore commonly associated with systems of 4 d and 5 d magnetic ions. Usually such anisotropic terms are nev-ertheless corrections to the spin Hamiltonian of Eq. (2),whereas they may be the leading dissipative terms; wediscuss this situation in more detail, and comment onthe case of spin-conserving bath operators, in Sec. VII.Thus the Lindblad operators, B k , that we consider aresimply ˜ t † k,α , and have damping coefficient˜ γ k = γ s n ( ω k ) , (39)while B † k has damping γ s (1 + n ( ω k )). We neglect a pos-sible dependence of γ s on the wave vector, k , along thechains. While one may ask whether this approximationconstitutes a severe omission, given that energy and mo-mentum conservation allow dissipation only for partic-ular combinations of both, we observe that energy con-servation as contained in Fermi’s Golden Rule does notimpose a strong constraint when one recalls that the 1Dchains are embedded in a three-dimensional (3D) crys-tal. Thus k -dependent damping coefficients, γ s ( k ), areaveraged over the transverse momentum, (cid:126)k ⊥ , and the as-sumption that energy conservation is satisfied at somevalue of (cid:126)k ⊥ is fully justified. While some dependence of γ s on the longitudinal momentum, k , may indeed remain,we proceed for the purposes of our present pedagogicalexposition with a single value of γ s for clarity.To deduce the equations of motion when A H ( t ) inEq. (34) is one of the bilinear operator combinations inEq. (25), we first consider the Hamiltonian parts of therespective expressions,[ H s , u k ] = 0 , (40a)[ H s , ˜ v k ] = 2 ω k ˜ v k , (40b)[ H sp,s , u k ] = gq ( t ) y (cid:48) k (cid:0) ˜ v † k − ˜ v k (cid:1) , (40c)[ H sp,s , ˜ v k ] = 2 gq ( t ) (cid:2) y k ˜ v k + y (cid:48) k (cid:0) u k + (cid:1)(cid:3) . (40d)Combining the unitary parts of Eqs. (40) with the dis-sipative part, T , from Eq. (36b), and taking the appro-priate expectation values, leads to the final expressions ddt u k ( t ) = 2 gq ( t ) y (cid:48) k w k ( t ) − γ s [ u k ( t ) − n ( ω k )] (41a) ddt v k ( t ) = − ω k + gy k q ( t )] w k ( t ) − γ s v k ( t ) (41b) ddt w k ( t ) = 2[ ω k + gy k q ( t )] v k ( t )+2 gq ( t ) y (cid:48) k (cid:2) u k ( t ) + (cid:3) − γ s w k ( t ) . (41c)In combination with Eqs. (38a)-(38c), these form theequations of motion for the coupled spin-lattice system.Regarding the structure of these equations, we commentonly that n ph ( t ) [Eq. (38c)] does not have any direct ef-fect on the evolution of the other coupled equations andhence it appears that this variable can be neglected fordynamical purposes, but we will continue to show n ph ( t )as a valuable diagnostic of the state of the driven phononsector.Regarding the solution of these equations, in order tostudy the steady-state and dynamical properties of thedriven and dissipative ensemble of Fig. 1, this will be ourtask in Secs. III and IV. In the majority of our calcula-tions, we will use a periodic chain of N = 400 dimers andhence by inversion symmetry will have 201 independentvalues of k , which we will consider both separately andin summed quantities such as Eqs. (29) and (30). Theequations of motion [Eqs. (38a)-(38c) and Eqs. (41a)-(41c)] have no lower or upper validity cutoff in time, andthus can be applied to discuss the formation, switching,and relaxation of quantum spin NESS from t = 0 to ∞ . III. NESS IN THE PHONON-DRIVEN SPINSYSTEM
We begin by choosing input parameters that establishquantum spin NESS, deferring a detailed analysis of thelimits to NESS formation until Sec. V. Our first aim is apreliminary characterization of the response of NESS tothe different factors influencing their driving. To reducethe space of possible driving parameters, in the presentanalysis we restrict our considerations to resonant exci-tation of the Einstein phonon mode, meaning that weselect the laser frequency such that ω = ω and hence E ( t ) = a cos( ω t ). From a driving standpoint, for theelectric-field intensities we wish to study and for a gener-ically weak spin-phonon coupling, off-resonant driving islargely just a less efficient means, by a factor propor-tional to [( ω − ω ) + ( γ/ ] − , of exciting a responseat frequency ω . However, in systems with stronger spin-phonon coupling, nontrivial phenomena are indeed foundby pumping and probing at frequencies ω (cid:54) = ω . We re-mind the reader that the minimal model of Sec. IIA wasnot designed to describe driving by any of the other phys-ical mechanisms summarized in Sec. I, all of which areless frequency-selective than phonon driving. It is easy toanticipate that the strongest effects of the driven phononon the spin system will be found when ω matches thespectrum of triplon excitations.We consider first the driven phonon system withoutcoupling to the spin chain, meaning with g = 0. Torepresent the phonons of a typical inorganic materialwe choose a damping coefficient γ = 0 . ω . FromEqs. (38a)–(38c) one observes that, up to a couplingto the spin system ( g ) that is typically below 10%, thephonon has the behavior of a classical damped harmonicoscillator with a driving term. This is borne out by thetime-dependence of the variables q , p , and n ph , shown inFig. 2. Figure 2(a) illustrates that the phonon number isdriven up to a finite average value, and the inset that itoscillates steadily around this constant value for all latertimes; this is the NESS of the driven phonon system.Figures 2(b) and 2(c) show the corresponding behaviorof the displacement and momentum, which have a rela-tive π/ FIG. 2. Response of the Einstein phonon to a resonant drivingfield. Here ω /J = 1, a = 0 . J , γ = 0 . ω , and g =0. (a) Phonon number, n ph ( t ), produced by switching on aconstant laser electric field at t = 0. The inset shows thesteady state of the driven phonon system at long times. (b)Phonon displacement, q ( t ), and momentum, p ( t ), shown from t = 0. (c) q ( t ) and p ( t ) at long times. the phonon number operator in this laser-pumped steadystate has been driven to a non-equilibrium average valueof approximately 0.04. Although this value appearssmall, it does constitute a macroscopic occupation of asingle mode. This driven ω phonon is the primary sourceof lattice excitations in the system, and all other phononmodes will have very low occupations at low tempera-tures. In all of the considerations to follow, we maintainthe value of n ph in this range, both for meaningful com-parisons as other parameters are varied and for a realisticaccount of the temperature of the steadily driven system,as we discuss in Sec. VI.Second, the frequency of the oscillations in the drivenphonon occupation, n ph ( t ), is twice that of q ( t ) and p ( t ),as expected from the number of nodes in the displace-ment cycle; the latter pair can be taken as the base fre-quency of the system, while the former is characteristic of FIG. 3. Average value of the driven phonon occupationnumber, n ph0 , in a NESS of the phonon system, displayedas a function of a for various driving frequencies at fixed γ = 0 . ω and γ s = 0 . J . (a) No coupling to the spinsystem ( g = 0). (b) g = 0 . J . ω , reflecting the fact that n ph is essentially the sum ofthe squares of q and p . Third, the characteristic timescalefor convergence of the average of n ph to the phonon NESSis 2 /γ for all three quantities [Figs. 2(a) and (b)]. For q ( t ) and p ( t ), this is to be expected from the correspond-ing equations of motion [Eqs. (38)], which contain ex-plicit terms with prefactor − γ/
2, while for n ph ( t ) it isthe behavior of p ( t ) on the right-hand side of Eq. (38c)that induces the same convergence rate. Because theconvergence is exponential, the actual establishment of aphonon NESS depends on the chosen accuracy criterion.To a good approximation, the phonon number in theNESS [inset, Fig. 2(a)] is given by n ph ( t ) = n ph0 + n ph2 cos(2 ω t ); we will investigate the corrections to thissituation, which arise due to coupling to the spin system,in Sec. IV. To study the quasi-stationary behavior of theNESS, we focus on the mean phonon occupation, n ph0 .Figure 3(a) shows that, for all driving frequencies, theaverage energy in the driven phonon mode rises with thedriving power, which is proportional to the square of the electric-field amplitude. In a fully classical treatment ofthe driven oscillator, n ph0 ∝ ( a/γ ) , and this result mayequally be understood from Fermi’s Golden Rule, wherethe flow of energy into the system is proportional to thesquare of the matrix element, and hence to a . We dis-cuss the topic of energy flow in detail in Sec. VI. Becausewe have chosen for realism to scale γ to the phonon fre-quency, a will also be scaled to ω in all of the studies tofollow, thereby maintaining a constant ( a/γ ) when ω is varied.The spin system is driven by the pumped phononthrough the coupling parameter g . Given that the am-plitude of the phonon oscillation, q ( t ), is proportional to a/γ , it follows that the amplitude of the induced driv-ing of the spin system is proportional to ga/γ . Figure 3compares the driven phonon system with g = 0 to thesituation with a finite value of g . Here we have chosendriving parameters suitable for the formation of NESS;those causing the spin system to inhibit NESS formationare the explicit focus of Sec. V. We observe in Fig. 3(b)that a generic spin-phonon coupling, g = 0 . J , results inonly small changes being induced by the spin system rel-ative to the isolated driven and damped phonons of mostfrequencies, but that some more significant alterationsare possible at specific phonon frequencies, for reasonswe investigate next.Turning now to the response of the driven spin system,it is necessary first to establish the nature and character-istic frequencies of the excitations created by the drivingphonon. Throughout the present study, we will considerthe dimerized S = 1 / λ = J (cid:48) /J = 1 /
2. Equation (16) statesthat the triplon modes of this chain form one triply de-generate branch dispersing from a value of ω min = J/ √ k = 0 to ω max = (cid:112) / J at k = π . However, by spinconservation it is not possible for a phonon to create a sin-gle spin excitation, and from the form of H sp in Eq. (17)it is evident that one phonon ( b † ) couples to two spinexcitations. One therefore anticipates that the strongesteffects of the driving phonon on the spin system will befound when ω is chosen to lie within the band of two-triplon excitations, namely when 2 ω min ≤ ω ≤ ω max (1 . J ≤ ω ≤ . J ). Thus an origin for the specialbehavior of the ω /J = 1 . g = 0 . J . We include a direct spin damping, γ s = 0 . J , which we scale to the energy of the spin sys-tem; to reflect the observed fact that the spin degreesof freedom are in general very weakly damped, we alsoadopt a value that is significantly lower than the phonondamping over most of the range of ω . Figure 4(a) showsthat laser driving at any frequency does create a responsein the spin system that is qualitatively similar to that in0 FIG. 4. (a) Response of the spin system, measured by n x ( t )[Eq. (30)], to driving phonon frequencies ω /J = 0 .
5, 1.0, 1.5,2.0, 2.5, and 3.0. The driving field ensures that a/γ = 0 . g = 0 . J , γ s = 0 . J . n x ( t ) in the spin NESS is shownat (b) ω /J = 0 .
5, (c) ω /J = 1 .
5, and (d) ω /J = 2 .
5, wherewe compare results in the time window 1160 ≤ t ≤ ≤ t ≤ the phonon system, namely that the spin occupation is“pumped” to a new average value, about which it oscil-lates. At constant ( a/γ ), the average triplon occupation, n x0 , displays a hierarchy of values as the NESS is ap-proached. While at frequencies far from the two-triplonband ( ω /J = 0 .
5, 1.0, and 3.0) this degree of drivingproduces only a very weak occupation, n x0 < .
001 [in-set, Fig. 4(a)], for frequencies in or near the band we finda state with n x0 (cid:39) .
05 at ω /J = 1 .
5, but also one with an occupation of only n x0 (cid:39) . ω /J = 2 . n x ( t ), is shown for three se-lected frequencies in Figs. 4(b) to 4(d). In each case wecompare the triplon occupation in a time window nearthe center of Fig. 4(a) with the long-time limit, for whichwe take the window 9960 ≤ t ≤ ≤ φ < π to start eachcycle at the same point. The most important result ofFigs. 4(b) to 4(d) is to prove that the driven model systemdamped by γ and γ s does indeed host spin NESS, in thatidentical periodic traces are obtained for arbitrarily longtimes. The subsidiary result is that, for most ω values,a good approximation to the NESS is reached alreadyat rather short times. Because convergence is exponen-tial, any meaningful accuracy criterion will be reachedafter a single-digit number of time constants, and thusfor quantitative purposes (Sec. VI), bearing experimentaluncertainties in mind, we define a NESS to exist using arelative criterion of 2% (corresponding to approximately4 /γ s ). According to this criterion, the driven state inFig. 4(c) is not yet a NESS, for reasons we will revisitbelow, but those shown in Figs. 4(b) and 4(d) are.As a benchmark for the meaning of the n x0 values inFig. 4, one may compare with the value n b0 = 0 . n x0 , one may state that the spinNESS established at low and high frequencies constituteonly a weak perturbation of the equilibrium state. Thisresult also implies that in the “Floquet” regime of fre-quencies above the two-triplon band, the spin state isnot altered qualitatively, although it may obtain a non-trivial phase structure. By contrast, for some frequen-cies in and around the two-triplon band, the quantumspin NESS can be altered significantly from the equilib-rium state, and our results for ω /J = 1 . n x = 1 on thetriplon occupation, and in fact such a situation wouldrepresent the most extreme out-of-equilibrium state pos-sible, at which many of the approximations in Sec. IIAwould no longer be valid. Anticipating the discussion ofSecs. IV and VA, we introduce an operational thresholdvalue of n x ( t ) in the driven spin state, such that our de-scription of the spin sector will remain appropriate, andwe set this to n maxx = 0 . n x0 as a consequence of the driving frequency,Fig. 4 invites two further remarks. First, we observe that1the time structure, n x ( t ), of the NESS in Figs. 4(b) to4(d) shows a rather complex form, with a definite su-perposition of different frequency harmonics in evidence.We will investigate this harmonic mixing, which appearsto be strongest at the below-band frequency of Fig. 4(b),in detail in Sec. IV. Second, the timescale over whichthe spin system reaches its NESS appears to be simi-lar at all frequencies, other than ω /J = 0 . t ≈ J − . This value corresponds to 4-5 time constantsof the spin system (1 /γ s ). Of the exceptional cases, at ω /J = 0 .
5, where γ s = γ , the process is somewhat de-layed by the phonon “switch-on” timescale (Fig. 2). At ω /J = 1 .
5, the process appears to be longer still, withthe NESS not yet fully established after t = 1200 J − [Fig. 4(c)]. We will investigate the transient behavior ofthe spin system at switch-on, and explain this curiouslyslow convergence, in Sec. VA.We conclude our initial survey of spin NESS in re-sponse to a driving phonon by showing the spin-systemanalog of Fig. 3. In the analysis of experiments, a keyquantity in characterizing any phenomenon is its depen-dence on the power, or fluence, of the laser, which isquite straightforward to measure. From elementary elec-trodynamics, the fluence is proportional to the squaredamplitude of the laser field, and hence in Fig. 5(a) weshow the dependence on a of the average triplon occu-pation, n x0 , in the NESS for the six representative driv-ing frequencies. As for the driven phonon, the depen-dence is clearly linear over the full range of γ -normalized a values for all driving frequencies, again except for ω /J = 0 . ω /J = 1 .
5. The latter shows a satura-tion as n x0 is driven towards unphysical values at verylarge a , while the former shows a crossover to a depen-dence that it as least quadratic in ( a/γ ) at strong driv-ing. Next (Sec. IV) we discuss the dynamical propertiesof the driven NESS, which will allow us to understandthe origin of this form, after which (Sec. V) we will ad-dress the issue of limits on ( a/γ ) for spin NESS to existat long driving times.In Fig. 5(b) we show the dependence of the driventriplon occupation on the spin-phonon coupling, g , forthe same six driving frequencies. At low values of g , n x0 shows a g form that is directly analogous to its depen-dence on a . However, at high g we observe a suppressionof n x0 below its expected value, the onset of which oc-curs at lower g for the phonons closest to resonance withthe two-triplon band, and find that the spin response caneven decrease as the coupling is increased. This onset ofmore complex behavior, which is also evident in the re-sponse of the ω /J = 1 . g = 0 . J , and a“strong-coupling” regime. Most magnetic quantum ma-terials do not show strong spin-phonon coupling at equi-librium, and thus for the purposes of the present analy-sis, which is to discuss the properties of a generic drivenquantum magnet, we will focus on the weak-couplingregime. Hence we adopt the value g = 0 . J to be rep- FIG. 5. (a) Dependence of the average triplon number, n x0 , inthe NESS on the fluence, shown as ( a/γ ) , at driving phononfrequencies ω /J = 0 .
5, 1.0, 1.5, 2.0, 2.5, and 3.0. The fixedsystem parameters are g = 0 . J , γ = 0 . ω , and γ s = 0 . J .Only the ω /J = 0 . n x0 onthe spin-phonon coupling constant, g , for driving phonons ofthe same six frequencies at fixed a/γ = 0 .
2. A well-defined g dependence at all small couplings gives way to a suppressionof n x0 at larger g values whose onset depends on ω . resentative of the class of magnetic materials in which toseek linear quantum spin NESS phenomena.In this weak-coupling regime, one may exploit theequivalence of a and g to define a dimensionless effec-tive driving parameter for the spin system, D = ga/ ( γJ ) , (42)which can be used to simplify the analysis, and we willemploy this parameterization in Sec. V. However, whenworking beyond this regime it is not possible to avoidstudying the full space of a/γ and g . Although we deferthe analysis of strong coupling to a later study, we stressthat all of the treatment in Sec. II remains fully validfor all the g values shown in Fig. 5(b). Nevertheless, as2we will mention in Sec. VII, values of g up to 0 . J areknown in some dimerized-chain compounds, and for suchextreme spin-phonon coupling one may not exclude thepossibility of a different type of physics at equilibrium,such as the formation of combined phonon-triplon enti-ties; we comment only that the formalism of Sec. II wouldnot be appropriate for such a situation. IV. DYNAMICAL PROPERTIES OF THEQUANTUM SPIN NESS
We turn now to a quantitative analysis of the dynam-ics of the spin NESS. It is already clear from Sec. III, andparticularly Fig. 4, that the superposition of frequenciespresent in the steady state can be complex. For full in-sight into the harmonic content of the spin NESS, we in-troduce the Fourier transform (FT) of the NESS signal,which we apply to n ph ( t ), to the individual spin com-ponents, u k ( t ) and v k ( t ), and to the summed quantities n x ( t ) [Eq. (30)] and V ( t ) = 1 N (cid:88) k v k ( t ) , (43)which characterize respectively the average triplon oc-cupation and the average behavior in the off-diagonaltwo-triplon sector. The definition of the FT is simpliedby making use of the results in shown in Fig. 4, wherewe demonstrated that NESS had been achieved at longtimes. We use one cycle of the signal taken from thetime window 9960 ≤ t/J − ≤ X ( t ) = (cid:88) m X m exp( imωt ) (44)for any quantity X appearing in a NESS driven by anyfrequency ω ; with this notation, any quantity with an in-teger subscript ( X m ) denotes a Fourier component, andthose with m = 0 are all real numbers. Without per-forming a detailed analysis beyond the level of Fig. 4,we comment that the system described by the model ofSec. II does not generate any significant dynamics at fre-quencies other than mω , where m is an integer. We alsocomment that there are no discernible extrinsic featuresarising in the FT as a consequence of the finite length ofthe chain on which we perform our calculations.Returning to the case of resonant driving ( ω = ω ),we illustrate the FT in Fig. 6 by showing in the leftpanels the time structure of n ph ( t ), n x ( t ), V ( t ), and thesingle- k components u k =0 ( t ) and u k = π ( t ) in the NESS ofFig. 4 at ω /J = 1 .
5; juxtaposed in the right panels arethe corresponding harmonic decompositions determinedfrom Eq. (44). We have chosen a relatively conventionalNESS trace [Fig. 6(c), similar to Fig. 4(c)], in which n x ( t )and u k ( t ) are dominated by the even Fourier components m = 0 and 2, while V ( t ) [and by extension v k ( t )] is dom-inated by m = 1. This result is quite natural if one FIG. 6. Illustration of the quantities (a) n ph ( t ), (c) n x ( t ), (e) V ( t ), (g) u k =0 ( t ), and (i) u k = π ( t ) in the NESS obtained withdriving parameters a/γ = 0 . g = 0 . J at frequency ω /J = 1 . γ s = 0 . J . Pan-els (b), (d), (f), (h), and (j) show the corresponding Fourierdecompositions. considers the equations of motion [Eqs. (41)], taking q ( t )to be a sinusoidal driving with small amplitude and afrequency ω . At leading order in q ( t ), all components v k ( t ) and w k ( t ) will also oscillate at this same frequency,giving a dominant m = 1 component, while the leading-order response in u k ( t ) oscillates at 2 ω and possesses aconstant offset (a zeroth harmonic). Because n x ( t ) is thesum over all u k ( t ), it therefore shows harmonic compo-nents primarily at m = 0 and 2. All of these featuresare evident in Figs. 6(c)-6(j). In addition, we observethat the different k -components of u k ( t ) display differentharmonic contributions, and because ω /J = 1 . m = 0 and 2coefficients are larger at k = 0 than at k = π ; one mayverify (data not shown) that the converse is true at adriving frequency of ω /J = 2 .
5, and we consider the k -dependence of the response in more detail below.Given this conventional behavior of the spin NESS,it is somewhat surprising to observe the presence of asignificant m = 1 harmonic in the phonon NESS, n ph ( t ),of Fig. 6(a). In fact n ph0 is suppressed by 14% comparedto its g = 0 value (Fig. 2), which is a weaker version of theeffect visible for the ω /J = 1 . m = 1 harmonic is another consequence ofstrong feedback from the spin system at this “resonant”(in-band) frequency, and arises from the term g U ( t ) p ( t )in Eq. (38c), where U ( t ) oscillates primarily at 2 ω and p ( t ) at ω . It is also clear that additional harmonics arepresent in the spin NESS analyzed in Fig. 6, including3 FIG. 7. Coefficients of the Fourier transforms of (a) n x ( t )and (b) V ( t ) in the NESS obtained with driving a/γ = 0 . g = 0 . J , shown as a function of the driving phononfrequency, ω , for damping parameters γ = 0 . ω and γ s =0 . J . at higher multiples of ω , and one may anticipate [notleast from Fig. 4(b)] that for certain frequencies they aresignificant.To investigate the effect of the frequency of the driv-ing phonon, in Fig. 7 we show the coefficients of n x ( t )and V ( t ) from m = 0 to 4 as a function of ω . Acrossthe full range of frequencies, n x ( t ) is indeed dominatedby the m = 0 and 2 coefficients [Fig. 7(a)] and V ( t ) by m = 1 [Fig. 7(b)], meaning that the case study of Fig. 6,performed for ω /J = 1 .
5, is in fact well representativeof the hierarchy of coefficient values, with only one sig-nificant exception. This is the frequency range around ω = ω min , where a clear peak appears in a number ofthe harmonic components. Although frequencies around ω /J = 0 . ω /J = 1 . J ) coincides with the peak density of states at the two-triplon band minimum. In-spection of Eqs. (41b) and (41c) reveals that oscillationsare indeed induced at 2 ω because q ( t ) is multiplied by v k ( t ) or w k ( t ). While this process appears at next-to-leading order in q ( t ), it is strongly enhanced when thesecond harmonic satisfies the resonance condition.The resonantly enhanced second harmonic of w k ( t ) inturn induces stronger first and third harmonic compo-nents in u k ( t ), as may be read from Eq. (41b), in which w k ( t ) is multiplied by q ( t ) where it acts as a driving termfor u k ( t ) at frequencies (2 ± ω . This type of harmonicmixing results, for the driving we consider in Fig. 7, in thecoefficient | n x1 | even exceeding | n x2 | around ω = ω min ,where | n x3 | is also strongly enhanced. Similarly, | V | and | V | are also enhanced over a wide frequency rangearound ω = ω min , where at its peak | V | approaches | V | .Thus the resonant enhancement of the second harmonicexplains why the temporal behavior of the spin NESSdisplays more and different features at below-band fre-quencies around ω /J = 0 . ω /J = 1 . ω /J = 2 . ω = ω min . First,it is important to stress that the response observed at2 ω in n x ( t ) is not a doubling phenomenon; it is merelya consequence of the fact that the triplon number is anoperator square of the triplon degree of freedom, andin this sense the behavior of u k ( t ), v k ( t ), and w k ( t ) isdirectly analogous to that of the driven phonon vari-ables discussed in Sec. III. By contrast, the frequency-doubling observed between phonon driving at ω = ω min and the strong response of the spin system at 2 ω min isa real effect, which at a “classical” level can be read di-rectly from the equations of motion. At a quantum level,this frequency-doubling requires the involvement of twophonons at frequency ω , taking part in off-shell phonon-triplon processes that are allowed in the strongly out-of-equilibrium system.We stress again that all physical processes of this type[meaning those contained in Eqs. (41)] do involve multi-ple driving phonons, as is standard in Floquet physics.Our treatment of the lattice system does not allow forthe creation of phonons with frequencies of 2 ω , 3 ω , orhigher due to anharmonicities in the lattice potential,as was discussed in Refs. [24, 25]. Because the factorsenhancing multi-phonon response and harmonic mixing(Figs. 6 and 7) are the same, it is no surprise to find thatboth phenomena are strongest in the same range of fre-quencies. Quantitatively, the strength of these subdom-inant signals at constant a/γ is a product of powers of g with the height of the density-of-states peak at 2 ω min ,and the enhancement can exceed an order of magnitudeat ω = ω min .Turning to the physical quantities characterizing theNESS, we have seen in Sec. III, and see again in Fig. 7,that the response of the spin system is very sensitively4 FIG. 8. (a) Average triplon occupation, n x0 , in the NESSobtained with driving a/γ = 0 . g = 0 .
05, shown onlogarithmic axes as a function of ω for different values of γ s . The band-edge features become increasingly prominentas γ s decreases, as does the peak at ω = ω min , but for mostother phonon frequencies n x0 is quite insensitive to the spindamping. (b) Corresponding off-diagonal response, shown bythe quantity | V | . dependent on the driving frequency, with clearly differ-ent adiabatic, antiadiabatic, and “resonant” (by whichis meant in-band) forms. However, some in-band fre-quencies are not particularly remarkable, due to smallmatrix elements or low densities of two-triplon states,and some adiabatic frequencies clearly have rather stronganomalous (multiphonon) enhancement. For a quantita-tive visualization of this response, in Fig. 8(a) we showthe mean amplitude, n x0 , of the driven triplon occupa-tion and in Fig. 8(b) the amplitude of the off-diagonalresponse, which we gauge using | V | . The rising linesindicate decreasing values of γ s , which we terminate at γ s = 0 . J to avoid having n x0 exceed n maxx = 0 . n x0 is surprisingly FIG. 9. Wave-vector-resolved average triplon occupation, u k ,in the NESS obtained with driving a/γ = 0 . g = 0 . J ,shown as a function of the driving frequency, ω , for k = 0, π/ π/ (cid:15) , 3 π/
4, and π . Black squares show the maxima, u max k , of the u k functions peaking at different energies acrossthe Brillouin zone, which defines the wave vectors k res . At allfrequencies ω < ω min the strongest peak is found in u k =0 and at ω > ω max in u k = π . (cid:15) = π/N is an offset from theband center, where u k = π/ = 0. insensitive to γ s [Fig. 8(a)]. However, as ω approaches2 ω min and 2 ω max , the driven n x0 varies strongly with γ s ,and the same is true around ω = ω min .As Fig. 8(b) makes clear, analogous effects are presentat ω = 2 ω min and 2 ω max in | V | , which also rises tovalues of order 0.1 at the lower band edge for γ s = 0 . J but is essentially independent of γ s for driving frequenciesmore than 0 . J outside the two-triplon band. Because wehave chosen | V | as the off-diagonal diagnostic, and this isa primary driving term in Eqs. (41) rather than a driventerm, there is no γ s -dependence around ω = ω min ; thisresponse of the off-diagonal sector is found rather in thecoefficients | V | and | V | in Fig. 7. The other differencesin the frequencies of characteristic features in | V | , mostnotably the in-band minimum occurring at ω /J = 1 . y k in Eq. (23a) as opposedto y (cid:48) k in Eq. (23b).To understand the degree to which individual k -components of the spin system are selected by the phonondriving, in Fig. 9 we show u k over the full range of driv-ing frequencies for selected values of k across the Brillouinzone. For k = 0 and π , it is no surprise that the re-spective u k functions peak strongly at ω = 2 ω min and2 ω max , because these are the dominant available wave-vector components; we note that there is no problemwith the fact that u k =0 ( t ) exceeds the threshold whenthe system is driven at ω = 2 ω min , because the triplonoccupation is determined by the average over all compo-nents [Eq. (30)]. For driving frequencies within the two-triplon band, one might expect a broad spin response on5the grounds that triplon pairs from a wide range of wavevectors may contribute. However, the response at eachfrequency remains dominated by the resonance condition ω = 2 ω k , and thus the components u k for k = π/ π/ k -selection on the basis of the drivingenergy is rather accurate and it is well justified to intro-duce a “resonant” wave vector, k res , selected by each ω .The black squares in Fig. 9 show the maxima, u max k , of asequence of u k res ( t ) functions selected in this way.In addition to this characteristic frequency, each u k shows a pronounced below-band two-phonon process, vis-ible at one half of the peak frequency, and it is only theact of averaging over all the k -components that disguisesthese features in our figures showing n x0 . For the drivingand damping parameters used in Fig. 9, no three-phononprocesses are discernible in the individual k -components.Nevertheless, a wealth of structure is revealed by consid-ering the FT of the different k -components on logarith-mic axes for a range of frequencies (analogies of Fig. 6,data not shown). The differential response of different k -components is also clearly visible when the drive isswitched on, leading to complex envelope oscillations atinitial times, and we will touch on these phenomena inSec. VA. We remind the reader that the structure of ourmodel ensures no interactions between triplons at differ-ent k , and so all u k ( t ) components evolve independentlyin time.We close our discussion of dynamical phenomena in theNESS by commenting on the possibility of new dynamicalmodes emerging in the driven system, for example wherethe pumped phonon is strongly dressed by triplons. Ex-citations with combined phononic and spin character areknown in a number of materials, including manganitesand “spin-Peierls” chains. In general these are a prop-erty of the equilibrium system arising for strong g and, asnoted at the end of Sec. III, their inclusion would requirean extension of the present treatment. While this treat-ment does reveal unconventional dynamical processes inthe driven system, specifically those involving multiplephonons, it is not designed to capture the formation ofbound states of these excitations at equilibrium. V. TRANSIENT AND RELAXATIONPROCESSES
Although the primary aim of our present study is todiscuss NESS themselves, clearly their short-time (tran-sient) behavior on “start-up” is a key to measurementwindows, as well as to analyzing switching processes ofthe type one may wish to use in logic operations. Despitethe clear presence of the timescales set by the lattice andspin dampings, respectively 1 /γ and 1 /γ s , we have al-ready observed in Figs. 4(a) and 4(c) that curiously slowconvergence to a NESS can take place. To shed lighton this result, we first analyze the convergence process and identify a further effective timescale arising from thedriving. This allows us to illustrate the nature of con-vergence within the spin system, given the narrow res-onance regimes of all the different k -components shownin Sec. IV. We then discuss the consequences of this re-lationship between driving and convergence for the pos-sibilities, both theoretical and practical, that NESS maynot be reached at all because the system is driven toostrongly. Finally, the long-time behavior of the NESS inthe absence of driving has both important benchmark-ing properties for theoretical purposes and a key role inthermal control for experimental implementations. Asnoted in Sec. IIC, the formalism we derived there has nolower or upper cutoff in time, and thus can be applied toaddress every aspect of switching on and off a quantumspin NESS. A. Transients at switch-on
In the introduction to NESS in our model (Sec. III), weshowed in Figs. 2(a) and 2(b) how the phonon variablesare “pumped up” on application of the electric field, with n ph ( t ) approaching its steady state, and thus becoming asteady drive for the spin system, after a time of approxi-mately 4 /γ . In Fig. 4(a) we showed how the spin systemreacts to this oscillatory driving, with n x ( t ) approach-ing its steady state after a time of approximately 4 /γ s at most of the driving frequencies in Fig. 4(a). However,it was clear from the n x ( t ) curve at ω /J = 1 .
5, whichtook significantly longer to reach its NESS [Fig. 4(c)],that this reasoning alone does not explain every aspectof the spin response at the onset of driving.For a quantitative analysis of the convergencetimescale, we focus first on in-band driving (2 ω min ≤ ω ≤ ω max ) and consider the process by which n x ( t )is “pumped up” by the driving phonon [Fig. 4(a)]. Wesimplify the analysis by using the fact that, in this rangeof ω , the driven phonon approaches its plateau of con-stant n ph0 more quickly than the spin system, because1 /γ < /γ s . Thus we take the phonon oscillations assinusoidal with a fixed amplitude [Fig. 2(c)], which tomatch the unit slope of Fig. 3 is given by q ( t ) = 2 aγ sin( ω t ) = 2 DJg sin( ω t ) , (45)where we reintroduce the driving parameter, D [Eq. (42)], of the weak-coupling regime. For any selected k -value we define the function f k ( t ) = 4 DJ sin( ω t ) y (cid:48) k [ u k ( t ) + 3 / , (46)which appears as an inhomogeneous term in the lineardifferential equation of Eq. (41c). We combine Eq. (41c)with Eq. (41b) by defining the variable z k ( t ) = v k ( t ) + iw k ( t ), which then obeys the inhomogeneous differentialequation dz k dt = 2 i [ ω k + 2 DJ sin( ω t ) y k ] z k − γ s z k + if k ( t ) . (47)6A suitable primitive of the prefactor of the first term onthe right-hand side is h k ( t ) = 2 (cid:90) [ ω k + 2 DJy k sin( ω t )] dt (48a)= 2 ω k t − DJy k ω cos( ω t ) , (48b)allowing the solution of Eq. (47) to be expressed in theform z k ( t ) = ie ih k ( t ) − γ s t (cid:90) t f k ( t (cid:48) ) e − ih k ( t (cid:48) )+ γ s t (cid:48) dt (cid:48) . (49)This is not yet an explicit expression, because the right-hand side depends on u k ( t ), which remains unknown,but can be related to z k ( t ) by an expression based onEq. (41a),˜ u k ( t ) = u k ( t ) e γ s t (50a)= − iDJy (cid:48) k (cid:90) t sin( ω t (cid:48) )[ z k ( t (cid:48) ) − z ∗ k ( t (cid:48) )] e γ s t (cid:48) dt (cid:48) , (50b)where z ∗ ( t (cid:48) ) denotes the complex conjugate. While thisgeneral expression still does not represent an explicitfunction, it can be used to identify the primary trendsin the response of the spin NESS.We focus on the slowly varying component of n x ( t ),and not on the rapidly oscillating ones. For this it issufficient to consider the slowly varying parts of eachmode occupation, u k ( t ), as may be verified by numer-ical integration of Eqs. (41). Figure 9 indicates thatthe dominant term will be the one at the resonant mo-mentum, k res , which is determined from the driving fre-quency by 2 ω k res = ω . The behavior of k -componentsaway from k res is discussed in App. A and the results aresummarized below. Henceforth we omit the subscript k res . The slowly varying component of the right-handside of Eq. (49) is obtained by averaging over one period, T = 2 π/ω , giving1 T (cid:90) T sin( ω t ) e − ih ( t ) dt = J ( β ) /β, (51)in which β = 4 DJy/ω and J ( β ) is the Bessel functionof the first kind. Replacing sin( ω t (cid:48) ) exp[ − ih ( t (cid:48) )] in theintegrand of Eq. (49) by its average taken from Eq. (51)leads to z ( t ) = i y (cid:48) y ω J ( β ) e ih ( t ) − γ s t F ( t ) (52a)with F ( t ) = (cid:90) t (cid:0) ˜ u ( t (cid:48) ) + e γ s t (cid:48) (cid:1) dt (cid:48) , (52b)which is a real quantity. We stress that the approxima-tions leading to this result are well justified because thedriving oscillations are much faster than the build-up in the triplon expectation values [Fig. 4(a)]. By insertingEq. (52a) into Eq. (50b) we obtain˜ u ( t ) = (2 y (cid:48) ) DJω J ( β ) y Re (cid:104)(cid:90) t sin( ω t (cid:48) ) e ih ( t (cid:48) ) F ( t (cid:48) ) dt (cid:48) (cid:105) (53a)= (cid:18) y (cid:48) ω J ( β ) y (cid:19) (cid:90) t F ( t (cid:48) ) dt (cid:48) , (53b)where again we have used Eq. (51) to obtain the lastexpression. Taking the second derivative yields d ˜ udt = Γ (cid:0) ˜ u ( t ) + e γ s t (cid:1) , (54)in which we have definedΓ = (cid:12)(cid:12)(cid:12)(cid:12) y (cid:48) ω y J (cid:18) DJyω (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) . (55)The differential equation is readily solved with the rel-evant initial conditions, ˜ u (0) = 0 and d ˜ u/dt (0) = 0, togive the final expression for u ( t ) as u ( t ) = 3Γ4 (cid:18) − e − ( γ s − Γ) t γ s − Γ − − e − ( γ s +Γ) t γ s + Γ (cid:19) . (56)This result makes the essential feature clear immedi-ately. At the level of the present analysis, the true con-vergence rate is given by the quantity˜ γ s = γ s − Γ , (57)which can become arbitrarily small when Γ approaches γ s . The qualitative situation is quite intuitive: γ s de-cribes the rate of relaxation of the system back to a statewith zero triplons at zero temperature, which is the caseconsidered here (and discussed in Sec. VC); the phonondriving acts in the opposite direction by creating pairsof triplons, and thus strong driving changes the effec-tive relaxation (damping) timescale. In fact it is clearthat Eq. (57) also specifies a regime where Γ exceeds γ s ,so that triplon creation outweighs the relaxation termand Eq. (56) specifies that the resonant triplon occupa-tion, u ( t ), will undergo an exponential divergence. Thissituation will be the focus of our attention in Sec. VB.Quantitatively, in most circumstances the argument of J will be non-negative and smaller than 1.84, which iswhere the function has its first maximum. In this in-terval, J is a monotonically increasing function of thedriving strength, D , and thus one expects that Γ canindeed be raised to values on the order of γ s .Before computing Γ as a function of ω , we make twofurther general remarks. First, in the qualitative viewof Γ as a driving rate, or excitation rate, that com-petes with the relaxation rate, γ s , one is tempted to in-terpret Γ in terms of Fermi’s Golden Rule. However,a conventional application of the Golden Rule gives arate proportional to the square of the matrix element,7whereas in the present analysis Γ = 2 DJy (cid:48) at small driv-ing ( D → u ( t )in the NESS is given from the long-time limit of Eq. (56)by 3Γ / [2( γ − Γ )] and thus is indeed proportional to Γ ,and hence to D , in accordance with the Golden Rule.However, the timescale of the transient behavior as thesystem approaches the NESS is governed by a differentcoherent mechanism that yields Γ ∝ D .Second, for quantitative purposes it is necessary to con-sider the effect of driving at frequency ω on the modesat k (cid:54) = k res , meaning the action of the driving phononas a “detuned” pump of all other triplon modes. Thealgebra of the detuned case is presented in App. A andwe summarize the results as follows. As a function of adetuning parameter we define as δ = 2 ω − ω , (58)there are two possible regimes. If | δ | < Γ, it is convenientto define the quantity (cid:101)
Γ = (cid:112) Γ − δ , (59)in terms of which u ( t ) = 3Γ γ − (cid:101) Γ ) (cid:20) − e − γ s t (cid:16) cosh( (cid:101) Γ t ) + γ s (cid:101) Γ sinh( (cid:101) Γ t ) (cid:17)(cid:21) . (60)Thus from the behavior of the hyperbolic function, (cid:101) Γadopts the role of Γ in Eq. (57) and the relevant con-vergence rate becomes ˜ γ s = γ s − (cid:101) Γ. By contrast, when | δ | > Γ, so that the detuning of the driving frequencyexceeds the driving threshold, it is convenient to definethe quantity ˜ δ = (cid:112) δ − Γ , (61)in terms of which u ( t ) = 3Γ γ + ˜ δ ) (cid:20) − e − γ s t (cid:16) cos(˜ δt )+ γ s ˜ δ sin(˜ δt ) (cid:17)(cid:21) . (62)Because all the hyperbolic functions become trigonomet-ric, the sole remaining exponential convergence is gov-erned by γ s , leading to the result that the convergence isconventional. We note in this case that slow oscillationsarise at frequency ˜ δ , which may cause the triplon num-ber to overshoot before it converges to its NESS limit(example data not shown).Although one might assume that the resonant case,2 ω k = ω , described by Eq. (56) will provide the high-est threshold value, making (cid:101) Γ max = Γ, the complicateddependence of Γ on k [Eq. (55)] makes it possible that,for a given ω , a slightly detuned mode at k (cid:54) = k res yieldsa higher (cid:101) Γ. In particular, for frequencies close to butoutside the two-triplon band, detuned driving will be ofprimary importance. To capture these possible effects,we compute (cid:101) Γ max by variation of k at each fixed ω , andthe results are shown in Fig. 10. FIG. 10. Dependence of the inverse driving timescale, (cid:101) Γ max ,on the frequency, ω , of the driving phonon, shown for fourvalues of the driving strength, D . Clearly (cid:101) Γ max is finite throughout the in-band regime,although it drops to zero at the band center ( ω = 2 J )due to a matrix-element effect ( y (cid:48) k | k = π/ = 0). Althoughthe dependence of (cid:101) Γ max on ω is both direct and indi-rect, occurring both through proximity to the resonancecondition (2 ω k res = ω ) and through the momentum-dependence of y k and y (cid:48) k , it shows an almost linear risewith frequency towards the two band edges. Because weconsider the linear driving regime of Sec. III, it is also alinear function of D . Importantly, (cid:101) Γ max is also finite out-side the two-triplon band as a consequence of detuneddriving, although for the parameters in Fig. 10 it fallsrapidly (in fact over a frequency window of order DJ )beyond the band edges. For any given D , the function (cid:101) Γ max ( ω ) indicates the values of the triplon damping, γ s ,for which unconventional convergence can occur, and itis no surprise to find that in-band frequencies near thetwo band edges are the most likely candidates [Fig. 4(a)].From Fig. 10, and specifically from the value of (cid:101) Γ max at ω = 2 ω min , one may read that, at the level of our anal-ysis, the value of γ s ensuring conventional or slow con-vergence at all frequencies for driving D = 0 .
01 (Sec. IV)is approximately 0.007. We comment on the minor dis-crepancy with our numerical findings in Fig. 8, whereNESS formation was verified at all ω with γ s = 0 . (cid:101) Γ max , isnot easily read from the external driving parameters,because it depends crucially on the phonon amplitude.Even in the weak-coupling regime, meaning small g as de-fined in Sec. III, we have seen that the oscillations of thedriven phonon are not entirely independent of the spinsystem for in-band driving frequencies. Second, we donot consider the additional complexity of a k -dependent γ s , although the framework developed here could be usedwithout alteration. Third, the effect of non-linear pro-cesses occurring at multiples of ω , is not included in our8 FIG. 11. Creation of the NESS established with a/γ = 0 . g = 0 . J (driving parameter D = 0 .
02) for driving fre-quency ω /J = 1 . γ s = 0 . J ; these arethe parameters of the green line in Fig. 4(a). (a) n ph ( t ). (b) n x ( t ). (c) u k =0 . (d) u k = π . (e) u k res . Also shown is the relax-ation of each variable when the driving is removed after 3000time steps. discussion of Γ and (cid:101) Γ, although it could be incorporatedby considering a very weak effective D .To illustrate the phenomenon of slow convergence atswitch-on, we consider driving field a/γ = 0 . g = 0 . J , which is the situation in Fig. 4(a). For thein-band driving frequency ω = 1 . J and spin damping γ s = 0 . J , we show in Figs. 11(a) and 11(b) the drivenphonon and triplon numbers. The driving strength is thesame as that in Fig. 2(a), and thus n ph ( t ) first rises to-wards the plateau value of 0.04 in a time dictated by 1 /γ ,but is pulled down again to an average value n ph0 (cid:39) . / ˜ γ s . This is a direct reflection of the “inertia” of the spin system as it beginsto absorb some of the input phonon energy, a topic weanalyze in more detail in Sec. VIA. The values of γ s and (cid:101) Γ (Fig. 10) place the system very close to the thresholdspecified by Eq. (57), with the result that the spin NESS[Fig. 11(b)] is reached only after approximately 1200 timesteps [Fig. 4(c)], indicating that 1 / ˜ γ s ≈ /γ s . For a morequantitative understanding of the transient phenomenain this regime, in Figs. 11(c) to 11(e) we show the k = 0-, π -, and k res -components of u k ( t ). It is not a surprise toconfirm that the majority of the slow-convergence behav-ior is indeed concentrated in u k res ( t ) [Fig. 11(e)], whichis both the largest and the most slowly converging com-ponent, apparently requiring 50% longer than n x0 ( t ) toconverge within 2% of its final value. However, it issomewhat surprising to find that the non-resonant u k ( t )components actually rise above their NESS values (on atimescale dictated by γ s ) before falling again as the driv-ing phonon amplitude reaches its final NESS value [onthe timescale dictated by u k res ( t )]. B. Existence of NESS
In Secs. III, IV, and VA, we have used parameters al-lowing the formation of NESS in order to analyze theirresponse to the driving parameters and their internal dy-namical properties. Having established this foundation,we now discuss the crucial issue of whether a NESS canexist at all for strong driving over long driving times.Clearly, unlimited driving would lead to heating of thesystem on a finite timescale, and we defer a discussion ofthis topic until Sec. VI; here we continue to assume thatthe heat sink represented in Fig. 1 maintains a steady,low system temperature despite the injection of energyfrom the laser. The focus of our present discussion is thepossibility that the lattice or spin system could be drivenso strongly that it breaks down rather than converge toa NESS.The integrity of the driven lattice is easy to establish.A straightforward application of the Lindemann crite-rion, whose details we present in App. B, leads to theresult that lattice melting due to phonon driving wouldbecome an issue for average phonon mode occupancieson the order of n ph0 = 3. Thus the driving parameterswe consider here, and the resultant n ph0 values, pose nothreat to the periodic lattice. By contrast, based on thediscussion of Sec. VA, one might expect that Eq. (57)represents a threshold of driving strength ( D ) beyondwhich triplon creation exceeds their relaxation and n x should diverge exponentially, meaning that NESS forma-tion is impossible. Here we discuss two criteria for theloss of NESS. The first is breaching of the condition onthe triplon occupation, n x ( t ) < n maxx = 0 . γ s as defined in Eq. (57).Considering first the maximum triplon occupation,9 FIG. 12. Threshold value, ( a/γ ) t , of the normalized laser elec-tric field strength required to achieve the maximum steady-state triplon occupation of n x = 0 .
2, shown as a function of γ s and ω for fixed g = 0 . J and γ = 0 . ω . We draw atten-tion to the 3 regimes of behavior demarcated by ω (cid:48) = 2 ω min and ω = 2 ω max . Also marked are the frequencies ω (cid:48) = ω min and ω = ω max , where unlike Fig. 8(a) no additional structureis visible in ( a/γ ) t . in Fig. 12 we show the threshold value of the drivingstrength, ( a/γ ) t , required to drive the triplon occupa-tion of the spin NESS above n maxx . Red colors are chosento represent regions of small ( a/γ ) t , because this indi-cates efficient triplon occupation, and these are found atdriving frequencies corresponding to the lower and upperedges of the two-triplon band, intensifying as γ s becomessmaller (Sec. IV). As in Fig. 8(a), it is evident that thesystem does not respond as efficiently for in-band fre-quencies near the band center, and that very strong driv-ing is required when ω lies above the two-triplon band.In contrast to Fig. 8(a), ( a/γ ) t does not reflect the pres-ence of the two-phonon response feature at and above ω = ω min , underlining that the n x0 values arising due tothese processes are indeed small.Nevertheless, one may worry that n maxx = 0 . n ph0 , ex-pected from Fig. 3(a) because of the effects of the spinsystem that it drives. While this result was visible at g = 0 . J for the ω /J = 1 . g regime, as onemay observe from the value of n ph0 in Fig. 6(a). This isone example of a feedback effect between the lattice andspin systems, which we will encounter again in Sec. VIA.Its consequence for the analysis in Sec. VA is that onemay no longer assume a fixed driving strength, D , butbecause of this downwards renormalization one shouldwork with an effective driving, D eff . The deviation be-tween D and D eff becomes larger as D is increased. Wehave also encountered strong feedback effects as γ s be-comes very small and as n ph0 becomes large enough toalter the dimerization, λ , of the spin system. What allof these feedback effects have in common is that theyare significant only when the triplon occupation is large,meaning n x > n maxx , and thus the constraint n maxx = 0 . g , but lies beyond the scope of our currentexposition.Even without a driving criterion, it is nonetheless in-structive to ask how the breakdown of NESS occurs. Weperform only a brief and numerical examination of howthe model of Sec. II behaves when NESS formation isprecluded, for which we set γ s = 0. In Fig. 13 we depictthe time-evolution of the lattice and spin systems for dif-ferent driving frequencies with a fixed driving strength, D = 0 .
01. For in-band driving at ω /J = 1 . γ s = 0 . J , Fig. 13(a) shows the NESS ofFig. 6. However, when γ s = 0, Fig. 13(b) shows how thetriplon number is driven rapidly to a regime well beyond n maxx , which in turn causes the phonon occupation to be-come unstable and creates complex, aperiodic feedbackphenomena.When ω lies above the two-triplon band, one mayread from the detuning discussion, and also directly fromFig. 12, that two possibilities exist. If ω is sufficientlyfar beyond 2 ω max , as shown in Fig. 13(c) for the case ω /J = 3 .
0, a NESS can be formed even with γ s = 0.In this case, the phonon driving cannot cause the directoccupation of triplons and the steadily driven state ofthe spin sector remains only very weakly excited. Anyfeedback from the spin to the lattice sector under thesecircumstances is negligible, and thus the latter is alsounaffected by the value of g . The beating envelope in n x ( t ) in Fig. 13(c) is a consequence of transient signalsin individual components of u k ( t ) that are never dampedwith γ s = 0. The second possibility arises when ω liesabove but very close to 2 ω max , in which case the drivingphonon acts as a detuned pump of the spin response atthe upper band edge and the physics is that of Figs. 13(a)and 13(b).Finally, the situation for driving frequencies below thelower two-triplon band edge is somewhat more compli-0 FIG. 13. Time-dependence of the phonon number, n ph ( t ), andtriplon number, n x ( t ), shown with a/γ = 0 . g = 0 . J ( D = 0 . ω /J = 1 . γ s = 0 . J (theparameters of Fig. 6), the system converges to a NESS ona conventional timescale. (b) When ω /J = 1 . γ s = 0, n x ( t ) increases rapidly beyond n maxx , destabilizing the phononoccupation. (c) When ω /J = 3 .
0, the driving frequency liessufficiently far above the two-triplon band that NESS existeven when γ s = 0. (d) When ω /J = 0 .
75, the drivingfrequency lies well below the two-triplon band but the sec-ond harmonic, 2 ω , lies within it. In this case, when γ s = 0the lattice approximates a NESS, but with this near-constantdriving of the spin system a NESS cannot be formed. cated. Once again there is a regime of potentially diver-gent behavior due to detuned driving when ω lies slightlybelow 2 ω min (Fig. 12). [This phenomenon also allowsone to understand why the lower and upper two-triplonband edges do not create extremely sharp features, oreven discontinuities, as a function of ω in the responseobserved in Figs. 7 and 8.] At frequencies below the de-tuned regime, the generic situation is that illustrated inFig. 13(d) for a frequency ω = 0 . J . The phonondoesindeed approach a NESS, but this essentially steady driv-ing does not create a spin NESS because high-order pro-cesses always exist that pump the undamped spin systemon some potentially very long timescale. In Fig. 13(d) thehigher-order process involves the second harmonic and itis necessary both to follow the spin dynamics to multiplesof 10 time steps and to use long chains (here N = 3000)to verify the situation. In general, driving of the systemby a multi-phonon process can be captured by the sameanalytical arguments applied for in-band frequencies, al-though the effective value of D should be replaced by theamplitude of the relevant higher harmonic. As a result,the qualitative situation at arbitrary below-band frequen-cies is that of Fig. 13(d), but quantitatively the requiredtimescale may extend to millions of steps. We concludethis analysis by stressing again that NESS formation isthe most natural behavior in the model of Fig. 1 at allfrequencies for realistic values of a and γ s , as shown inFig. 13(a), as well as throughout Sec. III. C. Relaxation at switch-off
The process of relaxation of a system with Lindbladdamping is the recovery of thermal equilibrium in theabsence of the drive. In our analysis, the system startedat temperature T = 0 before the drive was switchedon, and thus it relaxes back to this state. The presentanalysis is readily extended to finite T by including (i)a thermal phonon occupation, (ii) a more sophisticatedtreatment of the spin system [82], and (iii) the thermalfactors in the definition of the bath properties that arealready contained in the Lindblad formalism (Sec. IIC).However, this extension would not account for the factthat the driving introduces energy to the system, andhence causes heating; for this we appeal to the heat sinkrepresented in Fig. 1, which corresponds to the coolingapparatus in any condensed-matter experiment. We willdiscuss the issues associated with the system tempera-ture, particularly in the presence of the laser drive, moredeeply in Secs. VI and VII.For the purposes of this subsection, in Fig. 11 we haveswitched off the phonon drive after 3000 time steps. It isclearly visible in all cases that the phonon sector [charac-terized by n ph ( t )] relaxes to its equilibrium, n ph = 0, overa timescale governed by 1 /γ and the spin sector [char-acterized by n x ( t )] over a timescale governed by 1 /γ s .This behavior is independent of the value of ω at whichthe system was being driven and of the amplitude of1 FIG. 14. Schematic representation of energy flow into andout of the NESS of the combined lattice and spin system. the driving (data not shown). In this sense, relaxationcan be considered as similar to the process of “pump-ing up” the NESS with a very low drive, so that thesystem remains far from the driving-induced timescaleobtained in Sec. VA. Thus one may conclude that un-conventional transient processes appear only when thesystem is driven, and indeed driven near its band edges,whereas relaxation dynamics are straightforward.
VI. ENERGY FLOW AND SYSTEM HEATINGA. Energy flow
Particularly valuable for both conceptual and practi-cal purposes is to consider the energy flow through thespin-lattice system. For a true NESS, the rate of energythroughput should be constant from the driving to thefinal stage of dissipation. Figure 14 provides a schematicrepresentation of the situation, which we characterize us-ing seven separate stages of the flow process. The energyflow (energy per unit time) is a power and is definedto be positive in the direction of the arrows in Fig. 14.Clearly the input power is the uptake of laser energy bythe driven phonon, part of which also drives the spinsector through the effect of the spin-phonon coupling.Energy absorption by the lattice, which is also the bath, will determine the temperature of the system; while thisis limited by possible melting of the crystal, the loss ofcoherence in the quantum spin states is a much morestringent criterion. To avoid a monotonic rise in the lat-tice, or bath, temperature, we model the system as beingconnected to a large and efficiently conducting heat sink.To compute the different energy flows in Fig. 14, thefollowing two relations may be read directly from thedifferential equation of Eq. (38c), which describes thetime evolution of the number of energy quanta in thedriven phonon. The energy flow from the laser into thisphonon, normalized to the number of dimers, is given by P L,P ( t ) = − E ( t ) ω p ( t ) , (63)while the energy flowing from it and directly to the bathis given by P P,B ( t ) = γω [ n ph ( t ) − n ( ω )] . (64)The energy flowing out of the driven phonon due to thepresence of the spin-phonon coupling is given by usingthe same equation to obtain P P,SP ( t ) = gω [ U ( t ) + V ( t )] p ( t ) . (65)By considering energy conservation for the phonon weobtain the sum rule P L,P0 = P P,SP0 + P P,B0 (66)for the temporal averages of each power. We remark thatall of the powers in Fig. 14 have a temporal oscillationin the NESS at multiples of the driving frequency, in thesame way as all the other quantites discussed in Secs. IIIand IV. However, these oscillations are not very relevantto the overall energy flow or system temperature andwe focus on their average values, which are the m = 0harmonics of Sec. IV, so we denote them by P X,Y0 .Turning to the spin sector, Eq. (41a) gives the energyflow into the spin system as P SP,S ( t ) = 2 gq ( t ) 1 N (cid:88) k y (cid:48) k ω k w k ( t ) . (67)The same equation also states that the energy flow fromthe spin system into the bath is given by the decay rateof all the triplons, which yields P S,B ( t ) = γ s N (cid:88) k ω k [ u k ( t ) − n ( ω k )] . (68)Once again, energy conservation within the spin systemenforces the sum rule P SP,S0 = P S,B0 (69)on the time-averaged values. However, if one considersEq. (65) as the work done by the phonon on the spinsystem and Eq. (67) as the work received by the spinsystem due to the phonon, it is evident that there is2no mathematical reason for these two quantities to beequal. To obtain the physical sum rule, it is necessary toconsider in detail the spin-phonon coupling term, H sp inEqs. (1) and (21a). In the mean-field approximation, wehave by construction N (cid:104) H sp (cid:105) ( t ) = gq ( t )( U + V )( t ) , (70)and hence the time derivative N ∂ t (cid:104) H sp (cid:105) = g [( ∂ t q )( U + V ) + q∂ t ( U + V )] . (71)Using Eqs. (38a), (41a), and (41b) to evaluate the partialderivatives on the right-hand side yields N ∂ t (cid:104) H sp (cid:105) = gω ( U + V ) p − gγq ( U + V ) (72a)+ gq N (cid:88) k y k (cid:2) gqy (cid:48) k w k − γ s (cid:0) u k − n ( ω k ) (cid:1)(cid:3) (72b) − gq N (cid:88) k y (cid:48) k [2( gqy k + ω k ) w k + γ s v k ] . (72c)By inspection, the first term in Eq. (72a) is P P,SP ( t ) andthe second term in Eq. (72c) is P SP,S ( t ), while the firstterms in Eqs. (72b) and (72c) cancel, as a result of whichthe expression takes the form N ∂ t (cid:104) H sp (cid:105) = P P,SP ( t ) − P SP,S ( t ) (73a) − gq ( t )( γ + γ s )( U + V )( t ) . (73b)The second line suggests rather strongly the definition P SP,B ( t ) = gq ( t )( γ + γ s )( U + V )( t ) (74a)= ( γ + γ s ) (cid:104) H sp (cid:105) ( t ) , (74b)where − P SP,B ( t ) describes a relaxation of (cid:104) H sp (cid:105) ( t ) to-wards zero. This quantity corresponds to a flow of energyfrom the spin-phonon coupling towards the bath and itstemporal average completes the balance P P,SP0 = P SP,S0 + P SP,B0 , (75)which results from the fact that the time-average of thederivative ∂ t (cid:104) H sp (cid:105) must vanish in a NESS. Thus the defi-nition of Eq. (74a) and the additional sum rule of Eq. (75)provide the appropriate linkage to describe energy con-servation in the coupled spin-phonon system.In Fig. 15 we show how the energy flows of Eqs. (63),(64), (65), (67), (68), and (74a) depend on the drivingfrequency. All of the powers we compute obey the steady-state sum rules of Eqs. (66), (69), and (75), which de-scribe every stage of the process. This is clearest whenconsidering the spin system, shown in Fig. 15(b), wherethe forms of P SP,S0 and P S,B0 reflect the exponential in-crease in its sensitivity to driving at frequencies near theedges of the two-triplon band, which we have seen al-ready in Secs. IV and VB. This is particularly true inFig. 12, which can in fact be understood as a graph ofenergy absorption by the spin system (red colors being
FIG. 15. Average energy flows, P X , Y0 , through the combinedlattice and spin system depicted in Fig. 14, normalized to thenumber of dimers and shown as a function of ω for γ = 0 . J , a/γ = 0 . g = 0 . J , and γ s = 0 . J . (a) Power deliveredto the driving phonon ( P L , P0 ) and power dissipated directlyto the bath from this phonon ( P P , B0 ), whose difference is thepower transferred towards the spin system from this phonon( P P , SP0 ). For clarity we show P X , Y0 normalized by ω . (b)Work done by the phonon system due to the spin system( P P , SP0 ) and power delivered to the spin system ( P SP , S0 ), withtheir difference ( P P , SP0 ) represented on the logarithmic scaleas a modulus with a solid (dashed) line for a positive (nega-tive) power. P SP , S0 is identically equal to the power dissipatedby the effect of the bath on the spin system ( P S , B0 ). high values). Although both energy flows bear a closeresemblance to Fig. 8(a), we note from the latter thatthe quantity summed to obtain the net power containsan additional weighting factor of ω k ; among other effects,this acts to make the heights of the peaks at 2 ω min and2 ω max more symmetrical in Fig. 15(b) than in Sec. IV.Turning to the phonon system, in Fig. 15(a) we showthe input energy flow from the laser, P L,P0 , the energyflowing out of the driven phonon due to the spin system, P P,SP0 , and the output flow directly from this phonon3to the bath, P P,B0 . Our first observation is that, in theregime of weak spin-phonon coupling considered here, themajority of the laser energy flows directly to the bath,while the quantity central to our analysis, P P,SP0 , is al-ways relatively small. Next we observe that it peaksaround ω = 2 ω min and 2 ω max , as anticipated fromSec. IV. To illustrate the relative importance of the en-ergy in the spin-phonon coupling term, H sp , we show P P,SP0 once more as the green line in Fig. 15(b) for com-parison with P SP,S0 . Their difference, P SP,B , remains atthe percent level for all driving frequencies within thetwo-triplon band, indicating that H sp does not act tostore significant energy, but in essence transmits it fromthe phonon to the spin system as expected physically.At very high and very low frequencies, | P SP,B | becomesa more significant fraction of the energy in the spin sys-tem, but this energy is in turn a very small fraction of thetotal (laser) energy flowing through the system. We takethese results as evidence that treating the spin-phononterm as a perturbation in the mean-field approach is welljustified, and by extension that the neglect of higher spin-phonon correlations is appropriate for the relevant driv-ing frequencies.We comment in passing that P SP,B can in fact have anegative sign, implying a small energy flow from the bathdue to the spin-phonon coupling term. While this may atfirst appear counterintuitive, we stress that the splittingof the system Hamiltonian into the three parts H p , H s ,and H sp [Eq. (1)] is somewhat arbitrary, and combining H s and H sp would remove this feature. In total, thereis no violation of the fact that energy flows from thecombined spin-phonon system into the bath, and indeedone may compute this net power, P P,B0 + P SP,B0 + P S,B0 ,which by the sum rules at each step of Fig. 14 matches P L,P0 . We do not calculate P B,H0 , assuming simply thatit matches the power flowing into the bath.Because P L , P0 is the average power, or fluence, taken upby the combined spin-lattice system, it is closely relatedto quantities that might be measured in an absorption ex-periment. To make contact with experimental methodsit is necessary to generalize our treatment. In compari-son with a conventional pump-probe procedure, we haveconsidered only the pumping step, because in a NESSthere is no concept of a delay time before probing. Fur-ther, we have considered pumping only at the frequencyof one hypothetical phonon, lying at any value of ω ,which we have varied to probe the behavior of the spinsystem. By contrast, in a real material there is only one,or a small number of, phonon(s) coupled strongly to theprimary magnetic bonds, but it is relatively straightfor-ward to pump the system at all frequencies ω (cid:54) = ω . Thusin Fig. 16 we depart from the conventions used so far inour study and adjust the frequency of the driving laserin order to illustrate the fluence as a function of ω forsystems with one strongly-coupled phonon, which lies ata frequency below, in, or above the two-triplon band.In this type of experiment it is clear that the reso- FIG. 16. P L , P0 shown as a function of ω for γ = 0 . ω , a/γ = 0 . g = 0 . J , and γ s = 0 . J , for systems with onephonon (coupled to the J bond as in Sec. II) at a frequency(a) ω /J = 1 .
0, (b) ω /J = 2 .
0, and (c) ω /J = 3 . nant phonon at ω dominates the absorption. However,the fingerprints of the two-triplon band are visible as, forthe parameters chosen, 0.1%-level effects across the fullfrequency range of the band. This non-resonant absorp-tion is naturally stronger in a material where the relevantphonon lies close to the frequencies 2 ω min and 2 ω max . Werecall that the net fluence at the phonon peak increaseswith ω , and thus that pumping higher-lying phononsmay result in a stronger signal if these are suitably cou-pled to the spin system.4However, for driving a phonon that lies very close toa resonant frequency of the spin system, we draw at-tention to an additional phenomenon. The blue line inFig. 15(a) shows that the absorption peak at ω = ω is actually suppressed when ω lies in the spin band,most strongly so for phonons resonant with 2 ω min and2 ω max . This “self-blocking” effect appears initially to becounterintuitive, as one might expect stronger absorptionwhen more system degrees of freedom are at resonancewith the incoming laser. However, the spin system isnot coupled directly to the light, being excited only bythe driven phonon, and this situation suggests a heuristicimage of the spin system as an extra “inertia” that thedriven phonon must move. While we also used this wordSec. VA, a rather more specific description of the physicscan be read from the prefactor of p ( t ) in Eq. (38c), whereone observes that the spin system acts against E ( t ), mak-ing it more difficult for the phonon to draw energy fromthe laser electric field by oscillating maximally. Oncethis self-blocking effect is taken into account, the differ-ence between the blue and red dashed lines in Fig. 15(a)shows the additional absorption of the incoming fluenceactually taken up by the spin system [shown again inmore familiar form in Fig. 15(b)]. B. Heating
Both conceptually and experimentally, extended con-tinuous driving must inevitably lead to heating, whichwithout remediation would destroy the coherence of thesystem, and later the system itself. Throughout this workwe have assumed that the heat sink represented in Figs. 1and 14 will be able to maintain a constant, low systemtemperature despite the steady drive, and our brief anal-ysis of relaxation to equilibrium in Sec. VC was predi-cated on this assumption. We now turn to a quantita-tive investigation of the reality of the situation in drivencondensed-matter systems.We comment first on the physical meaning of the Lind-blad bath model. Because the energy flowing directlyfrom the driven phonon to the bath, P P,B in Eq. (64), isdirectly proportional to the phonon damping and phononoccupation, γn ph0 /ω of the phonon energy is transferredto the bath in every period. Thus for the parameterswe use, an energy of ω per dimer is dissipated afterapproximately 1500 cycles; we recall that in our modelthe Einstein phonon modes are present on every bondin the system, meaning that the laser driving is a bulkeffect. To introduce some typical numbers for quantummagnetic materials, we consider the inorganic compoundCuGeO , which forms a quasi-1D spin-1/2 system andhas been well characterized in the context of quantummagnetism at equilibrium. In fact CuGeO was studiedin detail [93, 94] due to its spin-Peierls behavior, by whichis meant that it shows a lattice transition from a uniformto an alternating chain that is driven by reducing theenergy in the magnetic sector. This type of transition is a ground-state phenomenon not related to the phonondriving we consider here, and our present analysis wouldin principle be applicable to the distorted state; however,CuGeO also possesses second-neighbor and interchaininteractions, and thus the nature of the transition hasremained the subject of some debate. Although CuGeO is known for its strong spin-phonon coupling, meaningthat its g value lies outside the weak-coupling regime weconsider here (for analyzing driven, out-of-equilibriumphysical properties), we borrow its thermal parametersto analyze energy transfer and heating.CuGeO has leading magnetic exchange constants ofapproximately 10 meV [95], and so we use this value forillustrative energy estimates. Taking the triplon bandwidth into account, a phonon driving the system near2 ω min could also be assumed to have (cid:126) ω = 10 meV ≈ . n ph0 = 0 .
04 and γ = 0 . ω ) that theenergy deposited in the bath by the driving phonon is P P,B0 = 1 . × − Js − per dimer= 5 . × W per mole of spins . (76)As noted in App. B, this value of n ph0 poses no riskof melting the lattice. However, to understand the effectof the energy in Eq. (76) on the lattice temperature, weuse the result [96] that the low-temperature specific heatof CuGeO is given approximately by the standard pure-phonon form, C = βT , with prefactor β ≈ . ). Thus the time required for the driven system toreach a temperature T max in the absence of any coolingapparatus would be t h = β P P,B0 [ T − T ] ≈ . × − [ T − T ] sK . (77)Starting at T init = 0 or 2 K, the time required to reacha temperature T max = 20 K is t h = 2 . × − s. Weassume that this T max is a realistic estimate of the tem-perature where one could no longer argue for quantumcoherence of spin processes taking place in a triplon bandwhose minimum lies at 5 meV. The resulting t h corre-sponds to only 5 cycles of the driving phonon and isclearly too short by a factor of several hundred whencompared with the results of Secs. III to V.Even allowing for considerable latitude with system pa-rameters and materials choices, it is clear that the studyof spin NESS in a quantum magnet is not realistic with-out an efficient heat sink attached to the sample (Figs. 1and 14). To address the effect of the heat sink, it is nec-essary to introduce further materials parameters, specif-ically for sample dimensions and the thermal conductiv-ities removing heat from the sample. As will shortly be-come clear, there are two reasons why an experiment ofthe type we analyze is relevant for a very thin sample,and thus we illustrate the heat flow for a thickness of 20nm. Using that the mass of one mole of spins in CuGeO is 184 g and the density is 5.11 g cm − [97], a sample ofarea A = 1 mm would be 5 . × − moles of CuGeO ,5meaning from Eq. (76) that a laser power P laser = 3 .
26 kW , (78)should be transported through this area to the heat sink.First for the sample itself, the thermal conductivity ofCuGeO at low temperatures is neither constant norisotropic, but an approximate value for the cross-chain( b -axis) direction is κ = 0 . P κ = κA ∆ T / ∆ l = 9 . , (79)where we have set ∆ T = T max − l = 20 nm). Thusthe qualitative conclusion from this crude estimate is thatthe thermal conductivity of the sample can match thepower to be dissipated if a sufficiently thin film can beprepared. In slightly more detail, an energy-flow balancewould dictate that ∆ T should stabilize around 8.5 K.This worked example illustrates that P laser is directlyproportional to ∆ l and P κ inversely proportional, makingthe film thickness a crucial parameter. Nevertheless, thepenetration depth of light into insulating matter is notwell characterized for frequencies where the light is reso-nant with phonon excitations, and further with the spinsector. As a consequence, a thin film is indeed the mostreliable geometry for ensuring that the bulk is uniformlyirradiated by the incident laser beam. Materials thatare difficult to prepare as thin films therefore suffer thetwin disadvantages that their thermal conductivity be-comes a bottleneck in the energy-flow process and thatattenuation of the laser electric field inside the samplebecomes a concern. While a significantly more detailedand materials-specific analysis is required for planning anexperiment, our considerations indicate that it it alwayspossible to study spin NESS in thin-film systems.Certainly an optimized cooling system is a prerequi-site for such studies, even at the nominally weak drivingstrengths ( n ph0 = 0 .
04) we have considered in Secs. IIIto V. The heat sink should be a highly conducting metalable to remove the input power efficiently, and thus nobottleneck should arise due to its thermal contact or thethermal conductivity. However, metals are not knownto have a high heat capacity, and thus we estimate thethermal energy that could be taken up by a metal block.We consider high-quality Al (residual resistivity ratioRRR = 30) and note first that κ Al = 1 W/(K cm) [99],which is well in excess of the value in CuGeO . The spe-cific heat has the form C = γ Al T with γ Al ≈ .
05 J/(kgK ) at low temperatures [99], and hence for a block witharea 1 cm and thickness 5 mm (giving a mass m = 1 .
35 g[99]), the energy absorbed by increasing the temperaturefrom T init = 2 K to T hs is∆ E = γ Al m [ T − T ] = 7 . × − J (80)if the temperature of the heat sink is limited to 5 K. Withan input power of 3.26 kW, the time to overheating of the block is t sink = 2 . × − s , (81)which corresponds to over half a million cycles of the 2.4THz driving phonon. Thus an Al heat sink has plenty ofreserve capacity for the purposes of a NESS experiment.Returning to the beginning of the experimental processdepicted in Figs. 1 and 14, we have assumed that theparameter a is freely variable. The maximum electricfield in a modern THz laser source is approximately E =3 × Vm − [100]. If one assumes that the field acts onan oxygen ion, one obtains a = √ eEq osc = 4 . q osc = 1 . × − m when computed using M O . Forcomparison, the value we have used to ensure n ph = 0 . (cid:126) ω = 10 meV to a = 0 .
16 meV. Thus even for predominantly reflectivesurfaces, values of a suitable for probing the energy rangeof J and ω typical in inorganic quantum magnets arereadily achievable.In summary, experiments of the type we discuss toestablish and to control bulk quantum spin NESS arepossible in real magnetic materials. The sole caveat isthat it should be possible to prepare the system with athickness in the range of tens of nanometres. Even atthe rather modest phonon occupations required to ob-serve nontrivial nonequilibrium spin states, maintainingthe spin system at a low temperature over a long periodof steady driving does pose a significant challenge to thecooling capacity of a conventional cold finger, which nor-mally is designed to control the system temperature withhigh precision using liquid He coolant, rather than func-tioning as an optimized heat sink. We assume that bothof the issues we have identified can be solved for a widerange of quantum magnets. However, in the event of amaterials system that does not allow the driving energyto be removed quickly enough to avoid heating, one solu-tion may lie in altering the experimental geometry awayfrom laser irradiation of the entire sample, as we discussin more detail in Sec. VIIB.
VII. DISCUSSIONA. Approximations: Time, Coupling, and IntensityScales
In constructing our description of the phonon-drivenand dissipative quantum magnetic system we have ap-pealed to a number of approximations. In fact establish-ing the validity of the framework presents an interlinkedproblem involving (i) the treatments we have adoptedfor the laser, for the spin and phonon sectors, and in themaster-equation method, (ii) the fast and slow timescalesof the spin-lattice system, (iii) the coupling constants,6and (iv) the intensities or mode occupations. As exam-ples, the magnetic interactions determine our treatmentof the dimerized chain, a relatively weak spin-lattice cou-pling is intrinsic to our treatment of both sectors, thetimescales of the dynamics in these sectors should al-low the Born-Markov and rotating-wave approximationswithin the quantum master equation, and if mode occu-pations are too high anharmonic or non-linear effects canset in.
1. Triplons as Bosons
A first approximation is that we treat the triplons asnon-interacting bosons, diagonalizing them by a standardBogoliubov transformation. The triplons in a system ofcoupled dimers are in fact hard-core bosons, because atmost one may be present at each site, and finite inter-triplon interactions are well known when the quasipar-ticles are adjacent to each other in real space. How-ever, for relatively low densities (below our threshold of n maxx = 0 .
2) and weak inter-dimer coupling, λ = J (cid:48) /J ≤ .
5, approximating the triplons as interaction-free bosonsis well justified [82], as discussed in Sec. IIA. To studythe regime of larger inter-dimer coupling, the standardBogoliubov transformation can be replaced by a unitarytransformation controlled to high orders in λ [89–91].
2. Laser and mean-field decoupling
As noted in Sec. II, we have described the laser fielddriving the optical phonon as a classical oscillating field.In view of the fluences commonly used in experiment,which make the quantum fluctuations of the laser fieldnegligible relative to its expectation value, this approxi-mation is perfectly justified. The time-dependent mean-field decoupling of the driven phonon and the spin system[Eq. (21a)] is a further approximation, although we havedemonstrated in Sec. VI that it is well justified at all rele-vant driving frequencies. From the definition of Eq. (37), O (10 − ) values n ph on every bond mean that the opticalphonon is macroscopically occupied (the phonon number,being proportional to the system size, is extensive). Thusthe relative size of the quantum fluctations is again neg-ligible, justifying a mean-field treatment of the phononicfield. While we cannot exclude completely that morecomplex physics occurs for particularly large spin-phononcoupling, such as triplon-phonon bound-state formation,this would need to be built first into the ground statesand then into the driven dynamics.
3. Lindblad damping of driven phonon
The damping rate, γ , of the driven phonon should besignificantly smaller than its energy, and we have set itto a value of order 2% of the phonon energy. The way that γ enters, which leads to the description of the drivenphonon as a classical damped harmonic oscillator, is welljustified for the reasons stated in the preceding para-graph. The Lindblad framework treats the relaxationof a degree of freedom by using its established dampingterm. More generally, fundamental theorems about theLindblad formalism [45] state that the dynamics of anopen quantum system can always be captured by decayrates for certain Lindblad operators, which we labelled A l in Eq. (34); at this level, the only open issue is which op-erators in system Hamiltonian are the relevant Lindbladoperators, but this is manifestly obvious for the drivenEinsein phonon considered here.A more physical question concerns the microscopicmechanism of this damping. Clearly, the only bath inan insulator at the energies we consider consists of theoptical and acoustic phonons. Due to anharmonic ef-fects, the driven phonon can decay into two (or more)other phonons. Compared to the single driven phonon,the large number of other phonons in a 3D material con-stitute a large bath, which is not strongly influenced bythe driven phonon, except in cases where the driving ex-ceeds the heat-sink capacities in the sense of Sec. VI andheating effects enter. The fact that the phonon lines ob-served in inelastic scattering studies are usually rathersharp demonstrates that the coupling of a given phononto the bath represented by the other phonons is weak, asa result of which the Born approximation is fully justi-fied.A further required property of the bath is to obey theMarkov approximation, that its correlations should de-cay significantly faster than the decay dynamics of thequantum system (specified by the Hamiltonian). Thisproperty is difficult to verify without a detailed knowl-edge of all the phonons and their anharmonicities, but anestimate is possible. In general the spectrum of phononscovers the energy range from zero to the Debye energy, (cid:126) ω D , and hence 1 /ω D sets the timescale for the decay ofbath correlations. The Debye energy is typically 50-100meV (12-24 THz). For driving frequencies in the 1-10THz regime and a decay rate, γ , which is 1/50 of these,it is clear that the correlation time scale, 1 /ω D , is indeedshorter than the time scale 1 /γ of the phonon damping(with the possible exception of very soft materials). Fi-nally, the rotating-wave approximation is the statementthat one may neglect fast oscillations to focus only onthe slow variables (as we did in Sec. VA) [45], and againthis is clearly justified because the oscillatory terms fora phonon driven at ω are fast on the timescale of thedamping.
4. Lindblad damping of triplons
The previous arguments can be repeated to justify theuse of t k as the Lindblad operators in the spin sector.We observed in Sec. IIC that these operators break spinconservation, so that our treatment is relevant for sys-7tems with finite spin-orbit coupling. If this coupling islow, spin conservation requires that one consider termsof the type C kq = t † k t q , which we discuss in the nextsubsection (Sec. VIIB). However, weak spin-orbit cou-pling also implies that the damping of spin excitationsdue to a phononic bath is weak, and thus it is not un-reasonable to treat any magnetic excitations, which hereare the triplons, as weakly damped oscillators. Becausethe strongest effects of driving the magnetic system occurwhen the driving frequency, ω , matches the magnetic en-ergies, 2 ω k , estimates for the validity of the Born-Markovand rotating-wave approximations [45] are the same asfor the driven phonon. We defer comments on momen-tum conservation in C kq to Sec. VIIB.In summary, while the validity issue is a complex one,all of the approximations we have made are appropriate,and in fact for a typical condensed-matter system thereis a reasonable amount of parameter space (Sec. VI). Abroader discussion of materials and experiment may befound in Sec. VIIC below. From a theory standpoint,our current approach is by design the simplest avail-able, whose explicit intent is to establish the basic phe-nomenology, and a more detailed discussion of any givenissue may require more sophisticated methodology. Oneexample of this would be the use of flow-equation meth-ods [92] to extend the regimes of validity, in the hierarchyof timescale approximations, of the equations of motion.We close this part of the discussion by recalling thatthe intrinsic properties of the phonon-driven spin systemlead to a number of phenomena occurring over a rangeof different frequencies and times. By frequency, the keyregimes of driving are (i) in the spin band, where the re-sponse is resonant, (ii) below the spin band, where it iscontrolled by multiphonon processes, and (iii) above thespin band, which is the Floquet regime, featuring weakenergy absorption and coherently superposed phase- andfrequency-shifted states. By time, transient phenom-ena at switch-on occur (mostly) on the scale of the in-verse damping (Sec. VA), drive-induced heating occurson a strongly ω -dependent timescale, and relaxationphenomena at switch-off follow the Lindblad form to re-store the starting state (Sec. VC). B. Bath Models and System Heating
The equations of motion whose solutions we have stud-ied in Secs. III to V are intrinsic to one type of bath. Asstated in Sec. IIC, the physical content of the Lindbladformalism is that the spin operators are damped by bathoperators that also appear in the spin Hamiltonian. How-ever, the nature of these terms must reflect the physics ofthe entire system, by which is meant the manner in whichenergy can be dissipated by spin and phononic processes.Here we comment briefly, and with one specific example,on how our analysis would be extended in the case ofmore complex bath terms.It is clear in Sec. IIC that our use of t k as the sole type of spin-bath operator delivers the most straight-forward equations of motion, and we have exploited thecomplete independence of all k states to explore a widerange of phenomena. However, it is also clear that damp-ing processes involving a single t k operator are spin-non-conserving, allowing a triplon to decay into a phonon. Ina truly 1D system, meaning a spin chain in a 1D lattice,momentum conservation would restrict the phase spaceavailable for such processes, raising problems with the ap-plicability of the Lindblad formalism (which requires thatthe bath contain a continuum of energy states [45]) thatwould at minimum mandate a significant k -dependence ofthe decay rate, γ s ( k ). However, as explained in Sec. IIC,such concerns are not relevant to the present analysisbecause the 1D triplons are embedded in a real lattice,meaning with 3D phonons. In this case, for each mo-mentum, (cid:126) k , along the spin chains, there is a contin-uum of perpendicular momenta, (cid:126) (cid:126)k ⊥ , and hence the an-nihilation of the triplon can occur for a wide range of (cid:126)k ⊥ values, with the bath phonons that are created cov-ering a broad energy range. This energetic continuumwill depend on k , but only weakly, which both justifiesusing the Lindblad formalism and indicates a constant,momentum-independent damping rate, γ s .Nevertheless, the conventional spin-phonon coupling inany 3 d transition-metal compound, and hence the spin-damping effect of its phononic modes, takes the form of H sp in Eq. (17), and the spin-isotropic nature of thisinteraction means that phonon modes cannot alter thespin state (the number of excited triplons) directly. Themost straightforward spin-conserving bath operators ap-propriate to this situation, based on the operators t k in the spin-system Hamiltonian, would be of the form C kq = t † k t q , and for the reason given above need not bemomentum-conserving. Bath operators with the form of C kq manifestly act to mix wave-vector states of the sys-tem and thus lead to a significantly more involved set ofequations of motion, with in general N coupled equa-tions rather than only N . We defer a detailed analysis ofthis case to a follow-up study.As already noted, the type of bath studied inthe present work provides a meaningful description ofsystems with spin-dependent phonon scattering pro-cesses. These can arise in systems with apprecia-ble spin-orbit coupling, meaning 4 d and especially 5 d magnetic ions, where the resulting anisotropic interac-tions include may Dzyaloshinskii-Moriya (DM), exchangeanisotropies (XXZ and XYZ), or even bond-selective in-teractions. However, only in rather exceptional circum-stances would these dissipative channels be stronger thanspin-conserving damping terms, and thus the considera-tion of more advanced bath operators is required to dis-cuss real experiments. In addition to the question ofspin conservation, it is also necessary to address the is-sue of spin-system dimensionality, which as in the presentstudy may be lower than the phonon-system dimension-ality (which is 3D), and hence to establish the level atwhich to enforce conservation of momentum.8The nature of the bath reflects directly on the heatingof the system, which was discussed for the simple, spin-non-conserving case in Sec. VI. In a more complex bath,one may expect the redistribution of energy through themodes of the spin system to be more efficient, althoughthis is in general a small contribution to the (phonon-dominated) system temperature. Of more practical rele-vance to the issue of quantum coherence is the fact thata momentum-mixing bath operator would also impactthe coherence of individual k -components of the spin sys-tem. From this standpoint one may consider “reservoirengineering” [84], meaning influencing the form of C kq (for example by promoting forward-, backward-, or skew-scattering in the bath), as an alternative to controllingthe system temperature only through the balance be-tween the laser driving strength and the cooling appa-ratus (Sec. VI).With a view to maintaining quantum coherence in thespin sector, we turn to some more general considerationsfor controlling the temperature of the system. We re-mark that in some classes of system it is possible to de-couple the degrees of freedom in such a way as to obtaineffective electronic or spin temperatures different fromthe lattice temperature. However, this is not an optionin a system of spins localized on the sites of a lattice,where the temperature controlling the response of thespin system is that of the lattice. Similarly, our modelis also far from the paradigm of a spatially separate sys-tem and bath, where different effective temperatures forthe two components are related by controllable couplingconstants. Within the confines of the situation we con-sider, we mention two approaches to temperature control,namely system geometry and the laser driving protocol.In the present work we have considered only bulk driv-ing by the electric field of the laser, meaning that the Ein-stein phonon of every bond is stimulated. In Sec. VIB weshowed that this “bulk” system should in fact be ratherthin (tens of nanometres). However, it is certainly pos-sible that a device of µ or mm length is illuminated onlyat one end, causing the phonon and spin excitations topropagate through the equilibrium material over a dis-tance far larger than the nonequilibrium irradiated por-tion, and possibly larger than the penetration depth ofthe light. Such a situation would require a model for spa-tial gradients of heat, magnetization, and temperature,which would certainly be of direct interest for switchingand transport in spintronic devices. To date some experi-ments already present this type of situation [68], and cer-tain theoretical discussions have also invoked the frame-work of driving only at the ends of the system [47, 48, 85].Finally, another means of controlling the system tem-perature lies in driving by repeated short pulses. In itssimplest form, this allows the system to relax back toits cold state by the action of the heat sink (Sec. VC),although the required pulse separation would be a veryslow timescale (Sec. VIB). At a more sophisticated level,pulsed driving processes allow new degrees of control overthe system, including the imposition of dynamics on new timescales quite separate from the driving frequency, asin the case of Floquet engineering of the electronic bandstructure [34]. While certain types of driving protocolhave already been proposed for controlling small num-bers of quantum spins [4] and ensembles of effectively S = 1 / C. Experiment
As stated in Sec. I, the last decade has seen an enor-mous expansion in the technological capabilities of lasersources, both in ultrafast timescales and in high intensi-ties, with applications both for pump-probe experimentsand for steady driving. Where in Sec. I the focus of ourremarks was the new physics made possible by these newsources, it is also worth commenting on the new tech-nologies that have led to such growth in the applicationof lasers to condensed matter. For decades this was lim-ited by the “Terahertz gap,” the problem that light atthe energies of most interest to the intrinsic processes incondensed matter was neither easily generated nor eas-ily guided or focused, but was easily absorbed and scat-tered. Starting with the initial compilation of methodsmaking it possible to engineer transient states of con-densed matter [86], further technological solutions havebeen developed and applied to frontier science challenges[87]. The best review of terahertz enabling technologies,both for generation and for beam control, may be foundin Ref. [88]. On the generation side, one has not onlynew free-electron laser sources but also a range of new“table-top” techniques, including plasma-based sourcesand high-harmonic generation, many made possible byexploiting new materials. On the control side, beamguiding, transport, focusing, and diagnostics have alsobenefited strongly from the optical properties of certainmaterials. Together this progress has led to a qualitativeexpansion in the type of physics that can be probed, orindeed created, in condensed matter, and the aim of thepresent study is to extend these capabilities to quantummagnetic materials.For the purposes of this preliminary analysis, we havefocused on well-dimerized (and thus robustly gapped)quantum spin chains, meaning that the system we con-sider does not, either at equilibrium or in its driven state,approach a phase transition to a magnetically ordered,to a gapless quantum disordered, or to any other differ-ent magnetic state. In a 1D system with only Heisen-berg spin interactions, the primary requirement is sim-ply that the spin gap (the one-triplon band minimum, ω min = ω k =0 ) does not close, including on laser driving ofa selected phonon. An excellent example of an inorganiccompound realizing quasi-1D alternating S = 1 / ) [101], which is thought to have noanisotropy, negligible second-neighbor interactions, anda gap of approximately 0.38 meV (compared to a band9width of 0.12 meV, yielding λ = 0 . ) are lower by a factor of20 than the test-case numbers presented in Sec. VI, pre-senting a different balance of slower heating rates, slowerconvergence to NESS, and altered damping ratios.A recently discovered class of alternating spin-chainmaterials includes AgVOAsO [102] and NaVOAsO [103], which have magnetic coupling constants in the5 meV range. Although here λ is at the upper valid-ity limit of our present simple treatment of the spinchain (Sec. IIA), as noted earlier, more sophisticated ap-proaches are available for this purpose and no part ofthe equations of motion (Sec. IIC) is invalidated. An-other class of candidate systems is the set of metal-organic TTF compounds [104, 105], and even purely or-ganic TCNQ compounds [106], in which the spin-Peierlstransition has been observed and the distorted (low-temperature) state is an alternating spin chain. A furthercategory of interest in quantum magnetism has been al-ternating antiferromagnetic-ferromagnetic (AF-FM) S =1 / Cu SbO ( λ = − .
79) [107] and (CH ) NH CuCl ( λ = − .
92) [108], our analysis remains fully applicableregardless of the signs of the interactions.In addition to CuGeO , which we introduced in Sec. VIto consider its thermal properties, (VO) P O [109, 110]constitutes a further system that in fact realizes al-ternating S = 1 / g values that made both of these compounds attractivefor equilibrium experiments in quantum magnetism doplace them outside the weak-coupling regime we ana-lyze here. In a later study we will extend our considera-tions to the regime of strong spin-phonon coupling, whichis also described by the equations of motion derived inSec. II. Here one may anticipate nonlinear driving effects(Figs. 3 and 5), which could allow experiments at lowerlaser intensities, stronger mixing of frequency harmonics(Figs. 6 and 7), stronger below-band multi-phonon pro-cesses (Fig. 12), more delicate driving-induced anomalousconvergence (Sec. VA), and stronger transfer of spectralweight between different frequencies (Fig. 16).Returning to the spin sector alone, our present formal-ism is readily extended to alternating chains with differ-ent gap-to-bandwidth ratios, although an accurate treat-ment of systems with small gaps would require a morenumerical approach for the systematic summation of per-turbative terms to high orders, as in the method of con-tinuous unitary transformations (CUTs) [89–91]. Other1D gapped spin systems to which our considerations areimmediately applicable include even-leg S = 1 / S = 1) chains. A parallel classof systems that could be treated by very similar con- siderations would be that of magnetically ordered quan-tum spin systems, which includes many 2D and 3D ma-terials of all spin quantum numbers; this would involvethe straightforward adoption of a (constrained) spin-waveframework to describe the spin sector. Ordered mag-netic systems also provide the simplest cases in whichto analyze the effects of anisotropic magnetic interac-tions, such as single-ion, DM, XXZ, and other terms,which have recently attracted intensive interest with aview to creating topological magnons [54–57], vortices[59], skyrmions [60], and other means of encoding pro-tected quantum information [2]. More complex spin sec-tors include anisotropic systems such as Ising and XYmodels without magnetic order, gapless spin chains, andgapped or gapless non-ordered states in higher dimen-sions, meaning in the former category Z quantum spinliquids and in the latter algebraic spin liquids and quan-tum critical systems. Here the challenge is not only tofind a suitable framework in which to describe the com-plex correlated spin sector, especially if this is changedby using laser driving to push it across a magnetic quan-tum phase transition, but also to deal with the situationwhere the excitations of the spin system extend to arbi-trarily low energies, thus interacting strongly with eventhe acoustic phonons.In addition to an adequate treatment of the spin sec-tor, the quantitative analysis of real materials will requireaccurate lattice dynamics calculations to give the phononmodes and frequencies, and the corresponding oscillatorstrengths. The normal modes can be used to estimatespin-phonon coupling strengths and the frequencies tochoose the laser excitation parameters. Typically, thephonon spectrum in inorganic materials extends up to theDebye energy, which rarely exceeds 100 meV (24 THz),while a lower limit for optical phonons is perhaps 5-10meV. Spin energy scales can extend up to a one-magnonband maximum of 300 meV (70 THz) in cuprates, andhave no lower limit. There is no established relationshipbetween the two, as materials with predominantly high-energy phonon modes can have a very low-energy spinsector, and those with high spin energies need have nospecial phonon properties. However, metal-organic sys-tems do tend to have a softer phonon sector, due to thenature of the interactions between weakly polar organicgroups, and very low magnetic energies due to the longpaths between magnetic ions. Thus it is safe to say thata very wide range of frequency scales and phononic mod-ulation possibilities is available for planning experimentsof the type we discuss.Finally, it is also necessary to consider how to measureall of the physical quantities characterizing the spin sys-tem using a terahertz laser. In Sec. VIA we presentedthe example of absorption of the incident laser fluence(Fig. 16), while further quantities that are also functionsof frequency and temperature include reflection, polariza-tion rotation (due to birefringence changes or the Faradayeffect), and the two-magnon response. It is not generallypossible to probe the wave-vector response of the spin0sector, except in special cases possessing a strong cou-pling to one well-characterized “probe phonon” mode.In most such measurements, the signal arising due to thespin system will be weak in comparison with the manyother contributions to the total response of a sample, andhere we point to the strong frequency-selectivity allowedby the phonon-coupled model, and visible in Figs. 15 and16, as the primary means of ensuring that the spin signalis detectable. VIII. SUMMARY
We have investigated the nonequilibrium steady quan-tum mechanical states of a lattice spin system under thecontinuous, coherent laser excitation of lattice phononsat a single frequency and in the presence of a realisticdissipation. In real materials, this dissipation is domi-nated by the many phonons of the lattice system, and itsinclusion at the operator level, which we effect within theLindblad formalism, goes well beyond much of the workon driven many-body systems presently in the literature.We have focused for illustrative simplicity on a dimerizedspin chain driven by a bulk Einstein phonon.By establishing and solving the quantum master equa-tions governing the time-evolution of this model system,we have demonstrated the establishment of quantum spinNESS and their dependence on the system parameters.We have performed a detailed analysis of the internal dy-namics of these states, and of the accompanying behaviorof the driven lattice system, in the presence of Lindbladdissipation. We find a dramatic but readily comprehen-sible sensitivity of the response to the driving frequencyand analyze the factors influencing the response of thequantum spin NESS at different harmonics of this fre-quency. Our equations of motion are valid at all timesand we use them to study both the transient behavior ofthe system when driving is switched on, where we finda complex phenomenology, and the relaxation to equilib-rium when the driving is removed. We have computedthe energy flow through the composite spin-lattice sys-tem, from its arrival as the driving laser light to its dissi-pative loss, which allows us to gauge the self-consistencyof our analysis by applying sum rules. The Lindblad for-malism allows direct access to the energy flow into thebath, from which we have estimated heating timescalesand hence the practical requirements, in the form of driv-ing limits, sample geometry, and cooling capacity, of aNESS experiment.The framework we establish makes it possible to per-form quantitative investigations of many different typesof spin system, including those with magnetic order, withsmall or vanishing gaps, with topological properties, orwith nontrivial quantum entanglement. It also enablesthe analysis of more complex types of bath, most notablyones describing spin-conserving dissipative processes, andhence the modelling of real materials with laser-drivenphonons. With appropriate treatment of spatial gradi- ents (in driving, magnetization, and temperature), onemay also model real device geometries, leading to spin-tronic applications where the challenge is to preservequantum coherence over the timescale required for read-ing and processing the quantum information encoded inthe spin sector.
ACKNOWLEDGMENTS
We thank D. Bossini, F. Giorgianni, S. Haas, C. Lange,Z. Lenarˇciˇc, A. Rosch, Ch. R¨uegg, and B. Wehinger forhelpful discussions. Research at TU Dortmund was sup-ported by the German Research Society (DFG) throughGrant No. UH 90-13/1 and together with the RussianFoundation of Basic Research through project TRR 160.Research at LANL was supported by the U.S. DOENNSA under Contract No. 89233218CNA000001 throughthe LANL LDRD Program.
Appendix A: Detuned Triplon Pair Creation
In Sec. VA we analyzed the gradual change of the diag-onal triplon components, u k ( t ), in the resonant case, bywhich is meant for ω = 2 ω k res . Henceforth we omit thesubscript k . If this resonance condition is not met, westate that there is a finite detuning, δ = 2 ω − ω . In thisappendix we present the differences between the deriva-tions for the detuned and resonant cases (Sec. VA). Inthe detuned case, Eq. (52a) still holds, but the definitionof Eq. (52b) changes to F ( t ) = (cid:90) t (cid:0) ˜ u ( t (cid:48) ) + e γ s t (cid:48) (cid:1) e − iδt (cid:48) dt (cid:48) , (A1)which implies that F ( t ) is no longer real, but complex.Equation (53a) still holds, but Eq. (53b) is modified to˜ u ( t ) = e γ s t u ( t ) = Γ Re (cid:34) (cid:90) t e iδt (cid:48) F ( t (cid:48) ) dt (cid:48) (cid:35) , (A2)with Γ as defined in Eq. (55). We stress that the validityof replacing the rapidly oscillating terms by their average,as performed in Eq. (51), is justified if | δ | (cid:28) ω . If thedetuning becomes too large, the deviations from the fullresult may become large, but comparing our analyticalapproximation to the results of a numerical integrationrevealed very good agreement over a broad range of pa-rameter space.Equation (A2) does not lead to a closed differentialequation by double differentiation because of the re-striction caused by taking the real part. It is neces-sary instead to take three derivatives of ˜ u ( t ) and define x ( t ) = d ˜ u/dt ( t ), for which we obtain d xdt ( t ) = (Γ − δ ) x ( t ) + Γ γ s e γ s t . (A3)1Clearly, the nature of the solutions to this differentialequation depends crucially on the sign of the prefactor,Γ − δ . As in the resonant case (Sec. VA), if it is positivethen exponentially increasing and decreasing functionsappear, whereas if it is negative then oscillating trigono-metric functions appear. Solving Eq. (A3) for the initialconditions x (0) = 0 and dx/dt (0) = 3Γ / u ( t ), from whichthe expressions given for u ( t ) in Eqs. (60) and (62) follow. Appendix B: Lindemann Criterion
Breakdown of the lattice is governed by the Linde-mann criterion [111], which dates to Lindemann’s intro-duction of the concept that a solid will begin to meltwhen (cid:104) q (cid:105) ≈ ρ a , meaning when the fluctuations, q ,of its atoms around their equilibrium positions exceed afraction, ρ , of the interatomic distance, a . While Lin-demann used the concept to relate (cid:104) q (cid:105) to the meltingtemperature, T m , for our purposes it is sufficient to re-late (cid:104) q (cid:105) to n ph .It has been established for a broad range of condensed-matter systems that Lindemann’s proposed relationholds, with the value of ρ being in the range 0.1-0.15 [112–114]. To relate this result with n ph , for a local phononicoscillation the potential energy is half of the total energy, M ω (cid:104) q (cid:105) = (cid:126) ω ( n ph + ) , (B1)which implies (cid:104) q (cid:105) q = n ph + , (B2) where q osc is the oscillator length, (cid:112) (cid:126) / ( M ω ). The Born-Oppenheimer approximation [115] implies that the char-acteristic lengthscale of the atomic motion, q osc , is a frac-tion ( m e /M ) / of the electronic lengthscale, which weidentify crudely with a ; here m e is the mass of an elec-tron and M the mass of the oscillating atom. Thus oneobtains from Eq. (B2) that (cid:104) q (cid:105) a (cid:114) Mm e = n ph + , (B3)and hence n ph + ≈ ρ (cid:114) Mm e . (B4)Using the mass of the oxygen atom as a generic value for M , and ρ ≈ .
02 for the Lindemann ratio, we concludethat the phonon number should not exceed n ph ≈ n ph is the to-tal number per atom of phonons polarized in one direc-tion and the conventional Lindemann criterion appliesbecause we assume the chains to be embedded in a 3Dcrystal (no low-dimensional instability need be consid-ered). From the estimates made above, it is clear thatthe Lindemann threshold is in no way threatened by thedriven phonon occupations we consider, and thus thatthe integrity of the periodic lattice is not an issue. [1] V. Baltz, A. Manchon, M. Tsoi, T. Moriyama, T.Ono, and Y. Tserkovnyak, Antiferromagnetic spintron-ics, Rev. Mod. Phys. , 015005 (2018).[2] A. V. Chumak, V. I. Vasyuchka, A. A. Serga, and B.Hillebrands, Magnon spintronics, Nature Phys. , 453(2015).[3] D.-Q. Jiang, M. Qian, and M.-P. Qian, MathematicalTheory of Nonequilibrium Steady States (Springer, Hei-delberg, 2004).[4] G. S. Uhrig, Quantum coherence from commensuratedriving with laser pulses and decay, SciPost Phys. ,040 (2020).[5] M. S. Rudner, N. H. Lindner, E. Berg, and M. Levin,Anomalous Edge States and the Bulk-Edge Correspon-dence for Periodically Driven Two-Dimensional Sys-tems, Phys. Rev. X , 031005 (2013).[6] L. DAlessio and A. Polkovnikov, Many-body energylocalization transition in periodically driven systems,Ann. Phys. , 19 (2013).[7] R. Nandkishore and D. A. Huse, Many-Body Localiza-tion and Thermalization in Quantum Statistical Me-chanics, Annu. Rev. Condens. Matter. Phys. , 15(2015). [8] K.-O. Chong, J.-R. Kim, J. Kim, S. Yoon, S. Kang, andK. An, Observation of a non-equilibrium steady state ofcold atoms in a moving optical lattice, Commun. Phys. , 25 (2018).[9] R. Labouvie, B. Santra, S. Heun, and H. Ott, Bistabil-ity in a Driven-Dissipative Superfluid, Phys. Rev. Lett. , 235302 (2016).[10] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Luschen,M. H. Fisher, R. Vosk, E. Altman, U. Schneider, and I.Bloch, Observation of many-body localization of inter-acting fermions in a quasirandom optical lattice, Science , 1547 (2015).[11] A. Eckardt, Colloquium: Atomic quantum gases in pe-riodically driven optical lattices, Rev. Mod. Phys. ,011004 (2017).[12] I. Bloch, J. Dalibard, and W. Zwerger, Many-bodyphysics with ultracold gases, Rev. Mod. Phys. , 885(2008).[13] I. Gierz, Probing carrier dynamics in photo-excitedgraphene with time-resolved ARPES, J. Electron Spec-troscopy and Related Phenomena , 53 (2017).[14] A. Cavalleri, Photo-induced superconductivity, Con-temp. Phys. , 31 (2018). [15] M. Buzzi, M. F¨orst, R. Mankowsky, and A. Cavalleri,Probing dynamics in quantum materials with femtosec-ond X-rays, Nature Reviews Materials , 299 (2018).[16] A. Zong, X. Shen, A. Kogar, L. Ye, C. Marks, D.Chowdhury, T. Rohwer, B. Freelon, S. Weathersby, R.Li, J. Yang, J. Checkelsky, X. Wang, and N. Gedik, Ul-trafast manipulation of mirror domain walls in a chargedensity wave, Science Adv. , eaau5501 (2018).[17] S. Iwai, M. Ono, A. Maeda, H. Matsuzaki, H. Kishida,H. Okamoto, and Y. Tokura, Ultrafast Optical Switch-ing to a Metallic State by Photoinduced Mott Transitionin a Halogen-Bridged Nickel-Chain Compound, Phys.Rev. Lett. , 057401 (2003).[18] A. Zong, A. Kogar, Y.-Q. Bie, T. Rohwer, C. Lee, E.Baldini, E. Erge¸cen, M. B. Yilmaz, B. Freelon, E. J. Sie,H. Zhou, J. Straquadine, P. Walmsley, P. E. Dolgirev,A. V. Rozhkov, I. R. Fisher, P. Jarillo-Herrero, B. V.Fine, and N. Gedik, Evidence for topological defects ina photoinduced phase transition, Nature Phys. , 27(2019).[19] T. Kitagawa, T. Oka, A. Brataas, L. Fu, and E. Demler,Transport properties of nonequilibrium systems underthe application of light: Photoinduced quantum Hall in-sulators without Landau levels, Phys. Rev. B , 235108(2011).[20] N. H. Lindner, G. Refael, and V. Galitski, Floquet topo-logical insulator in semiconductor quantum wells, Na-ture Phys. , 490 (2011).[21] P. Rodriguez-Lopez, J. J. Betouras, and S. E. Savel’ev,Dirac fermion time-Floquet crystal: Manipulating Diracpoints, Phys. Rev. B , 155132 (2014).[22] T. Iadecola, D. Campbell, C. Chamon, C.-Y. Hou,R. Jackiw, S.-Y. Pi, and S. V. Kusminskiy, Materi-als Design from Nonequilibrium Steady States: DrivenGraphene as a Tunable Semiconductor with TopologicalProperties, Phys. Rev. Lett. , 176603 (2013).[23] M. A. Sentef, M. Claassen, A. F. Kemper, B. Moritz,T. Oka, J. K. Freericks, and T. P. Devereaux, Theory ofpump-probe photoemission in graphene: Ultrafast tun-ing of Floquet bands and local pseudospin textures, Na-ture Commun. , 7047 (2015).[24] M. F¨orst, C. Manzoni, S. Kaiser, Y. Tomioka, Y.Tokura, R. Merlin, and A. Cavalleri, Nonlinear phonon-ics as an ultrafast route to lattice control, Nature Phys.
854 (2011).[25] A. Subedi, A. Cavalleri, and A. Georges, Theory ofnonlinear phononics for coherent light control of solids,Phys. Rev. B , 220301 (2014).[26] A. von Hoegen, R. Mankowsky, M. Fechner, M. F¨orst,and A. Cavalleri, Probing the Interatomic Potential ofSolids with Strong-Field Nonlinear Phononics, Nature , 79 (2018).[27] A. L¨auchli and C. Kollath, Spreading of correlationsand entanglement after a quench in the one-dimensionalBose-Hubbard model, J. Stat. Mech. P05018 (2008).[28] S. A. Hamerla and G. S. Uhrig, One-dimensionalfermionic systems after interaction quenches and theirdescription by bosonic field theories, New J. Phys. ,073012 (2013).[29] S. Paeckel, B. Fauseweh, A. Osterkorn, T. K¨ohler, D.Manske, and S. R. Manmana, Detecting superconduc-tivity out of equilibrium, Phys. Rev. B , 180507(R)(2020).[30] L. Schwarz, B. Fauseweh, N. Tsuji, N. Cheng, N. Bit- tner, H. Krull, M. Berciu, G. S. Uhrig, A. P. Schnyder,S. Kaiser, and D. Manske, Classification and characteri-zation of nonequilibrium Higgs modes in unconventionalsuperconductors, Nature Commun. , 287 (2020).[31] B. Fauseweh and J.-X. Zhu, Laser pulse driven controlof charge and spin order in the two-dimensional Kondolattice, to appear in Phys. Rev. B (arXiv:2002.03023).[32] M. Rigol, Quantum quenches and thermalization in one-dimensional fermionic systems, Phys. Rev. A , 053607(2009).[33] J. Berges, Sz. Bors´anyi, and C. Wetterich, Prethermal-ization, Phys. Rev. Lett. , 142002 (2004).[34] T. Oka and S. Kitamura, Floquet Engineering of Quan-tum Materials, Annu. Rev. Condens. Matter Phys. ,387 (2019).[35] E. Arrigoni, M. Knap, and W. von der Linden, Nonequi-librium Dynamical Mean-Field Theory: An AuxiliaryQuantum Master Equation Approach, Phys. Rev. Lett. , 086403 (2013).[36] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, andP. Werner, Nonequilibrium dynamical mean-field theoryand its applications, Rev. Mod. Phys. , 779 (2014).[37] T. Qin and W. Hofstetter, Nonequilibrium steady statesand resonant tunneling in time-periodically driven sys-tems with interactions, Phys. Rev. B , 125115 (2018).[38] M. Eckstein and P. Werner, Ultrafast Separation ofPhotodoped Carriers in Mott Antiferromagnets, Phys.Rev. Lett. , 076405 (2014).[39] Y. Murakami and P. Werner, Nonequilibrium steadystates of electric field-driven Mott insulators, Phys. Rev.B , 075102 (2018).[40] A. Herrmann, Y. Murakami, M. Eckstein, and P.Werner, Floquet prethermalization in the resonantlydriven Hubbard model, Europhys. Lett. , 57001(2018).[41] S. A. Sato, P. Tang, M. A. Sentef, U. De Giovannini, H.H¨ubener, and A. Rubio, Light-induced anomalous Halleffect in massless Dirac fermion systems and topologi-cal insulators with dissipation, New J. Phys. , 15767(2017).[43] Z. Lenarˇciˇc, E. Altman, and A. Rosch, Activating many-body localization in solids by driving with light, Phys.Rev. Lett. , 267603 (2018).[44] G. Lindblad, On the generators of quantum dynamicalsemigroups, Comm. Math. Phys. , 119 (1976).[45] H.-P. Breuer and F. Petruccione, The Theory of OpenQuantum Systems , 2nd Ed. (Oxford University Press,Oxford, 2007).[46] U. Weiss,
Quantum Dissipative Systems , 2nd Ed.(World Scientific, Singapore, 2012).[47] T. Prosen, Open XXZ Spin Chain: NonequilibriumSteady State and a Strict Bound on Ballistic Transport,Phys. Rev. Lett. , 217206 (2011).[48] M. ˇZnidariˇc, T. Prosen, and P. Prelovek, Many-bodylocalization in the Heisenberg XXZ magnet in a randomfield, Phys. Rev. B , 064426 (2008).[49] J. Z. Imbrie, On Many-Body Localization for QuantumSpin Chains, J. Stat. Phys. , 998 (2016).[50] S. A. Weidinger and M. Knap, Floquet prethermaliza-tion and regimes of heating in a periodically driven, in-teracting quantum system, Sci. Rep. , 45382 (2017). [51] B. ˇZunkoviˇc, M. Heyl, M. Knap, and A. Silva, Dynam-ical Quantum Phase Transitions in Spin Chains withLong-Range Interactions: Merging Different Conceptsof Nonequilibrium Criticality, Phys. Rev. Lett. ,130601 (2018).[52] A. Greilich, A. Shabaev, D. Yakovlev, Al. L. Efros, I. A.Yugova, D. Reuter, A. D. Wieck, and M. Bayer, Nuclei-Induced Frequency Focusing of Electron Spin Coher-ence, Science , 1896 (2007).[53] I. Kleinjohann, E. Evers, P. Schering, A. Greilich, G.S. Uhrig, M. Bayer, and F. B. Anders, Magnetic fielddependency of electron spin revival amplitude in peri-odically pulsed quantum dots, Phys. Rev. B , 155318(2018).[54] R. Shindou, R. Matsumoto, S. Murakami, and J.-i. Ohe,Topological chiral magnonic edge mode in a magnoniccrystal, Phys. Rev. B , 174427 (2013).[55] K. Nakata, S. K. Kim, J. Klinovaja, and D. Loss,Magnonic topological insulators in antiferromagnets,Phys. Rev. B , 224414 (2017).[56] X. S. Wang, H. W. Zhang, and X. R. Wang, Topologicalmagnonics: A Paradigm for Spin-Wave Manipulationand Device Design, Phys. Rev. App. , 024029 (2018).[57] M. Malki and G. S. Uhrig, Topological magnon bandfor magnonics, Phys. Rev. B , 174412 (2019).[58] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera,M. R. Norman, and T. Senthil, Quantum Spin Liquids,Science , eaay0668 (2020).[59] S. Gao, O. Zaharko, V. Tsurkan, Y. Su, J. S. White,G. S. Tucker, B. Roessli, F. Bourdarot, R. Sibille, D.Chernyshov, T. Fennell, A. Loidl, and Ch. R¨uegg, Spiralspin-liquid and the emergence of a vortex-like state inMnSc S , Nature Phys. , 157 (2017).[60] S. A. D´ıaz, J. Klinovaja, and D. Loss, Topologi-cal Magnons and Edge States in AntiferromagneticSkyrmion Crystals, Phys. Rev. Lett. , 187203(2019).[61] P. Wadley, B. Howells, J. ˇZelezn´y, C. Andrews, V. Hills,R. P. Campion, V. Nov´ak, K. Olejn´ık, F. Maccherozzi,S. S. Dhesi, S. Y. Martin, T. Wagner, J. Wunderlich, F.Freimuth, Y. Mokrousov, J. Kuneˇs, J. S. Chauhan, M.J. Grzybowski, A. W. Rushforth, K. W. Edmonds, B.L. Gallagher, and T. Jungwirth, Electrical switching ofan antiferromagnet, Science , 587 (2016).[62] R. Lebrun, A. Ross, S. A. Bender, A. Qaiumzadeh, L.Baldrati, J. Cramer, A. Brataas, R. A. Duine, and M.Kl¨aui, Tunable long-distance spin transport in a crys-talline antiferromagnetic iron oxide, Nature , 222(2018).[63] H. Yu, S. D. Brechet, and J.-P. Ansermet, Spincaloritronics, origin and outlook, Phys. Lett. A ,825 (2017).[64] M. Gibertini, M. Koperski, A. F. Morpurgo, and K. S.Novoselov, Magnetic 2D materials and heterostructures,Nature Nanotech. , 408 (2019).[65] T. Li, A. Patz, L. Mouchliadis, J. Yan, T. A. Lograsso,I. E. Perakis, and J. Wang, Femtosecond switching ofmagnetism via strongly correlated spin-charge quantumexcitations, Nature , 69 (2013).[66] D. Bossini, S. Dal Conte, Y. Hashimoto, A. Secchi, R.V. Pisarev, Th. Rasing, G. Cerullo, and A. V. Kimel,Macrospin dynamics in antiferromagnets triggered bysub-20 femtosecond injection of nanomagnons, NatureCommun. , 10645 (2016). [67] D. Bossini, S. Dal Conte, G. Cerullo, O. Gomonay, R.V. Pisarev, M. Borovsak, D. Mihailovic, J. Sinova, J.H. Mentink, Th. Rasing, and A. V. Kimel, Laser-drivenquantum magnonics and terahertz dynamics of the or-der parameter in antiferromagnets, Phys. Rev. B ,024428 (2019).[68] M. J¨ackl, V. I. Belotelov, I. A. Akimov, I. V. Savochkin,D. R. Yakovlev, A. K. Zvezdin, and M. Bayer, Exci-tation of magnon accumulation by laser clocking as asource of long-range spin waves in transparent magneticfilms, Phys. Rev. X , 021009 (2017).[69] S. Manipatruni, D. E. Nikonov, C.-C. Lin, T. A. Gosavi,H. Liu, B. Prasad, Y.-L. Huang, E. Bonturim, R.Ramesh, and I. A. Young, Scalable energy-efficient mag-netoelectric spinorbit logic, Nature , 35 (2019).[70] T. S. Seifert, S. Jaiswal, J. Barker, S. T. Weber, I. Raz-dolski, J. Cramer, O. Gueckstock, S. F. Maehrlein, L.Nadvornik, S. Watanabe, C. Ciccarelli, A. Melnikov,G. Jakob, M. M¨unzenberg, S. T. B. Goennenwein, G.Woltersdorf, B. Rethfeld, P. W. Brouwer, M. Wolf, M.Kl¨aui, and T. Kampfrath, Femtosecond formation dy-namics of the spin Seebeck effect revealed by terahertzspectroscopy, Nature Commun. , 2899 (2018).[71] B. I. Shraiman and E. D. Siggia, Two-Particle Excita-tions in Antiferromagnetic Insulators, Phys. Rev. Lett.. , 740 (1988).[72] J. H. Mentink, K. Balzer, and M. Eckstein, Ultrafastand reversible control of the exchange interaction inMott insulators, Nat. Commun. , 6708 (2015).[73] T. F. Nova, A. Cartella, A. Cantaluppi, M. Frst, D.Bossini, R. V. Mikhaylovskiy, A. V. Kimel, R. Merlin,and A. Cavalleri, An effective magnetic field from opti-cally driven phonons, Nature Phys. , 132 (2017).[74] M. Fechner, A. Sukhov, L. Chotorlishvili, C. Kenel,J. Berakdar, and N. A. Spaldin, Magnetophononics:Ultrafast spin control through the lattice, Phys. Rev.Mater. , 064401 (2018).[75] G. Mazza and A. Georges, Superradiant Quantum Ma-terials, Phys. Rev. Lett. , 017401 (2019).[76] M. A. Sentef, M. Ruggenthaler, and A. Rubio, Cav-ity quantum-electrodynamical polaritonically enhancedelectron-phonon coupling and its influence on supercon-ductivity, Sci. Adv. , eaau6969 (2018).[77] S. Sachdev and R. Bhatt, Bond-operator representationof quantum spins: Mean-field theory of frustrated quan-tum Heisenberg antiferromagnets, Phys. Rev. B ,9323 (1990).[78] S. Gopalan, T. M. Rice, and M. Sigrist, Spin ladderswith spin gaps: A description of a class of cuprates,Phys. Rev. B , 8901 (1994).[79] M. Matsumoto, B. Normand, T. M. Rice, and M.Sigrist, Field- and pressure-induced magnetic quantumphase transitions in TlCuCl , Phys. Rev. B , 054423(2004).[80] K. P. Schmidt and G. S. Uhrig, Excitations in one-dimensional S = 1 / , 227204 (2003).[81] Ch. R¨uegg, B. Normand, M. Matsumoto, Ch. Nieder-mayer, A. Furrer, K. W. Kr¨amer, H. U. G¨udel, Ph.Bourges, Y. Sidis, and H. Mutka, Quantum Statistics ofInteracting Dimer Spin Systems, Phys. Rev. Lett. ,267201 (2005).[82] B. Normand and Ch. R¨uegg, Complete bond-operatortheory of the two-chain spin ladder, Phys. Rev. B , , 214417 (2015).[84] S. Diehl, A. Micheli, A. Kantian, H. P. B¨uchler andP. Zoller, Quantum states and phases in driven openquantum systems with cold atoms, Nature Phys. , 878(2008).[85] W. Berdanier, M. Kolodrubetz, R. Vasseur, and J. E.Moore, Floquet Dynamics of Boundary-Driven Systemsat Criticality, Phys. Rev. Lett. , 260602 (2017).[86] T. Kampfrath, K. Tanaka, and K. A. Nelson, Reso-nant and nonresonant control over matter and light byintense terahertz transients, Nature Photonics , 680(2013).[87] X. C. Zhang, A. Shkurinov, and Y. Zhang Extreme ter-ahertz science, Nature Photonics , 16 (2017).[88] P. Sal´en, M. Basini, S. Bonetti, J. Hebling, M. Krasil-nikov, A. Y. Nikitin, G. Shamuilov, Z. Tibai, V.Zhaunerchyk, and V. Goryashko, Matter manipulationwith extreme terahertz light: Progress in the enablingTHz technology, Phys. Rep. , 1 (2019).[89] C. Knetter and G. S. Uhrig, Perturbation theory by flowequations: dimerized and frustrated S = 1 / , 209 (2000).[90] C. Knetter, K. P. Schmidt and G. S. Uhrig, The struc-ture of operators in effective particle-conserving modelsJ. Phys. A: Math. Gen. , 7889 (2003).[91] H. Krull, N. A. Drescher, and G. S. Uhrig, Enhancedperturbative continuous unitary transformations, Phys.Rev. B , 125113 (2012).[92] M. Vogl, P. Laurell, A. D. Barr and G. A. Fiete, FlowEquation Approach to Periodically Driven QuantumSystems, Phys. Rev. X , 021037 (2019).[93] M. Hase, I. Terasaki, and K Uchinokura, Observationof the spin-Peierls transition in linear Cu (spin-1/2)chains in an inorganic compound CuGeO , Phys. Rev.Lett. , 3651 (1993).[94] R. Werner, C. Gros, and M. Braden, Microscopic spin-phonon coupling constants in CuGeO , Phys. Rev. B , 14356 (1999).[95] G. S. Uhrig, Symmetry and Dimension of the Dispersionof Inorganic Spin-Peierls Systems, Phys. Rev. Lett. ,163 (1997).[96] T. Lorenz, H. Kierspel, S. Kleefisch, B. B¨uchner, E.Gamper, A. Revcolevschi, and G. Dhalenne, SpecificHeat, Thermal Expansion, and Pressure Dependen-cies of the Transition Temperatures of Doped CuGeO ,Phys. Rev. B , R501 (1997).[97] C. Ecolivet, M. Saint-Paul, G. Dhalenne, and A.Revcolevschi, Brillouin scattering and ultrasonic mea-surements of the elastic constants of CuGeO , J. Phys.Condens. Matter, , 4157 (1999).[98] M. Hofmann, T. Lorenz, A. Freimuth, G. S. Uhrig, H.Kageyama, Y. Ueda, G. Dhalenne, and A. Revcolevschi,Heat transport in SrCu (BO ) and CuGeO , PhysicaB , 597 (2002).[99] P. Duthil, Materials Properties at Low Temperature, CERN Yellow Report , 77 (2014).[100] F. Giorgianni, B. Wehinger, S. Allenspach, N. Colonna,C. Vicario, P. Puphal, E. Pomjakushina, B. Normand,and Ch. R¨uegg, Nonlinear Quantum Magnetophononicsin SrCu (BO ) , unpublished.[101] D. A. Tennant, B. Lake, A. J. A. James, F. H. L. Essler,S. Notbohm, H.-J. Mikeska, J. Fielden, P. K¨ogerler, P.C. Canfield, and M. T. F. Telling, Anomalous dynamicalline shapes in a quantum magnet at finite temperature,Phys. Rev. B , 014402 (2012).[102] N. Ahmed, P. Khuntia, K. M. Ranjith, H. Rosner, M.Baenitz, A. A. Tsirlin, and R. Nath, Alternating spinchain compound AgVOAsO probed by As NMR,Phys. Rev. B , 224423 (2017).[103] U. Arjun, K. M. Ranjith, B. Koo, J. Sichelschmidt, Y.Skourski, M. Baenitz, A. A. Tsirlin, and R. Nath, Sin-glet ground state in the alternating spin-1/2 chain com-pound NaVOAsO , Phys. Rev. B , 014421 (2019).[104] I. S. Jacobs, J. W. Bray, H. R. Hart, L. V. Interrante, J.S. Kasper, G. D. Watkins, D. E. Prober, and J. C. Bon-ner, Spin-Peierls transitions in magnetic donor-acceptorcompounds of tetrathiafulvalene (TTF) with bisdithio-lene metal complexes, Phys. Rev. B , 3036 (1976).[105] M. C. Cross and D. S. Fisher, A new theory of the spin-Peierls transition with special relevance to the experi-ments on TTFCuBDT, Phys. Rev. B , 402 (1979).[106] S. Huizinga, J. Kommandeur, G. A. Sawatzky, B. T.Thole, K. Kopinga, J. M. de Jonge, and J. Roos, Spin-Peierls transition in N-methyl-N-ethylmorpholinium-ditetracyanoquinodimethanide [MEM-(TCNQ) ], Phys.Rev. B , 4723 (1979).[107] Y. Miura, R. Hirai, Y. Kobayashi, and M. Sato, Spin-Gap Behavior of Na Cu SbO with Distorted Honey-comb Structure, J. Phys. Soc. Jpn. , 084707 (2006).[108] M. B. Stone, W. Tian, M. D. Lumsden, G. E. Granroth,D. Mandrus, J.-H. Chung, N. Harrison, and S. E. Na-gler, Quantum Spin Correlations in an OrganometallicAlternating-Sign Chain, Phys. Rev. Lett. , 087204(2007).[109] G. S. Uhrig and B. Normand, Magnetic properties of(VO) P O from frustrated interchain coupling, Phys.Rev. B , R14705 (1998).[110] G. S. Uhrig and B. Normand, Magnetic properties of(VO) P O : two-plane structure and spin-phonon in-teractions, Phys. Rev. B , 134418 (2001).[111] F. A. Lindemann, ¨Uber die Berechnung molekularerEigenfrequenzen, Phys. Z. , 609 (1910).[112] J. J. Gilvarry, The Lindemann and Gr¨uneisen LawsPhys. Rev. , 308 (1956).[113] J. Dudowicz, K. F. Freed, and J. F. Douglas, Gener-alized entropy theory of polymer glass formation, Adv.Chem. Phys. , 125 (2008).[114] R. Gross and A. Marx, Festk¨orperphysik (Oldenbourg,Munich, 2012).[115] M. Born and J. R. Oppenheimer, Zur Quantentheorieder Molek¨ulen Ann. Phys.84