Pumping approximately integrable systems

Weak perturbations can drive an interacting many-particle system far from its initial equilibrium state if one is able to pump into degrees of freedom approximately protected by conservation laws. This concept has for example been used to realize Bose–Einstein condensates of photons, magnons and excitons. Integrable quantum systems, like the one-dimensional Heisenberg model, are characterized by an infinite set of conservation laws. Here, we develop a theory of weakly driven integrable systems and show that pumping can induce large spin or heat currents even in the presence of integrability breaking perturbations, since it activates local and quasi-local approximate conserved quantities. The resulting steady state is qualitatively captured by a truncated generalized Gibbs ensemble with Lagrange parameters that depend on the structure but not on the overall amplitude of perturbations nor the initial state. We suggest to use spin-chain materials driven by terahertz radiation to realize integrability-based spin and heat pumps.

A simple classical example for a weakly driven system is a well-insulated greenhouse. Due to the approximate conservation of the energy within the greenhouse, even weak sunlight can lead to high temperatures in its interior, which can be computed from the simple rate equation for the energy transfer. Similarly, large spin accumulation can be achieved in systems with approximate spin conservation 1 . Using approximate conservation of the number of photons, magnons or exciton polaritons, one can use pumping by light to reach densities, which allow for the realization of Bose-Einstein condensates [2][3][4] . Number-conserving collisions induce a quasi-equilibrium state in these systems, which can be efficiently described by introducing a chemical potential whose value is determined by balancing pumping and decay processes. Related theoretical approaches that describe electron-phonon systems far from equilibrium are so-called two-temperature models 5 : here one uses that the energy of the electrons and phonons are approximately separately conserved to introduce two different temperatures for the subsystems.
Integrable many-particle systems, like the one-dimensional (1D) fermionic Hubbard model or the XXZ Heisenberg model, are described by an infinite number of (local or quasi-local) conservation laws [6][7][8][9][10][11] . In closed integrable systems those prevent the equilibration into a simple thermal state, for example, after a sudden change of parameters. Instead the system can be described by a generalized Gibbs ensemble (GGE) [12][13][14][15][16][17][18][19][20][21] where C i are the conserved quantities and l i the corresponding Lagrange parameters. It has also been shown experimentally 22 that GGEs for a Lieb-Liniger model can provide highly accurate descriptions of interacting bosons in 1D. Many materials are described with high accuracy by integrable models 23 , however, weak integrability breaking terms and the coupling to thermal phonons imply that in equilibrium these systems are described by simple thermal states, r 0 $ e À bH , instead of GGEs. The proximity to the integrable point and the presence of approximate conservation laws leads to enhanced spin or heat conductivities (within linear-response theory) [24][25][26] and also to a slow relaxation after a quantum quench (via GGE-prethermalization) towards the equilibrium state 27 .
We will show that-as in the greenhouse example, see Fig. 1-such an approximately integrable system can be driven far from its thermal equilibrium by weak perturbations arising, for example, from a driving periodic in time or from coupling to a non-thermal bath. To balance the constant heating due to driving the system has to be weakly open, for example, by coupling to a phonon bath. As we will demonstrate this mechanism can be used for example to create large spin and heat currents. Besides the quasi-1D systems considered by us, also approximately many-body localized systems are characterized by infinitely many approximate conservation laws which may lead to a strong response to driving 28,29 .

Results
Weakly driven system. We consider an interacting many-body system that is approximately described by Hamiltonian H 0 and characterized by a finite or infinite number of (quasi-)local conserved quantities C i , [H 0 , C i ] ¼ 0, one of them being H 0 . Energy and other conservations are weakly broken by coupling to thermal or non-thermal baths and/or perturbations periodic in time. For simplicity, we assume periodic boundary conditions and a (discreet) translational invariance. We describe the system with density matrix r whose dynamics is governed by the Liouvillian super-operatorL, whereL can be split into the dominant unitary Hamiltonian evolutionL 0 r¼ À i H 0 ; r ½ and perturbationL 1 of strength E. We are interested in the limit of small E for t-N, where a unique (Floquet) steady state r N is obtained. The general structure of perturbation theory in this case has, for example, been discussed in refs 30-32. In this limit, r N can be approximated by r 0 ¼lim E!0 lim t!1 r withL 0 r 0 ¼0 according to equation (2). We assume and later support numerically that r 0 is approximately described by a GGE, see equation (1).
Here it is essential to note that-as in the greenhouse example discussed above-the parameters l i are not determined by the initial state but by the form of the weak perturbationsL 1 . Our central goal is to compute the l i . We first discuss the case of Lindblad dynamics, where perturbation theory linear in E can be used, and then focus on Hamiltonian dynamics where we have to consider E 2 contributions.
Markovian perturbation. Within the Markovian approximation one can use the Lindblad form forL 1 (ref. 33). Note that Lindblad dynamics is considered here mainly for pedagogical purposes (formulas are simpler) while no Lindblad approximation is used for the models studied below. The coefficients l i that fix the GGE are determined from the condition that the change of the approximately conserved quantities has to vanish in the steady state where we used thatL 0 r 0 ¼ À i H 0 ; r 0 ½ ¼0. Relation (3) yields a set of coupled equations for l i , where the number of equations is equal to the number of conserved quantities. We define the super-projectorP onto the tangential space of GGE density matrix,P using w ii 0 ¼ À Tr(C i @r 0 /@l i 0 ). Then the conditions for r 0 can be compactly written asL This equation can also be derived by considering higher order perturbations in E, see Methods for details.

Hamiltonian perturbation. For Hamiltonian dynamics
where H 1 may be a sum of several integrability breaking perturbations. Perturbation theory linear in E vanishes, TrðC i EL 1 r 0 Þ¼0 for all l i . Therefore one has to expand to order E 2 and equation (5) is replaced bŷ SincePðL 1 r 0 Þ¼0,L 1 r 0 is not in the kernel ofL À 1 0 . For periodic driving this equation has to be interpreted within the Floquet formalism, see Methods.
Model. As discussed in the introduction, our goal is to describe a situation which can be realized experimentally in spin-chain materials driven by lasers operating in the terahertz regime. We assume that spin chains are approximately described by a spin-1/2 XXZ Heisenberg model, possibly in the presence of an external magnetic field B, The system is driven out of equilibrium by a weak (integrability breaking) time-dependent perturbation with driving frequency o. This specific term has been chosen because it can induce heat and spin currents (as can be shown by a symmetry analysis), and because it can be realized experimentally. Such staggered exchange couplings and staggered magnetic fields arise naturally in certain compounds with (at least) two magnetic atoms per unit cell when coupled to uniform electric and magnetic fields, respectively [34][35][36][37] . See Fig. 1 for a schematic drawing of such a compound and Methods for concrete experimental suggestions. Therefore H d can be realized by shining a laser (typically at terahertz frequencies) onto the sample. In this case E 2 d is proportional to the laser power. Note that for T ¼ 0 and B ¼ 0 in the adiabatic limit, o-0, equations (7) and (8) realize an adiabatic Thouless pump, where per pumping cycle one spin is transported by one unit cell 38 . We will be interested in the opposite regime of large o and large (effective) temperatures.
Formally, the periodic perturbation H d would drive the system to infinite temperature [39][40][41][42] (up to remaining conservation laws 43 , possibly through a prethermal-like regime 44 ). In a solid state experiment this is prohibited by the coupling to phonons and, ultimately, to the thermal environment of the experimental set-up. We mimic this effect by coupling the spin system to a bath of Einstein phonons, H ph 0 ¼o ph P j a w j a j þ . . . , where dots stand for the couplings to further reservoirs which guarantee that the phonon system is kept at fixed temperature T ph , r ph $ e À H ph 0 =T ph . See Methods for details on finite size calculation using a broadened distribution of phonon energies. The (weak) coupling to the spin system is described by To obtain a unique steady state it is essential to break all symmetries, including the S z conservation. Relativistic effects which relax S z are mimicked by g m in our approach. We expect g m ( 1 in materials without heavy elements. For simplicity, we set g m ¼ 1 within our numerics as this is found to minimize finite size effects, without a qualitative influence on the results. Besides phonons also other integrability breaking perturbations exist in real materials, including defects, which typically dominate at the lowest temperatures. For high temperatures of the order of J (relevant for the considered set-up) it is realistic to assume that phonon coupling dominates.
In the presence of a periodic perturbation, equation (8), in the long-time limit the density matrix is changing periodically, r(t-N) ¼ P n e À iont r (n) with r ðnÞ w ¼r ð À nÞ , n 2 Z. Within the Floquet formalism one therefore promotes the steady-state density matrix to a vector and Liouville operator to a matrix, see Methods. For weak driving, E d -0, only the n ¼ 0 sector remains and the GGE ansatz, equation (1), simply reads r ðnÞ 0 ¼ðr 0 r ph Þd n;0 , where we included also the phonon density matrix, see above.
Steady state. We will use two different approaches to determine an approximate solution for the steady-state density matrix. First, we will parametrize r 0 , equation (1), with a small number of (quasi-)local conserved quantities, In an alternative approach, feasible for small systems, we take all conserved quantities into account: local and non-local, commuting and non-commuting. While the second approach is formally exact in the limit E d ,E ph -0, the first one is, perhaps, more intuitive and can be computed for larger system sizes.
For the XXZ Heisenberg model an infinite set of mutually commuting local conserved quantities C i is known, see Methods.
In addition there also exist (infinite) sets of quasi-local commuting conserved quantities [8][9][10] . As shown in refs 8,46 the spin-reversal parity-odd family has an overlap with the spin current J S at DoJ. Therefore both heat and spin current could show a large response to a weak perturbation. For our analysis, we choose three or five (N C ¼ 4, N C ¼ 6) most local conserved quantities C i , i ¼ 1, y, N C À 1. From the quasi-local sets, we include as a single (effective) operator the conserved part of spin current J S c , computed numerically 25,47 . For details see Methods. In the presence of an external magnetic field, equation (7), the heat current also has, in addition to C 3 , a spin current component, For the visualization of our results it is useful to define generalized forces F i in the space of Lagrange parameters by rewritingP _ r¼ P i @r 0 computed using exact diagonalization, see Methods. The vector F is a function of the Lagrange parameters l i , which points into the direction of the steady-state stable fixed point obtained from F i ¼ 0. In the absence of driving (Fig. 2a) one obtains the expected thermal state with T ¼ T ph while all other Lagrange parameters l i vanish. For finite driving the GGE is activated and the l i become finite (Fig. 2b). To obtain the steady state, we solve wF ¼ 0 using Newton's method. For the second approach, performed on small N-site systems, we first numerically construct a basis in the set of all (local and non-local) conserved operators, Q¼ n Due to degeneracies we find (for finite B and DaJ) about 2 Á 2 N elements Q i 2 Q. In the limit E d , E ph -0 the steady-state density matrix r N has to fulfil L 0 r 1 ¼ 0 and therefore can be exactly written as a linear combination of the Q i , Using equation (6), we therefore find that the steady-state density matrix for E d , E ph -0 is exactly given by the unique eigenvector with eigenvalue zero of the matrix In Fig. 3, we show the expectation value of the energy and of the heat current densities as functions of E d /E ph taking into account N C ¼ 4, N C ¼ 6, and all conserved quantities. The energy density expectation value is already obtained with good accuracy for N C ¼ 4 and even better for N C ¼ 6. The heat current vanishes both in thermal equilibrium, E d -0, and for E ph -0, where the system is described by an infinite temperature state with finite magnetization, r 0 $ e À l 1 S z and H 0 h i¼À B S z h i. It takes its largest value for E ph $ E d . For the currents a description in terms of N C ¼ 4 or 6 is qualitatively but not quantitatively accurate. Our study strongly suggests that further quasi-local conserved quantities contribute, as discussed in quench protocols [15][16][17] , see also ref. 25. For the chosen parameters our results depend only weakly on the system size N, see inset of Fig. 3. System size analysis is performed for N C ¼ 6 since the solution based on all conservations cannot be obtained for larger systems.
Our set-up can also be used to create spin currents. Whilst, by symmetry (bond-centered rotation in real and spin space by p around y axis), a finite external field B is needed to obtain a finite heat current, this is not the case for the spin current. Figure 4 displays the spin current density as a function of E d /E ph for B ¼ 0. Qualitatively one obtains a behaviour rather similar to the results for the heat current shown in Fig. 3 with a maximum in the spin current for E ph $ E d .
The external magnetic field B is a parameter which can easily be tuned experimentally. Figure 5 shows heat and spin current densities as a function of external magnetic field B for E d =E ph À Á 2 ¼2:5. Note that the sign of the magnetic field determines the sign of the heat current J H h i¼ C 3 h iÀ B J S h i. All main features of the B-dependence are semi-quantitatively reproduced by the truncated GGE with N C ¼ 6. For very large magnetic fields the convergence to the steady-state fixed point becomes slow as transitions rates connecting sectors with different magnetization are strongly suppressed, see Methods for further details.

Discussion
We have demonstrated that driving approximately integrable systems activates and pumps into approximately conserved quantities. Perhaps the most simple experimental set-up to measure the pumping effect predicted in this work, is to use a terahertz laser that excites a spin-chain material like Cu-benzoate where by symmetry staggered terms of the form (8) are expected 34,35 . As a consequence of the induced heat currents it is anticipated that the system cools down on one side while it heats upon the other. The direction of the effect can be controlled  For vanishing magnetic field a spin current (but no heat current) is generated within our model for finite ratios of E d /E ph . The expectation value of spin current density is again maximal for E d /E ph E1. Parameters: Effective force F in the space of Lagrange parameters (b, l 3 ) using e À bH0 À l3C3 as an ansatz for the generalized Gibbs ensemble. Parameters: For the chosen parameters the spin and heat currents expressed in dimensionless units appear to be rather small, of the order of 10 À 3 . While these values can definitely be increased by tuning parameters, for example the external magnetic field, it is important to note that the currents are actually quite large compared to the typical heat or spin currents obtained in bulk materials. To create a heat current of similar size in a good heat conductor like Cu (assuming J $ k B Á 100 K, 5 Å for the distance of the spin chains, and k Cu E400 Wm À 1 K À 1 ) one would need a temperature gradient of several 10 5 Km À 1 . Similarly, to create a (transversal) spin current of comparable size in a heavy element like Pt using the spin Hall effect (assuming r Pt E10 mO cm and a Pt s % 10% for the spin Hall angle 48 ) one needs electric fields of the order of 10 4 V m À 1 or sizable current densities of the order of 10 11 Am À 2 . These numbers are even more remarkable when one takes into account that the electron densities in Cu or Pt are at least an order of magnitude higher than the spin density for spin chains with a distance of 5 Å.
While our study has focused on the steady state, it is instructive to discuss the relevant timescales for its buildup. For this argument, we consider a quench where at time t ¼ 0 an initial state is perturbed both by the integrable part of the Hamiltonian and by small non-integrable perturbations. At short times of the order of several 1/J the initial state will prethermalize 27,49-52 into a GGE, where the values of the conserved quantities, C i h i, are set by the initial conditions (with small corrections from the perturbations 52,53 ). Further time evolution can be approximately described by a GGE with time-dependent Lagrange parameters. Their time-dependence is determined by perturbations which assert forces F i $ E 2 , such that dl i /dtEF i . Governed by the perturbations the system will loose the memory of its initial condition on a timescale of order 1/E 2 and relax to the steady state (obtained from F i ¼ 0) which is, in general, completely unrelated to the prethermalized state. Note that the same approach predicts ordinary thermalization in the absence of external driving. Our results suggest that the concept of generalized GGEs has a much broader range of application than previously anticipated, now extended to open systems where symmetries are not exact and integrability is weakly broken. A truncated GGE proved to be useful for qualitative description, however, it showed quantitative discrepancies most probably due to disregarded quasi-local conserved quantities, as observed already in quench protocols 15,16 . We are planning a future study tailored to address this issue systematically. It would be interesting to develop integrability-based methods similar to the quench-action approach 15,16,54,55 to treat such situations.
Most important for applications is that the integrability is not required to be realized exactly but only approximately. Efficient pumping requires only that the pumping rates are of the same order of magnitude as the loss rates arising from integrability breaking terms. Especially the integrability-based creation of large spin currents could find its application in future spintronics devices.

Methods
Perturbing around q 0 . The central equations (5) or (6), used to determine the density matrix r 0 in the limit E-0, have to be consistent and can also even be derived by considering perturbations around r 0 , r N ¼ r 0 þ dr.
First, the leading dr correction to _ C i , equation (3), arising from Tr C iL0 dr which is nominally of the same order as Tr C i EL 1 r 0 vanishes trivially as For arbitrary r 0 , dr is exactly given by dr¼ ÀL À 1 EL 1 r 0 , whereL À 1 is a short-hand notation for lim Z!0 ðL À Z1Þ À 1 with the infinitesimal regularizer Z. The correct expansion point r 0 is found if lim E!0 dr¼0. Below we show that for the projection operatorP, equation (4), which would yieldL À 1P EL 1 r 0 $ Oð1Þ. This contradicts our perturbative approach unlessPL 1 r 0 ¼0, as set by our condition equation (5). Equation (12) is a consequence of the fact thatP projects onto the tangential space to GGE density matrix. In this spaceL 0 vanishes by definition, L 0 @r 0 =@l i ð Þ ¼ 0, andL¼L 0 þ EL 1 is therefore of order E. Technically, this can be seen by using the general relation forX ¼PEL 1P ; withQ¼1 ÀP andX þŶ¼L. Then The second term is O(1) asL 0P ¼0 and thereforeŶP $ OðEÞ. The divergence ofL À 1P for E-0 can be directly related to the fact that integrable systems are characterized by infinite conductivities (finite Drude weights) at finite temperatures 56 as can, for example, be seen 24 within the memory matrix formalism 57 . All arguments given above can be generalized to situations where leading corrections arise from second-order perturbation theory in which case one obtains equation (6) instead of equation (5). Staggered hopping and magnetic field modulation. Sizable staggered g-tensors leading to staggered B-fields have been observed in a number of different compounds [34][35][36][37] . Similarly an external electric field will distort the crystalline structure in these materials, leading to staggered exchange couplings linear in homogeneous electric fields. An example of such a material is Cu-benzoate 34 with the above modulations allowed by symmetry for electric (magnetic) fields applied in the 010 (001) crystallographic direction. In this system the staggered g-tensor has been measured to be B0.08 (ref. 35), the size of the staggered exchange coupling is unknown. For simplicity, we assume in equation (8) that the two staggered terms are of the same size.
Conservation laws of the XXZ Heisenberg model. An infinite set of local conserved quantities C i of the Heisenberg model . In general, C i are operators involving maximally i neighbouring sites. Importantly, C 3 in the absence of external magnetic field equals the heat current with rescaled spin operators S 0a In the presence of external magnetic field, equation (7), heat current has in addition to C 3 also a spin current component, As understood recently there also exist families of quasi-local conserved quantities [8][9][10] , which are mostly disregarded in our study with the exception of a spin-reversal parity-odd operator, J S c . The latter is constructed as the conserved part of the spin current operator J S , whereñ j i are simultaneous eigenstates of the C i . Since it is known that the spin current has an overlap with the quasi-local family 45 for DoJ, the conserved J S c contains quasi-local components (and, possibly, non-local components not contributing in the thermodynamic limit).
Floquet formulation. For a periodically driven system described by _ r¼LðtÞr witĥ Lðt þ TÞ¼LðtÞ the density matrix changes periodically in the long-time limit. Therefore it is useful to split it into Floquet components, r¼ X n e À inot r ðnÞ ; n 2 Z ð19Þ with r ð À nÞ ¼r ðnÞ w and o ¼ 2p/T. The Floquet components are combined into the vector r¼ . . . r ð À 1Þ ; r ð0Þ ; r ð1Þ ; . . . À Á . The Liouvillian is promoted to a (static) matrixL nm ¼inod nm þL n À m withL n À m ¼ 1 T R T 0L ðtÞe ioðn À mÞt dt. Using this notation, all results obtained for static Liouvillian super-operators directly translate to the time-periodic case. Within our set-up H 0 , all approximate conservation laws C i and the GGE density matrix r 0 are static and therefore the projection operatorP, equation (4), projects onto the n ¼ 0 Floquet sector only. The steady-state condition, equation (6), thus means that the approximately conserved quantities do not grow after averaging over an oscillation period. To second order in E d only transitions from the n ¼ 0 to the n ¼ ±1 Floquet sector and back contribute to equations (6) or (10) asL n ¼0 for n j j41. For the generalized force due to the periodic driving, we obtain from (10) where we used H 0 eigenstates m

Þ.
Note that equation (20) contains-as expected-transition rates well-known from Fermi's golden rule. Equation (20) is evaluated for finite systems of size N by replacing the d function by a Lorentzian (1/p)Z/(o 2 þ Z 2 ) (Z ¼ 0.1J for N ¼ 12).
Equation (20) is only valid for situations where all conservation laws commute with each other, with C i ¼ P m m j iC i;m m h j, see below for a brief discussion of the non-commuting case.
Phonon coupling. As written in the main text, we assume that the phonon system always remains at equilibrium, r ph $ e À H ph 0 =Tph . Using equation (10), after tracing over phonons, we obtain for the generalized force where n B ðEÞ¼1=ðe E=Tph À 1Þ is the equilibrium Bose distribution evaluated at the temperature T ph and A (ph) (o) is the phonon spectral function. For our finite size calculation, we broaden the spectral function of the Einstein phonons using A ðphÞ ðoÞ¼YðoÞ o ophZ ffiffi p p e À ðo À ophÞ 2 =Z 2 . This choice of broadening ensures detailed balance relations (necessary to obtain a thermal state in the absence of driving) and the positivity of phonon frequencies (necessary for stability). For all plots we use Z ¼ 0.4J. However, we have checked that similar results are obtained, for example, for Z ¼ 0.1J for magnetic fields up to B j j¼ 2J. For larger fields Z ¼ 0.1J does not provide a sufficient amount of relaxation between sectors with different magnetization and convergence becomes slow and unstable. For Z ¼ 0.4J larger fields, B j jt5J, can be reached.
Implementation of non-commuting conservation laws. As discussed in the main text, a complete basis of all non-local commuting or non-commuting conserved quantities is given by Q¼ n j i m h jwith E 0 m ¼E 0 n È É which solve the equationL 0 Q i ¼0 for Q i 2 Q. Using the exact eigenstates of H 0 it is straightforward to evaluate equation (11), where we use for our finite size calculations the broadening procedures described above. As a technical detail, we note that, when one follows this procedure, one has to evaluate in the phonon sectors integrals of the type ð Þdo 0 numerically. For efficient evaluations, we use interpolating functions for these integrals. GGE estimation for other conserved quantities. To provide further support for our claim that truncated GGEs give a semi-quantitative description of our weakly open system we show in Fig. 6 additional comparison of the H 0 h i and C 4 h i as a function of magnetic field B at E d =E ph À Á 2 ¼2:5, comparing as in the main text the exact calculation including all conserved quantities and the truncated GGE with N C ¼ 6 (quasi-)local conserved quantities. The GGE ansatz captures the right magnitude and the correct behaviour in the dependence on B also for more complicated 4-spin operators like C 4 . We use same parameters as for the Fig. 5 in the main text: E d =E ph À Á 2 ¼2:5, J ¼ 1, Code availability. Custom computer codes used in this study are available from the corresponding author upon request. Documentation of the codes is not available.
Data availability. Data is available from the corresponding author upon reasonable request.