Perturbative approach to weakly driven many-particle systems in the presence of approximate conservation laws

TL;DR

This paper introduces a perturbative method to analyze weakly driven quantum many-particle systems with approximate conservation laws, enabling efficient description of their steady states via generalized Gibbs ensembles.

Contribution

It develops a Liouville perturbation theory that accounts for conservation laws and their violations, providing a systematic way to find steady states in weakly driven quantum systems.

Findings

Perturbative expansion of steady states using Lagrange multipliers.

Application to models like electron-phonon systems and Bose condensates.

Validation through study of interacting fermions with non-thermal reservoirs.

Abstract

We develop a Liouville perturbation theory for weakly driven and weakly open quantum systems in situations when the unperturbed system has a number of conservations laws. If the perturbation violates the conservation laws, it drives the system to a new steady state which can be approximately but efficiently described by a (generalized) Gibbs ensemble characterized by one Lagrange parameter for each conservation law. The value of those has to be determined from rate equations for conserved quantities. Remarkably, even weak perturbations can lead to large responses of conserved quantities. We present a perturbative expansion of the steady state density matrix; first we give the condition that fixes the zeroth order expression (Lagrange parameters) and then determine the higher order corrections via projections of the Liouvillian. The formalism can be applied to a wide range of problems including two-temperature models for electron-phonon systems, Bose condensates of excitons or photons or weakly perturbed integrable models. We test our formalism by studying interacting fermions coupled to non-thermal reservoirs, approximately described by a Boltzmann equation.

Figures

Tables

Table 1. Table 1: Comparison of semi-classical open Boltzmann dynamics for level occupation function and quantum Liovillian formulation for density matrix

Boltzmann	Liouvillian
occupation function $f_e$	density matrix $\rho$
$\frac{d{f}_e}{dt}=M{\left[f\right]}_e+ϵD{\left[f\right]}_e$	$\frac{d\rho }{dt}={\widehat{ℒ}}_0\rho +{\widehat{ℒ}}_1\rho$
$f_e\left(t\to \infty \right)=f_e^0+\delta f_e$	$\rho \left(t\to \infty \right)={\rho }_0+\delta \rho$
Fermi function: $f_e^0=\frac{1}{1+e^{\beta \left(e-\mu \right)}}$	GGE: ${\rho }_0=\frac{e^{-{\lambda }_i{C}_i}}{\text{Tr}\left[e^{-{\lambda }_i{C}_i}\right]}$
conservation laws and scalar product
$c_i\left(e\right)$ and $q_i^{\beta ,\mu }$	$C_i$ and $\partial {\rho }_0/\partial {\lambda }_i$
$\int 𝑑e{c}_i\left(e\right)f_e$	$\text{Tr}\left[C_i\rho \right]=⟨C_i⟩$
${\chi }_{ij}(\beta ,\mu )=-\int 𝑑e{c}_i\left(e\right)q_j^{\beta ,\mu }\left(e\right)$	${\chi }_{ij}=-\text{Tr}\left[C_i\partial {\rho }_0/\partial {\lambda }_j\right]$
zeroth order perturbation theory
$\int 𝑑e{c}_i\left(e\right)ϵD{\left[f^0\right]}_e=0$	$\text{Tr}\left[C_i{\widehat{ℒ}}_1{\rho }_0\right]=0$
$\hat{P}\left[D\left[f^0\right]\right]=0$	$\hat{P}\left({\widehat{ℒ}}_1{\rho }_0\right)=0$
first order corrections
${𝒇}_{⟂}^1=-{\left(Q{M}^{\left(0\right)}Q\right)}^{-1}Q𝑫\left[{𝒇}^0\right]$	$\delta {\rho }_{1,⟂}=-{\left(\hat{Q}{\widehat{ℒ}}_0\hat{Q}\right)}^{-1}\hat{Q}{\widehat{ℒ}}_1{\rho }_0$
${𝒇}_{\parallel }^1={\left(P{D}^{\left(1\right)}P\right)}^{-1}P{D}^{\left(1\right)}Q\times$	$\delta {\rho }_{1,\parallel }={\left(\hat{P}{\widehat{ℒ}}_1\hat{P}\right)}^{-1}\hat{P}{\widehat{ℒ}}_1\hat{Q}\times$
$\times {\left(Q{M}^{\left(0\right)}Q\right)}^{-1}Q𝑫\left[{𝒇}^0\right]$	$\times {\widehat{ℒ}}_0^{-1}\hat{Q}{\widehat{ℒ}}_1{\rho }_0$

Equations

$$\displaystyle\hat{\mathcal{L}}_{0}\rho=-i[H_{0},\rho],\quad\hat{\mathcal{L}}_{1}\rho=\left\{\begin{array}[]{ll}-i[\epsilon H_{1},\rho],\\ \epsilon\hat{\mathcal{D}}\rho.\end{array}\right.$$

$$\hat{\mathcal{D}}\rho=\sum_{\alpha}\big{(}L_{\alpha}\rho L_{\alpha}^{\dagger}-\frac{1}{2}\{L_{\alpha}^{\dagger}L_{\alpha},\rho\}\big{)}$$

$$\underset{t\to \infty }{\lim }\rho (t)={\rho }_{\infty }={\rho }_0+\delta \rho .$$

$$\delta \rho =-{\hat{\mathcal{L}}}^{-1}{\hat{\mathcal{L}}}_1{\rho }_0,{\hat{\mathcal{L}}}^{-1}\to \underset{\eta \to 0}{\lim }(\hat{\mathcal{L}}-\eta \hat{1}{)}^{-1}.$$

$${\rho }_0=\frac{e^{-{\sum }_i{\lambda }_i{C}_i}}{\text{Tr}[e^{-{\sum }_i{\lambda }_i{C}_i}]},$$

$$⟨{\dot{C}}_i⟩$$

$$=\text{Tr}[C_i{\hat{\mathcal{L}}}_1{\rho }_0]+\text{Tr}[C_i{\hat{\mathcal{L}}}_1\delta \rho ]\approx \text{Tr}[C_i{\hat{\mathcal{L}}}_1{\rho }_0]$$

$$\text{Tr}[C_i{\hat{\mathcal{L}}}_1{\rho }_0]\stackrel{!}{=}0$$

$$(X+Y{)}^{-1}=X^{-1}-(X+Y{)}^{-1}Y{X}^{-1},$$

$$\text{Tr}[C_i{\hat{\mathcal{L}}}_1{\hat{\mathcal{L}}}_0^{-1}{\hat{\mathcal{L}}}_1{\rho }_0]\stackrel{!}{=}0.$$

$$\rho (t)=\sum _n{e}^{-in\omega t}{\rho }^{(n)},n\in \mathbb{Z},$$

$$\hat{\mathcal{L}}$$

$${\hat{\mathcal{L}}}_1^{nm}$$

$${\hat{\mathcal{L}}}_0^{nm}$$

$$\hat{P}X$$

$$\hat{Q}X$$

$${\hat{P}}^{†}A=-\sum _{ij}{C}_i({\chi }^{-1}{)}_{ji}\text{Tr}\left[A\frac{\partial {\rho }_0}{\partial {\lambda }_j}\right].$$

$$\hat{P}\dot{\rho }$$

$${\dot{\lambda }}_i$$

$$\hat{P}({\hat{\mathcal{L}}}_1{\rho }_0)=0.$$

$$\hat{P}({\hat{\mathcal{L}}}_1{\hat{\mathcal{L}}}_0^{-1}{\hat{\mathcal{L}}}_1{\rho }_0)=0.$$

$$\mathcal{Q}=\{|n⟩⟨m|\mathrm{with}{E}_n^0=E_m^0\}.$$

$${\hat{\mathcal{L}}}_{nm}^{\mathcal{Q}}=\text{Tr}[{\mathcal{Q}}_n^{†}{\hat{\mathcal{L}}}_1{\mathcal{Q}}_m].$$

$${\hat{\mathcal{L}}}_{nm}^{\mathcal{Q}}=-\text{Tr}[{\mathcal{Q}}_n^{†}{\hat{\mathcal{L}}}_1{\hat{\mathcal{L}}}_0^{-1}{\hat{\mathcal{L}}}_1{\mathcal{Q}}_m].$$

$$\tilde{P}X=\sum _{n,m}⟨n|X|m⟩{\delta }_{E_n^0,E_m^0}|n⟩⟨m|$$

$$\displaystyle\langle C_{i}\rangle=\text{Tr}\big{(}C_{i}(\rho_{0}+\delta\rho)\big{)}=\text{Tr}\big{(}C_{i}(\rho_{0}+\delta\rho_{\parallel})\big{)}$$

$$\delta \rho =-{\hat{\mathcal{L}}}^{-1}\hat{Q}{\hat{\mathcal{L}}}_1{\rho }_0.$$

$${\hat{\mathcal{L}}}^{-1}\hat{Q}=({\hat{\mathcal{L}}}_0+{\hat{\mathcal{L}}}_1{)}^{-1}\hat{Q}=({\hat{\mathcal{L}}}_0+\hat{P}{\hat{\mathcal{L}}}_1\hat{P}{)}^{-1}$$

$$\displaystyle\times\sum_{n=0}^{\infty}\big{[}-(\hat{P}\hat{\mathcal{L}}_{1}\hat{Q}+\hat{Q}\hat{\mathcal{L}}_{1}\hat{P}+\hat{Q}\hat{\mathcal{L}}_{1}\hat{Q})(\hat{\mathcal{L}}_{0}+\hat{P}\hat{\mathcal{L}}_{1}\hat{P})^{-1}\big{]}^{n}\hat{Q}.$$

$$({\hat{\mathcal{L}}}_0+\hat{P}{\hat{\mathcal{L}}}_1\hat{P}{)}^{-1}\hat{P}$$

$$({\hat{\mathcal{L}}}_0+\hat{P}{\hat{\mathcal{L}}}_1\hat{P}{)}^{-1}\hat{Q}$$

$$\delta \rho$$

$$\delta {\rho }_{1,\parallel }$$

$$\delta {\rho }_{1,\perp }$$

$$(\hat{P}{\hat{\mathcal{L}}}_1\hat{P}{)}^{-1}\hat{P}X$$

Full text

In equilibrium many particle systems can be efficiently described by Gibbs ensembles, characterized by one Lagrange parameter (inverse temperature, chemical potential) for each exactly conserved quantity (energy, particle number). Such a simple but powerful description is, in general, absent for driven non-equilibrium systems. As long as weak perturbations, which drive the system out of equilibrium, have small effects one can resort to perturbation theory. Kubo formulas, for example, describe the response to time-dependent Hamiltonian perturbations. Also for open systems described by Lindblad operators similar formulations of perturbation theory can be developed benatti11 ; li14 ; li16 ; maizelis14 ; garcia09 ; kessler12 ; znidaric15 ; Reiter12 ; lemos17 . If the density matrix of the unperturbed system is not unique (as is the case for all Hamiltonian systems) one has to use degenerate Liouville perturbation theory benatti11 ; znidaric15 .

There are many situations where even weak perturbations of interacting many-particle systems can have large effects. A famous example is the Bose-Einstein condensation of excitons and polaritons butov01 ; kasprzak06 ; balili07 : irradiation by light creates more and more excitons which equilibrate approximately and form a Bose-Einstein condensate. Importantly, the exciton number is approximately conserved thus even weak pumping can compensate for exciton losses and leads to large densities of these excitations. Similarly, condensations of photons or magnons have been observed klaers10 ; klaers12 ; schmitt14 ; demokritov06 ; chumak09 .

This is a general phenomenon: whenever approximate conservation laws exist, small perturbations which weakly break those conservation laws can drive the system out of equilibrium, to a steady state completely different from the initial one. In the example given above, the approximately conserved quantity is the exciton number. Other examples use the approximate conservation of spin to induce large nuclear spin polarization walker97 for medical applications or use the weak coupling of electrons and phonons to induce vastly different temperatures in the two subsystems anisimov74 ; armstrong81 ; eesley86 ; allen87 ; singh10 . In all cases, the long-time steady state – a Bose-Einstein condensate for the exciton example – is very different from any thermal state. Nevertheless, simple approximate theoretical descriptions can be found, e.g., by introducing an effective chemical potential for the excitons.

As we have argued recently in Ref. lange17 , this approach works quite generally: generalized Gibbs ensembles characterized by Lagrange parameters for approximately conserved quantities give an approximate description of the steady state of weakly perturbed driven systems. They describe situations where the driven system is completely different from the initial one, with Lagrange parameters determined by the perturbations. The nearly conserved quantities – e.g. the exciton number in the example discussed above – show a strong (non-linear) response to weak perturbations. In this paper we develop a perturbation theory around these generalized Gibbs ensembles and show that deviations remain small as long as the perturbations driving the system out of equilibrium are small. This is, however, only true if the correct zeroth order approximation, i.e. the correct generalized Gibbs ensemble, is chosen as the reference point of the perturbative expansion. The formalism is developed both for Lindblad and unitary time-periodic (Floquet) perturbations.

From this perspective our investigation is related to other Floquet studies in the presence of interactions and/or coupling to the environment dehghani14 ; shirai15 ; iadecola15 ; liu15 ; iwahori16 ; shirai16 ; seetharam15 ; iwahori17 ; genske15 . For example, in Ref. genske15 experiments of the Esslinger group on a Floquet-realization of the Haldane model esslinger have been modeled based on the derivation of a Floquet-Boltzmann equation genske15 ; seetharam15 . The concept of a (time-dependent) temperature was used to describe the heating of interacting Floquet systems for situations where energy-conserving scattering processes dominate scattering events, where the energy changes. A related approach is taken when describing weakly-coupled electron-phonon systems out of equilibrium. Here one often associated two different temperatures to the two subsystems with energies that are separately approximately conserved anisimov74 ; armstrong81 ; eesley86 ; allen87 ; singh10 ; sentef13 . In Ref. lange17 we used generalized Gibbs ensembles to describe the steady state of spin-chains coupled to phonons and driven by an external laser, modeled by a periodic driving. In these systems an infinite number of approximately conserved quantities exists due to the integrability of the underlying spin-1/2 Heisenberg model. As the heat current is also approximately conserved, one can realize large heat currents by weak driving.

From a more general point of view, our theoretical approach is based on the idea to project the dynamics of the density matrix on a few relevant degrees of freedom – in our case the approximate conservation laws. Such projection formalisms are widely used by many authors Robertson66 ; Robertson68 ; kawasaki73 ; Ochiai73 ; Willis74 ; Grabert74 ; Koide02 ; breuer02 . In contrast to these studies, we will use generalized Gibbs ensembles as a reference point and we focus our calculation on the computation of the steady state.

The article is organized in the following way: Sec. I prepares the perturbative setup, Sec. II gives the prescription on how to find the zeroth order approximation to the steady state via the introduction of a super-projector. The latter is used in Sec. III to determine corrections around the approximate description of the steady state both for Lindblad and unitary time-periodic (Floquet) perturbations. The applicability of our degenerate perturbation theory is tested in Sec. IV on an example of interacting fermions coupled to non-thermal reservoirs which induce weak gain and loss of particles. We first perform an exact calculation on small system and then use approximate description by Boltzmann dynamics in the thermodynamic limit.

I Model

We consider a many-particle quantum system described by the density operator $\rho$. Its dynamics is determined from the Liouville equation $\dot{\rho}=\hat{\mathcal{L}}\rho$, where $\hat{\mathcal{L}}$ is the Liouvillian super-operator. We consider situations where a Hamiltonian system is weakly perturbed. Therefore we split $\hat{\mathcal{L}}=\hat{\mathcal{L}}_{0}+\hat{\mathcal{L}}_{1}$ into two parts, where $\hat{\mathcal{L}}_{0}$ describes the dominant unitary Hamiltonian evolution and $\hat{\mathcal{L}}_{1}$ a weak perturbation of strength $\epsilon$ which drives the system out of equilibrium,

[TABLE]

Most importantly, the unperturbed system described by $H_{0}$ has the property of having additional conserved quantities $C_{i}$, $[H_{0},C_{i}]=0$, $i=1,..,N_{\textrm{c}}$. We are interested both in cases where $N_{\textrm{c}}$ is small or when $N_{c}$ is extensive as is the case when $H_{0}$ is integrable. The perturbation $\hat{\mathcal{L}}_{1}$ is of the unitary or/and of the Markovian Lindblad form. In $H_{1}$ we consider both static and time-dependent perturbations which are periodic in time. The Lindblad dynamics is described by

[TABLE]

with Lindblad operators $L_{\alpha}$. We consider only translational invariant perturbations and situations where a unique (Floquet-) steady state is obtained. We assume that the perturbations break all conservation laws considered above. Remaining exactly conserved quantities are fixed by initial conditions and can easily be included in the theoretical description but are omitted in the following to simplify notations.

We are mainly interested in the nonequilibrium steady state. We split its density operator $\rho_{\infty}$ into the zeroth order approximation $\rho_{0}$ and corrections $\delta\rho$,

[TABLE]

with $\rho_{0}=\lim_{\epsilon\to 0}\lim_{t\to\infty}\rho(t)$. Note that the limits $t\to\infty$ and $\epsilon\to 0$ do not commute: as we will show, small perturbations can completely change the density matrix in the long-time limit.

For perturbations periodic in time, the density matrix is periodically oscillating in the long time limit but one can still use the above formulas by interpreting them in Floquet space, see Section II.2. The higher order corrections $\delta\rho$ are formally given by

[TABLE]

where we used $\hat{\mathcal{L}}\rho_{\infty}=\hat{\mathcal{L}}_{1}\rho_{0}+\hat{\mathcal{L}}\delta\rho=0$ and $\hat{\mathcal{L}}_{0}\rho_{0}=0$. The inverse $\hat{\mathcal{L}}^{-1}$ should be interpreted using infinitesimal regularization $\eta$, see App. A. Due to the conservation laws of $H_{0}$, $\hat{\mathcal{L}}_{0}\rho_{0}=0$ has no unique solution. While Eq. (6) is formally valid for arbitrary $\rho_{0}$ with $\hat{\mathcal{L}}_{0}\rho_{0}=0$, the correction $\delta\rho$ will only be small if $\rho_{0}$ is correctly chosen as discussed in the next section.

II Zeroth order: generalized Gibbs ensemble

We have defined $\rho_{0}$ to fulfill the equation $\hat{\mathcal{L}}_{0}\rho_{0}=0$. In the thermodynamic limit, generic steady states of interacting many-particle system with (quasi-local) conservation laws $C_{i}$ approach states which can be described by a (generalized) Gibbs ensemble (GGE)

[TABLE]

More precisely, the system as a whole may be in a different ensemble (e.g. a canonical ensemble) but for the computation of local observables one expects that the different ensembles are generically equivalent in the thermodynamic limit. One manifestation of this is the eigenstate thermalization hypothesis polkovnikov11 ; nandkishore15 which argues that a generic pure state becomes equivalent to a Gibbs state in the thermodynamic limit. Eq. (7) has also been extensively tested for integrable models rigol07 ; calabrese07gge ; cramer08gge ; barthel08gge ; essler12gge ; pozsgay13gge ; fagotti13gge ; fagotti13gge2 ; fagotti14gge ; fagotti14agge ; wouters14gge ; pozsgay14gge ; sotiriadis14gge ; goldstein14gge ; brockmann14gge ; pozsgay14agge ; rigol14gge ; mestyan15gge ; ilievski15agge ; eisert15gge ; ilievski16gge ; piroli16gge ; calabrese11gge ; calabrese12gge ; calabrese12agge ; fagotti13gge ; caux13gge ; essler16gge ; essler17gge ; ilievski16 which appear to approach a state described by $\rho_{0}$ after a quench. For certain local initial states (satisfying cluster decomposition property) this has been shown rigorously ilievski15agge . It has been argued recently that there exist a protocol based on truncated GGEs, convergent in number of included conservation laws pozsgay17 .

Note that the parameters $\lambda_{i}$ in Eq. (7) have not yet been determined. Most importantly, they are not fixed by initial conditions as the $C_{i}$ are not conserved in the presence of the perturbations described by $\hat{\mathcal{L}}_{1}$. Instead the $\lambda_{i}$ or, equivalently, the expectation values of the $C_{i}$ have to be determined from rate equations governed by the weak perturbations. As we will show next, these will lead to changes of the $\lambda_{i}$ of order $\mathcal{O}(1)$.

II.1 Determination of $\lambda_{i}$

All GGEs satisfy $\hat{\mathcal{L}}_{0}\rho_{0}=0$ by definition. Therefore the perturbation $\hat{\mathcal{L}}_{1}$ fixes the steady state parameters $\{\lambda_{i}\}$. Technically, we determine the $\{\lambda_{i}\}$ from the condition that the state is stationary in combination with the condition that $\delta\rho$ should be small in the limit $\epsilon\to 0$. The first condition ensures that $\langle\dot{C}_{i}\rangle=0$, which is evaluated using straightforward perturbation theory

[TABLE]

for $i=1,...,N_{\textrm{c}}$. Above we used that $\text{Tr}[C_{i}\hat{\mathcal{L}}_{0}\delta\rho]=0$ due to $\hat{\mathcal{L}}_{0}^{\dagger}C{}_{i}=i[H_{0},C_{i}]=0$ and that $\delta\rho$ is small. Note that the adjoined of a Liouvillian is defined by the equation $\text{Tr}[A\hat{\mathcal{L}}\rho]=\text{Tr}[(\hat{\mathcal{L}}^{\dagger}A)\rho]$. We therefore obtain from $\langle\dot{C}_{i}\rangle=0$ the condition fixing $\rho_{0}$

[TABLE]

If $\text{Tr}[C_{i}\hat{\mathcal{L}}_{1}\rho_{0}]=0$ is trivially fulfilled for all $\rho_{0}$, one has to consider higher-order perturbation theory, see below.

Assuming for the moment that $\text{Tr}[C_{i}\hat{\mathcal{L}}_{1}\rho_{0}]\neq 0$ for all $i=1,...,N_{\textrm{c}}$ and generic $\rho_{0}$, then the $N_{\textrm{c}}$ equations (9) can be used to fix the $N_{\textrm{c}}$ Lagrange parameters $\lambda_{i}$. These equations can be viewed as rate equations, which describe how the perturbations change the approximately conserved $C_{i}$. Within the approximation used above, the size of the perturbation, i.e. the value of $\epsilon$, completely drops out of the equation for the steady state. Therefore it fixes the Lagrange parameters (and, accordingly, also $\rho_{0}$) to order $\epsilon^{0}$. While to this order the size of the perturbation is not important, its structure determines the Lagrange parameters and induces changes of order $1$.

We now consider situations where $\text{Tr}[C_{i}\hat{\mathcal{L}}_{1}\rho_{0}]=0$ for all $\rho_{0}$. In this case one has to use perturbation theory at least to order $\epsilon^{2}$ to fix $\rho_{0}$ to order $\epsilon^{0}$. This is always the case for Hamiltonian perturbations where $\hat{\mathcal{L}}_{1}\rho_{0}=-i[\epsilon H_{1},\rho_{0}]$ (this follows from $[C_{i},\rho_{0}]=0$ and the cyclic property of the trace). It corresponds to the well-known fact that transition rates, e.g., calculated from Fermi’s golden rule, are always quadratic in the perturbation. To calculate $\langle\dot{C}_{i}\rangle$ to order $\epsilon^{2}$, $\delta\rho$ has to be expanded up to $\mathcal{O}(\epsilon)$ using Eq. (6) and the relation

[TABLE]

with $X=\hat{\mathcal{L}}_{0}$ and $Y=\hat{\mathcal{L}}_{1}$. All inverses here and in the following are regularized as in Eq. (6). From this equation we find that $\delta\rho=-\hat{\mathcal{L}}_{0}^{-1}\hat{\mathcal{L}}_{1}\rho_{0}+\hat{\mathcal{L}}^{-1}\hat{\mathcal{L}}_{1}\hat{\mathcal{L}}_{0}^{-1}\hat{\mathcal{L}}_{1}\rho_{0}$. A priori, it is not obvious that the first term linear in $\hat{\mathcal{L}}_{1}$ in this expression is large compared to the second one quadratic in $\hat{\mathcal{L}}_{1}$ due to possible singularities in $\hat{\mathcal{L}}^{-1}$ and we will indeed show that this is in general not the case. Nevertheless, using the machinery developed in Sec. III, we can show that for the rate equation under discussion, one can make this approximation, see App. E.2. Therefore we obtain from $\text{Tr}[C_{i}\hat{\mathcal{L}}_{1}\delta\rho]=0$ the condition

[TABLE]

As above, the size of $\epsilon$ drops out of this equation which therefore fixes $\rho_{0}$ to order $\epsilon^{0}$ in situations where (9) is trivially fulfilled. For combined unitary and Markovian perturbations terms from Eq. (9) and Eq. (11) have to be considered simultaneously.

There can be situations where there is no contribution of order $\epsilon^{2}$ to the decay rate of a conservation law such that the left-hand side of both Eq. (9) and Eq. (11) vanishes. This happens, for example, for an integrable Heisenberg chain perturbed by next-nearest neighbor interactions as discussed in Ref. jung06 many years ago. Also in Refs. lenarcic15 ; benenti09 where perturbed $L$-site Ising model and Heisenberg model with strong magnetic field were studied, the boundary Lindblad operators fix the density matrix completely only after terms up to order $L$ or $L-1$ were included, for odd and even system sizes, respectively. In such cases, one has to use Eq. (10) recursively to obtain higher-order corrections.

A non-trivial test of the order of perturbation that fixes $\rho_{0}$ can be obtained by exact (numerical) calculation of the Liouville gap, $\Lambda\equiv\min(|Re\ \lambda|;Re\ \lambda<0),\ \hat{\mathcal{L}}\rho=\lambda\rho$, on a system of smaller size. If scaling $\Lambda\sim\epsilon^{k}$ is obtained then $\rho_{0}$ is determined by the condition of order $\epsilon^{k}$, which is for example $k=1$ for Eq. (9) and $k=2$ for Eq. (11).

II.2 Periodic driving

Above we have discussed stationary states obtained for $t\to\infty$ in systems with time-independent Liouvillians. The same approach can, however, also be used for systems with periodic driving. Here we focus on a case where $\hat{\mathcal{L}}_{0}$ and the associated conservation laws $C_{i}$ are time-independent, while $\hat{\mathcal{L}}_{1}(t)=\hat{\mathcal{L}}_{1}(t+T)$ with period $T$. In this case only minor modifications of the formulation given above are necessary. In the long-time limit the density matrix always has oscillatory (Floquet) components,

[TABLE]

where $\rho^{(-n)}=\rho^{(n)\dagger}$, $\omega=2\pi/T$. The analog of the unique stationary state of Eq. (5) is a state with time-independent Floquet components $\rho^{(n)}$. Due to the limit $\epsilon\to 0$ the zeroth order $\rho_{0}$ contains only time independent $\rho^{(n=0)}$ component. It is useful to organize the Floquet components into a vector $\bm{\rho}=\{\dots\rho^{(-1)},\rho^{(0)},\rho^{(1)},\dots\}$ and promote the Liouvillian into a (static) matrix

[TABLE]

Then the condition for $\rho_{0}$, Eq. (11), has to be reformulated in the above sense: with $\rho_{0}$ and $C_{i}$ interpreted as vectors with non-zero $n=0$ component and $\bm{\hat{\mathcal{L}}_{1}}$, $\bm{\hat{\mathcal{L}}_{0}}$ as matrices. $\bm{\hat{\mathcal{L}}_{1}}$ contains off-diagonal terms due to periodic driving and diagonal terms in case of static perturbations. $\bm{\hat{\mathcal{L}}_{0}}$ contains in addition to diagonal terms due to $\hat{\mathcal{L}}_{0}$ also $in\omega\delta_{nm}$ due to explicit time dependence of $\rho(t)$, $\dot{\rho}=\sum_{n}e^{-in\omega t}(-in\omega\rho^{(n)}+\dot{\rho}^{(n)})$. Stationarity of conserved quantities $C_{i}$, Eq. (11), in case of periodic driving applies to their time averages.

II.3 Projection operators and effective forces

When investigating corrections $\delta\rho$ of the steady state to the zeroth order approximation $\rho_{0}$ (or when investigating the time-dependence of $\rho(t)$), one has to distinguish corrections within the slow subspace with $\hat{\mathcal{L}}_{0}\delta\rho_{\|}=0$ from those perpendicular to this subspace with $\hat{\mathcal{L}}_{0}\delta\rho_{\perp}\neq 0$. Therefore it is useful to introduce a super-projector $\hat{P}_{\rho_{0}}$ which projects on the (density-matrix) space tangential to the GGE manifold at the expansion point $\rho_{0}$. In the following, we will omit the argument $\rho_{0}$, using $\hat{P}=\hat{P}_{\rho_{0}}$ to avoid a cluttering of notations.

Small changes within the manifold around $\rho_{0}$ can be parametrized by $\delta\rho_{\|}=\sum_{i}\delta\lambda_{i}\,\partial\rho_{0}/\partial\lambda_{i}$ where $\delta\lambda_{i}$ are arbitrary infinitesimal changes of the Lagrange parameters. The super-operator $\hat{P}$, projecting on these density matrices, and its complement $\hat{Q}$ are uniquely given by Robertson68

[TABLE]

where $\chi_{ij}=-\text{Tr}[C_{i}\,\partial\rho_{0}/\partial\lambda_{j}]=\langle C_{i}C_{j}\rangle_{{}_{0,\textrm{c}}}$ is a matrix of generalized susceptibilities. We use the notation $\text{Tr}[A\rho]=\langle A\rangle$ and $\text{Tr}[A\rho_{0}]=\langle A\rangle_{{}_{0}}$ for expectation values with respect to $\rho_{0}$. In addition, $\langle AB\rangle_{{}_{0,c}}=\langle AB\rangle_{{}_{0}}-\langle A\rangle_{{}_{0}}\langle B\rangle_{{}_{0}}$ stands for connected correlation function. $\hat{Q}$ and $\hat{P}$ have the property $\hat{Q}^{2}=\hat{Q}$ and $\hat{P}^{2}=\hat{P}$ with $\hat{P}\delta\rho_{\|}=\delta\rho_{\|}$ and are therefore projectors. By construction, $\hat{P}\delta\rho$ is for arbitrary $\delta\rho$ a linear combination of $\partial\rho_{0}/\partial\lambda_{i}$ and therefore an element of the tangential space. Note that $\hat{P}\rho_{0}\neq\rho_{0}$ as $\hat{P}$ is a not a projector on the space of GGE density matrices but instead a projector on the tangential space at $\rho_{0}$. Note that Ref. kawasaki73 uses an alternative projector, which adds an extra term to guarantee $\hat{P}\rho_{0}=\rho_{0}$.

The super-operator $\hat{P}^{\dagger}$, adjoint to $\hat{P}$, has also a direct physical interpretation. The natural scalar product within our approach is $\text{Tr}[A\delta\rho]$, where $A$ is an operator and $\delta\rho$ a density matrix, and therefore $\text{Tr}[(\hat{P}^{\dagger}A)\delta\rho]=\text{Tr}[A(\hat{P}\delta\rho)]$. We obtain

[TABLE]

$P^{\dagger}$ acts on operators and maps each operator onto the space of conserved operators. It gives the part of an operator $A$ which does not decay when the dynamics of the system is described by $\hat{\mathcal{L}}_{0}$ only. The super-operator $\hat{P}^{\dagger}$ naturally shows up when studying the dynamics of systems with conservation laws $C_{i}$. For example, the projection operator used in the memory matrix formalism forster75 ; mori65 ; zwanzig60 ; jung06 ; jung07a can (in this case) be identified with $\hat{P}^{\dagger}$. The operator can also be used to express the Drude weight of conductivities using the seminal results of Mazur mazur69 and Suzuki suzuki71 . The Drude weight $D(T)$ is defined as the prefactor of a $\delta$-function in the optical conductivity, $\text{Re}[\sigma(\omega)]=\pi D(T)\delta(\omega)+\sigma_{\rm reg}(\omega)$. At finite temperatures $T$, the Drude weight is finite in situations where conservation laws $C_{i}$ prohibit the decay of the current. Therefore the Drude weight in a thermal state mazur69 ; suzuki71 ; zotos97 can simply be written in terms of the static cross susceptibility of $J$ and $\hat{P}^{\dagger}J$, $D(T)=\frac{\beta}{L}\langle(\hat{P}^{\dagger}J)J\rangle_{c}$ where $\beta=1/T$ and $L$ is the system size.

Below, we will heavily rely on $\hat{P}$ when deriving perturbation theory for the stationary state. $\hat{P}$ can, however, also be used to define generalized forces allowing to track the changes of the Lagrange parameters during time evolution Robertson66 . This allows to calculate, e.g., the heating of a driven system, and thus to obtain an intuitive picture on the dynamics. Assuming that a state $\rho_{0}(t)$ with time-dependent Lagrange parameters $\lambda_{i}(t)$ describes the system approximately, we can use

[TABLE]

to obtain generalized forces $F_{i}$ governing to leading order the dynamics of the Lagrange parameters. Depending on the studied case it is enough to include only the dominant contribution, e.g., $F_{i}=-\sum_{j}(\chi^{-1})_{ij}\text{Tr}[C_{j}\hat{\mathcal{L}}_{1}\rho_{0}]$. In this paper, we will mainly focus on corrections to the stationary state, but we will discuss dynamics briefly in the context of the Boltzmann equation.

The stationarity condition (9) can be then rewritten as

[TABLE]

Geometrically it means that $\hat{\mathcal{L}}_{1}\rho_{0}$ must be perpendicular ($\perp$) to the slow manifold. Similarly the condition for unitary perturbation, Eq. (11), transcripts into

[TABLE]

In both cases we require that the forces driving the system have effectively to be perpendicular to the GGE manifold.

Note also that in case of periodic driving $\hat{P}$ should be understood as a projector on the slow modes within $n=0$ Floquet sector, Eq. (12). Therefore $\hat{P}\hat{\mathcal{L}}_{1}\rho_{0}=0$ and $\hat{P}\hat{\mathcal{L}}_{1}\hat{\mathcal{L}}_{0}^{-1}\hat{\mathcal{L}}_{1}\rho_{0}=\hat{P}\hat{\mathcal{L}}_{1}\hat{Q}(\hat{Q}\hat{\mathcal{L}}_{0}\hat{Q})^{-1}\hat{Q}\hat{\mathcal{L}}_{1}\rho_{0}$.

II.4 Numerical construction of $\rho_{0}$

In cases where the relevant conservation laws are not known, one can construct for finite systems also the zeroth order result by brute force using exact diagonalization lange17 .

Using the exact eigenstates of $H_{0}$, $H_{0}|n\rangle=E_{n}^{0}|n\rangle$, the set of conservation laws is given by

[TABLE]

To use these conservation laws in Eq. (16) which is written for Hermitian conservation laws, one has to construct the corresponding Hermitian operators $|n\rangle\langle m|+|m\rangle\langle n|$ and $(|n\rangle\langle m|-|m\rangle\langle n|)/i$ for $n\neq m$.

Denoting elements of $\mathcal{Q}$ by $\mathcal{Q}_{i}\in\mathcal{Q}$, one can write $\rho_{0}=\sum\alpha_{i}\mathcal{Q}_{i}$ which fulfills by constuction the condition $\hat{\mathcal{L}}_{0}\rho_{0}=0$. In this case Eq. (9) is a linear equation for the parameters $\alpha_{i}$ which can be solved by finding the kernel of the matrix

[TABLE]

Similarly, the exact solution of Eq. (11) is obtained by finding the eigenvector with eigenvalue 0 of

[TABLE]

Note that the dimension of $\hat{\mathcal{L}}^{\mathcal{Q}}_{nm}$ is much smaller than the dimension of the full Liouvillian super-operators. For a spin-1/2 XXZ Heisenberg chain with $L$ sites in the presence of an external magnetic field, we found in Ref. lange17 that the dimension of $\mathcal{Q}$ is approximately $2\cdot 2^{L}$ (the factor $2$ arises from degeneracies in the spectrum), to be compared with $(2^{L})^{2}=4^{L}$, the dimension of density matrices on which the Liouville super-operators act.

A tricky, unresolved question is under what condition the construction above converges in the thermodynamic limit. For example, there could in principle be situations where a conservation law exists only in the thermodynamics limit but not in the finite size system. A related question is whether the limits $\epsilon\to 0$ and $L\to\infty$ commute or not.

When using the set $\mathcal{Q}$ of conserved quantities, the projector $\hat{P}$ has to be replaced by $\tilde{P}$,

[TABLE]

projecting on the density operator subspace spanned by $|n\rangle\langle m|,\ E_{n}^{0}=E_{m}^{0}$. This type of projector has been used also in other perturbative studies of Liouville dynamics benatti11 ; znidaric15 implementing degenerate perturbation theory for super-operators in case of weak openness.

III Perturbation theory

Projectors $\hat{P},\hat{Q}$ are necessary to develop the perturbation theory of corrections to $\rho_{0}$, given to all orders by $\delta\rho$ in Eq. (6). Below we present a controlled expansion of the inverse $\hat{\mathcal{L}}^{-1}$, which removes possible singularities and separates different orders $\mathcal{O}(\epsilon^{n})$ within the steady state. With $\parallel$ and $\perp$ we denote the density matrix subspaces, which are in the image of $\hat{P}$ and $\hat{Q}$, respectively, $\delta\rho=\hat{P}\delta\rho+\hat{Q}\delta\rho=\delta\rho_{\|}+\delta\rho_{\perp}$. Note that only the tangential component is relevant for the expectation values of conserved operators $C_{i}$,

[TABLE]

as $\hat{Q}^{\dagger}C_{i}=0$ for conserved quantities. $\delta\rho_{\perp}$ does, however, affect expectation values of other operators.

III.1 Markovian perturbation

First we consider Liouvillians where the perturbations break all relevant conservation laws already to linear order in $\hat{\mathcal{L}}_{1}$ such that $\hat{P}\hat{\mathcal{L}}_{1}\hat{P}$ has no zero modes. Otherwise the procedure explained in Sec. III.3 is needed.

Our goal is to expand $\delta\rho=-\hat{\mathcal{L}}^{-1}\hat{\mathcal{L}}_{1}\rho_{0}$, Eq. (6), using that we determined $\rho_{0}$ from the condition $\hat{P}\hat{\mathcal{L}}_{1}\rho_{0}=0$, Eq. (20). With this definition, we obtain

[TABLE]

where the Liouville inverse has to be interpreted using regulator as in Eq. (6), see also App. A. Due to its conservation laws, $\hat{\mathcal{L}}_{0}^{-1}$, is singular. To avoid these singularities, we expand $\hat{\mathcal{L}}^{-1}$ around $\hat{\mathcal{L}}_{0}+\hat{P}\hat{\mathcal{L}}_{1}\hat{P}$ using $\hat{\mathcal{L}}_{1}=(\hat{P}+\hat{Q})\hat{\mathcal{L}}_{1}(\hat{P}+\hat{Q})$ and that $\hat{P}\hat{\mathcal{L}}_{1}\hat{P}$ is invertible in P subspace,

[TABLE]

For power counting in $\epsilon$, it is important to realize that

[TABLE]

where we have used $\hat{P}\hat{\mathcal{L}}_{0}=\hat{\mathcal{L}}_{0}\hat{P}=0$. This reflects that the dynamics in the subspace of approximately conserved quantities is slow, while it is fast in the perpendicular space. Note that $(\hat{P}\hat{\mathcal{L}}_{1}\hat{P})^{-1}$ and $(\hat{Q}\hat{\mathcal{L}}_{0}\hat{Q})^{-1}$ in Eq. (29) should be interpreted as inverses within the $P$ and $Q$ subspaces, respectively, where $\hat{\mathcal{L}}_{0}$ and $\hat{\mathcal{L}}_{1}$ are invertible. One first performs projection and then the inversion within the subspace that dynamics was projected on. Moreover, Eq. (29) gives an alternative derivation of Eq. (20): the steady state $\rho_{0}$ has to fulfill $\hat{P}\hat{\mathcal{L}}_{1}\rho_{0}=0$, otherwise corrections $\delta\rho\sim\hat{\mathcal{L}}^{-1}\hat{P}\hat{\mathcal{L}}_{1}\rho_{0}\sim 1$ which contradicts our perturbative approach in the limit $\epsilon\to 0$, see App. E.2.

The combination of Eq. (27), (28) and (29) allows to obtain a straightforward expansion of $\delta\rho$ in powers of $\epsilon$. The steady state density matrix has a distinct structure in the tangential ($\parallel$) and the perpendicular ($\perp$) subspace in all orders, therefore one has to consider the contributions from the two subspaces separately. To linear order in $\epsilon$ we obtain

[TABLE]

For brevity we use notation $\hat{Q}(\hat{Q}\hat{\mathcal{L}}_{0}\hat{Q})^{-1}\hat{Q}=\hat{Q}\hat{\mathcal{L}}_{0}^{-1}\hat{Q}$, where $\hat{\mathcal{L}}_{0}^{-1}$ should be understood as an inverse within $Q$ subspace only. We would like to stress again, that $\rho_{0}$ already contains effects of order $\epsilon^{0}$. But after $\rho_{0}$ has been correctly chosen, $\delta\rho$ does indeed vanishes for $\epsilon\to 0$, $\delta\rho\sim\epsilon$. In the $\|$ space, one has, however, to expand to second order in $\hat{\mathcal{L}}_{1}$ to obtain the correction of order $\epsilon$ due to the presence of the $(\hat{P}\hat{\mathcal{L}}_{1}\hat{P})^{-1}$ term. This contributes as $1/\epsilon$ factor, resulting in a term proportional to $\epsilon^{2}/\epsilon=\epsilon$. In the perpendicular space no such issue arises. The fact that perturbation theory for steady states has a different structure compared to perturbation theory in thermal equilibrium is well known, see, e.g. Ref. rosch03 , and can simply be understood from the fact that rate equations determining the steady state have to be readjusted in the presence of perturbations.

Note that the change $\delta\rho_{\|}$ can in principle be absorbed in a redefinition of the Lagrange parameters, while $\delta\rho_{\perp}$ describes contributions which cannot be captured by a GGE ensemble.

For concrete calculations one needs to compute

[TABLE]

where $M^{-1}$ is the inverse of matrix $M$ with components $M_{pr}=-\text{Tr}[C_{p}\,\hat{\mathcal{L}}_{1}(\partial\rho_{0}/\partial\lambda_{r})]$, see App. B. Then, for example, the change of the expectation value of a conserved operator $C_{i}$ is given linearly in $\epsilon$ by, $\delta\langle C_{i}\rangle\approx\langle C_{i}\rangle_{1}$

[TABLE]

Higher order corrections can be obtained in a straightforward way from the Taylor expansion, Eq. (28), in combination with the rules of Eq. (29) which show how the inverse $(\hat{\mathcal{L}}_{0}+\hat{P}\hat{\mathcal{L}}_{1}\hat{P})^{-1}$ has to be evaluated. Eq. (33) and its higher order generalization can be rewritten in terms of standard correlation functions as we discuss in App. D where concrete formulas for $\langle C_{i}\rangle_{1}$ and $\langle C_{i}\rangle_{2}\sim\mathcal{O}(\epsilon^{2})$ are given.

In Fig. 1 we give a graphical representation of the relevant terms in the tangential ($\parallel$) and the orthogonal ($\perp$) subspace. As written on the left of the figure, the two types of inverse Liouville super-operators in $P$ and $Q$ space are shown as white and black circles, respectively. There are three types of links, $\hat{Q}\hat{\mathcal{L}}_{1}\hat{Q}$, $\hat{P}\hat{\mathcal{L}}_{1}\hat{Q}$ and $\hat{Q}\hat{\mathcal{L}}_{1}\hat{P}$, connecting black and black, white and black, and black and white circles, respectively. This simply follows from the rule that $\hat{P}\hat{Q}=0=\hat{Q}\hat{P}$.

Starting from $\hat{Q}\hat{\mathcal{L}}_{1}\rho_{0}$ (small square with thin line) to the right, one attaches circles and lines in all allowed combinations. The color of the last circle on the left, determines whether the density matrix operator is in $P$ or $Q$ space. The power in $\epsilon$ is determined by the number of lines minus the number of white circles as $(\hat{P}\hat{\mathcal{L}}_{1}\hat{P})^{-1}\sim 1/\epsilon$. The sign is simply given by $(-1)^{N_{L}}$ where $N_{L}$ is the number of lines.

We have checked that the formulas discussed above can alternatively be derived starting from projection-operator based time-dependent perturbation theories in the version discussed by Breuer and Petruccione breuer02 ; Koide02 .

As a remark we should point out that there is no need for any additional normalization since $\text{Tr}\,\delta\rho=0$ which is guaranteed by the regularization of the inverse operators, see Appendix A.

III.2 Missing conservation laws

Above, we have assumed that all relevant conservation laws of the unperturbed system are known and have been included in the construction of the GGE (or that $\rho_{0}$ has been constructed numerically, see Sec. II.4). It may, however, happen that either not all conservation laws are known or that it is just technically impossible to include them all. This might, for example, be the case in integrable systems where an infinite number of conservation laws exists.

In this case one may try to approximate system by a truncated GGE,

[TABLE]

considering only a finite subset of conserved operators. We denote by $\hat{P}_{t}$ the projector on the space tangential to the truncated GGE, and by $\chi^{[tt]}$ the matrix of susceptibilities, all defined as above.

A selection criterion for the truncated space could be, for example, to consider only the most local conservation laws, containing less than a certain number of derivatives in a continuum model or having a support of less than a certain number of sites for a lattice model. A truncated GGE was, for example, used in Ref. lange17 . In the following we will consider two questions: (i) How can one compute perturbatively the effects of the ’other’ conservation laws, and (ii) how does the perturbation theory developed above signal the presence of missing conservation laws.

In the following we denote the missing conservation laws not included in the truncated GGE by $C_{i}^{[m]}$. Similarly, we define susceptibility matrices $\chi^{[mm]}_{ij}=\langle C_{i}^{[m]}C_{j}^{[m]}\rangle_{{}_{0,\textrm{c}}}$ and $\chi^{[mt]}_{ij}=\langle C_{i}^{[m]}C_{j}^{[t]}\rangle_{{}_{0,\textrm{c}}}$, where $\langle\cdot\rangle_{{}_{0}}$ stands for expectation values with respect to $\rho_{0}^{[t]}$. Without loss of generality, we assume that they are orthogonal to the $C_{i}^{[t]}$, $\chi^{[mt]}_{ij}=0$ for all $i,j$ (if this is not the case, one can simply replace them by ${\hat{Q}_{t}}^{\dagger}C_{i}^{[m]}$ where $\hat{Q}_{t}=\hat{1}-\hat{P}_{t}$). We define the projector on the tangential space corresponding to the missing conservation laws by

[TABLE]

with (\partial\rho_{0}/\partial\lambda_{i}^{[m]})\big{|}_{\rho_{0}^{[t]}}=-(C^{[m]}_{i}-\langle C^{[m]}_{i}\rangle_{{}_{0}})\rho_{0}. The sum $\hat{P}=\hat{P}_{t}+\hat{P}_{m}$ is then the projector on the tangential space spanned by all conservation laws.

Assuming that the missing conservation laws give only a small correction to the Lagrange parameters $\lambda_{i}$, we can Taylor expand in $\delta\lambda_{k}=\lambda_{k}-\lambda_{k}^{[t]}$ using

[TABLE]

where the $k$ sum includes both the $C_{i}^{[t]}$ and $C_{i}^{[m]}$. From Eq. (9), we obtain directly a matrix equation $A\bm{\delta\lambda}=\bm{a}$ for $\delta\lambda_{i}$ solved by $\bm{\delta\lambda}=A^{-1}\bm{a}$ which is written in components as

[TABLE]

Note that $A^{[mm]}$ is equivalent to matrix $M$ in Eq. (32).

From the change of the Lagrange parameters, one can also directly calculate the change of observables. We are mainly interested in the change of $\langle\delta C_{i}^{[t]}\rangle_{{}_{0}}$

[TABLE]

In the last line we did an extra approximation, assuming that the relevant matrix elements between included and missing conservation laws, $C^{[t]}_{i}$ and $C^{[m]}_{j}$, are small, allowing to expand $A^{-1}$ in $A^{[tm]}$. We conclude that the effect of the missing conservation laws is small if either $\langle\dot{C}_{i}^{[m]}\rangle_{{}_{0}}$, the changes of missing conservation laws, are small and/or if the dynamical coupling $A^{[tm]}_{ij}=\langle\dot{C}^{[t]}_{i}\bar{C}^{[m]}_{j}\rangle_{{}_{0,\textrm{c}}}$ is small. The latter susceptibility describes how a change of the Lagrange parameter $\lambda^{[m]}_{j}$ induces a finite $\langle\dot{C}^{[t]}_{i}\rangle_{{}_{0}}$.

Fully equivalent to Eq. (45), one can write everything using projection super-operators $\hat{P}_{t},\hat{P}_{m}$

[TABLE]

Above, we have calculated effects which occur in cases when the problem was not solved accurately to order $\epsilon^{0}$ because not all conservation laws of $H_{0}$ had been considered. If a theory is not treated correctly to order $\epsilon^{0}$, any perturbative treatment in $\epsilon$ should signal this by the presence of divergencies. In the following, we will show that this is the case within our approach and we will identify directly the origin of these divergencies.

Technically, the divergence in the perturbative formula (33) arises from the inverse super-operator $\hat{Q}_{t}(\hat{Q}_{t}\hat{\mathcal{L}}_{0}\hat{Q}_{t})^{-1}\hat{Q}_{t}$ when only a truncated set of conservation laws is used. As $\hat{P}_{m}$ is the projection operator on the missing conservation laws, we can calculate the divergent contribution from

[TABLE]

where $\eta$ is the regulator used to define $\hat{\mathcal{L}}_{0}^{-1}$, see appendix A. We note that exactly the same divergences lead to the occurrence of infinities in the current response of, e.g., integrable systems. As already discussed in Sec. II.3, the latter are described by Drude weights mazur69 ; suzuki71 ; zotos97 obtained from $\langle(\hat{P}^{\dagger}J)J\rangle_{{}_{c}}$.

The divergent contribution to $\langle\delta C^{[t]}_{i}\rangle$ within $O(\epsilon)$ perturbation theory around the truncated GGE is obtained from Eqs. (33) and (47) as

[TABLE]

This contribution is of order $\epsilon/\eta$. Comparing Eq. (48) and the second line of Eq. (46) one observes that the two terms are identical if one replaces $(-\eta)$ by $\hat{P}_{m}\hat{\mathcal{L}}_{1}\hat{P}_{m}$ which is linear in $\epsilon$. Divergencies in perturbation theory can therefore be used to detect missing conservation laws. If, however, the weight of the divergent contribution is small, one can expect that adding the missing conservation laws will have little effect.

III.3 Unitary driving,

$\hat{P}\hat{\mathcal{L}}_{1}\hat{P}=0$

In the discussion given above, we assumed that $\hat{P}\hat{\mathcal{L}}_{1}\hat{P}$ is finite and invertible within $P$ subspace. For an important class of perturbations, arising from a (Floquet) Hamiltonian $H_{1}$, this is not the case as $\text{Tr}(C_{i}[H_{1},\partial\rho_{0}/\partial\lambda_{j}])=0$ due to the cyclic property of the trace and one therefore finds $\hat{P}\hat{\mathcal{L}}_{1}\hat{P}=0$. In this case, the rates with which the $C_{i}$ change are second order in $\epsilon$, as is well-known from Fermi’s golden rule.

For $\hat{P}\hat{\mathcal{L}}_{1}\hat{P}=0$ the exact inverse of the Liouvillian in the $\hat{P}$ sector is given by $\hat{P}\hat{\mathcal{L}}^{-1}\hat{P}=-\hat{P}(\hat{P}\hat{\mathcal{L}}_{1}\hat{Q}(\hat{Q}\hat{\mathcal{L}}\hat{Q})^{-1}\hat{Q}\hat{\mathcal{L}}_{1}\hat{P})^{-1}\hat{P}$, see appendix E.1. In the limit of small $\epsilon$, we therefore find

[TABLE]

where

[TABLE]

$\hat{P}\hat{\mathcal{L}}_{2}\hat{P}\propto\epsilon^{2}$ is an effective Lindblad super-operator acting in the $P$ space. With this notation the condition for $\rho_{0}$, Eq. (11), takes the form $\overline{\text{Tr}[C_{i}\hat{\mathcal{L}}_{2}\rho_{0}]}=0$, or, equivalently, Eq. (21), is written as

[TABLE]

where one should keep in mind that $\hat{P}$ project onto non-oscillatory density matrices only. Perturbation theory can now be derived in a straightforward way by a Taylor expansion of

[TABLE]

in the second term, with the inverse applied as in equation (6) for $\delta\rho$.

One can, again, derive diagrammatic rules to calculate corrections to order $\epsilon^{n}$, see Fig. 2. Compared to the previous case, shown in Fig. 1, there are only two changes. First, $(\hat{P}\hat{\mathcal{L}}_{1}\hat{P})^{-1}\sim\mathcal{O}(\epsilon^{-1})$ (open circle) is replaced by $(\hat{P}\hat{\mathcal{L}}_{2}\hat{P})^{-1}\sim\mathcal{O}(\epsilon^{-2})$ and, second, two sets of diagrams do not contribute any more as they either cancel due to the $-\hat{P}\hat{\mathcal{L}}_{2}\hat{P}$ in Eq. (52) or are set to zero by the condition (51) for $\rho_{0}$. Explicitly, there are no diagrams with neighboring open circles, with the combination open-filled-open, and finally also the combination open-fillled-$\rho_{0}$ (small square). The cancellation of these diagrams is also a necessary condition for the series expansion in $\epsilon$ to be valid. In Fig. 2, we show all remaining diagrams up to $\mathcal{O}(\epsilon^{3})$ and describe in the figure caption the corresponding diagramatic rules.

We would like to finish the section by pointing out that there can be more complex situations jung06 ; lenarcic15 where perturbation theory to order $\epsilon^{2}$ does not fix $\rho_{0}$ and where accordingly $\hat{P}\hat{\mathcal{L}}_{2}\hat{P}$ vanishes, see discussion at the end Sec. II.1. Or even more general, there can be cases where Hamiltonian and Lindblad perturbations occur on equal footing or where some approximate conservation laws change by processes to order $\epsilon$ and others by order $\epsilon^{2}$. We believe that in all these cases, one can generalize the approach described above: one defines projectors on subspaces governed by similar time-scales (same power of $\epsilon$), identifies the leading-order dynamics ($\hat{P}\hat{\mathcal{L}}_{2}\hat{P}$ in the case discussed above) in each subsector perturbatively, uses this to fix $\rho_{0}$ in analogy to Eq. (51), and performs then a Taylor expansion using the analog of Eq. (52) as a starting point.

IV Example: Interacting Fermions with particle-gain and loss

In the following we will give two examples where the perturbation theory developed above applies. We consider interacting fermions driven out of equilibrium by a weak coupling to the environment. As perturbation we choose processes which lead to a particle loss with rate $\epsilon l_{n}$ and a gain of particles with rate $\epsilon g_{n}$. In quantum optics experiments, the quasiparticles could be excitations of atoms (or cavities) and gain and loss is realized by emission and absorption of light. The same situation arises in experiments on exciton condensation, magnon condensation or in photon BECs butov01 ; kasprzak06 ; balili07 ; klaers10 ; klaers12 ; schmitt14 ; demokritov06 ; chumak09 . Also in experiments with ultracold atoms, loss processes arise when atoms absorb photons from external laser beams, kicking them out of their trap. However, most cold-atom experiments do not include processes where the lost atoms are replenished and therefore a true steady state, the main focus of our study, cannot be reached.

We first study in Sec. IV.1 a small finite-size system where we can compare directly perturbation theory and exact solution. Here it is important to realize that the applicability of perturbation theory for small systems (in the limit where perturbations are small compared to the level spacing) does not guarantee the validity of perturbation theory for systems close to the thermodynamic limit (perturbations large compared to the level spacing but small compared to internal equilibration rates). We therefore study in Sec. IV.2 the thermodynamic limit by considering a regime where the Boltzmann equation can be applied. A numerical investigation of the thermodynamic limit in a quantum approach is beyond the scope of the present paper and is left for future studies inpreparation .

IV.1 Lindblad dynamics in a small system

We consider a model where the Hamiltonian dynamics of the unperturbed system is given by

[TABLE]

with $e_{n}=n/L$. Gain and loss processes are described by the Lindblad operators $L^{g}_{n}=c^{\dagger}_{n}$ and $L^{l}_{n}=c_{n}$ such that

[TABLE]

with

[TABLE]

The loss rate $\epsilon l_{n}$ depends on energy. This mimics a situation where high-energy states have a higher probability to evaporate than low-energy states, thus implementing a cooling mechanism. Note that $\epsilon$, the common prefactor of $\hat{D}_{i}$, controls the overall strength of both the heating and cooling terms.

The conservation laws of the small system are just projectors on the eigenstates of $H_{0}$, $|n\rangle\langle n|$ with $H_{0}|n\rangle=E_{n}^{0}|n\rangle$. Using the formulas of Secs. II,III we can easily solve for the 0th order steady state and determine the perturbations around it. In Fig. 3 we show the expectation values of particle and energy density as a function of perturbation strength $\epsilon$, calculated from the exact steady state density matrix (solid line) or using our perturbation theory up to $k$-th order in $\epsilon$ on systems size $L=4$, $k=2,4,10$, and interaction strength $U=0.3$.

The figures show that for this small system with finite level spacing the perturbation theory works for small $\epsilon$. For the chosen model we find both from the exact result and from the perturbation theory that the variations of particle number and energy are tiny. The perturbation theory breaks down for rather small $\epsilon$ of the order of the level spacing. This is a phenomenon well known from standard perturbation theory where perturbative corrections are inversely proportional to the level spacing. This is, however, dfferent when correlations functions are evaluated in the thermodynamics limit by taking first the limit $L\to\infty$ and then the limit $\eta\to 0$, where $\eta$ is the regulator defined in Eq. (6), see also App. D. While the perturbative correction linear in $\epsilon$ vanishes for the finite system ($\eta$ smaller than level spacing), it is finite in the thermodynamics limit ($\eta$ larger than level spacing).

IV.2 Boltzmann dynamics

Now we consider a similar example in the thermodynamic limit where the unperturbed system $H_{0}$ has only two local conservation laws, energy and particle number. Interactions lead to an equilibration of all other conservation laws. We assume that the system can be described by weakly interacting quasiparticles with energies $e_{n}$, such that the collision term of a Boltzmann equation captures their equilibration dynamics. As above we consider fermions with a constant density of states, which we discretize using $L$ single-particle states with energies $e_{n}$ equally spaced between [math] and $1$ with $e_{n}=n/L$, $n=1,\dots,L$. The Boltzmann equation takes the form,

[TABLE]

where $f_{e_{n}}$ is the occupation function as function of the energy $e_{n}$ and $\bar{f}_{e}=(1-f_{e})$. The Kronecker-$\delta$ guarantees energy conservation in each collision process of the discretized model. Note that we consider a model where momentum conservation does not play a role. For simplicity, we also set all transition rates due to collisions to unity.

We keep the same type of perturbations as in the previous example, i.e. particle gain and loss. Instead of Lindblad operators these are now encoded directly in Boltzmann equation through $D[f]_{e}$, Eq. (57). Rates for particle gain and loss are the same as in Eq. (55), i.e., $g_{e}=1/4,l_{e}=e$. Note that $\epsilon$, the prefactor of $D[f]_{e}$, controls the overall strength of both the heating and cooling terms. Numerical calculations are performed for $L=41$.

The perturbation theory derived in the main text can be with some straightforward modifications also applied to the open Boltzmann dynamics. Before we proceed we collect in Table 1 all analogies with the formulation for the Liouvillian dynamics, developed in the main part of the article. Note that some objects are defined later within this section.

The collision integral $M$ preserves the total energy $E=\sum_{n}e_{n}f_{e_{n}}$ and the particle number $N=\sum_{n}f_{e_{n}}$ and these are the only conserved quantities in the absence of loss/gain term. In the following we expand the level occupations $f_{e}$ in orders of $\epsilon$,

[TABLE]

In correspondence with $\hat{\mathcal{L}}_{0}\rho_{0}=0$ now any Fermi-Dirac distribution $f_{e}^{0}$ satisfies $M[f^{0}]=0$,

[TABLE]

and only the perturbation $D$ fixes the parameters $\beta$ and $\mu$, where $\beta$ is the Lagrange parameter of the energy and $-\beta\mu$ is the Lagrange parameter of the particle number.

Zeroth order— The parameters $\beta$ and $\mu$ of the steady-state in the $\epsilon\to 0$ limit, $f^{0}_{e}(\beta,\mu)$, are determined by the stationarity of conserved quantities (total energy and particle number) from a set of coupled equations corresponding to Eq. (9),

[TABLE]

For our model we find $\beta=2.328$ and $\mu=0.288$.

Projection operator — The two slow modes $q_{i}$, corresponding to $\partial\rho_{0}/\partial\lambda_{i}$ in the main text, are defined as

[TABLE]

where we identified $\lambda_{2}=-\beta\mu$. Denoting the conserved quantities by

[TABLE]

the analog of the projector super-operator defined by Eq. (16) is simply given by

[TABLE]

With these notations, the steady state condition, Eq. (IV.2), is the analog of Eq. (9)

[TABLE]

or, equivalently,

[TABLE]

Relaxation towards the steady state — While the main focus of this paper is the computation of the steady state, we briefly discuss the relaxation towards the steady state, which can easily be computed exactly by solving the Boltzmann equation.

Similarly as GGEs are expected to be a fairly good description shortly after the system has prethermalized berges04 ; Aarts2000 ; Gring2012 ; Langen2015 ; Langen2016 ; moeckel08 ; Kollath2007 ; kollar11 ; Marino2012 ; Mitra2013 ; Marcuzzi2013 ; bertini15 ; Demler2015 ; Buchhold2014 ; Marcuzzi2016 ; Barnett2011 ; NessiIucci2014 ; Nowak2014 ; canovi16 ; golez16 ; Chiocchetta2017 to the GGE manifold, also here the Fermi-Dirac distribution describes approximately not only the steady state but also the approach to it. After a few collisions level occupations can be approximately described by a Fermi distribution with time-dependent parameters $\{\beta(t),\mu(t)\}$, whose evolution is determined by the forces introduced in Eq. (19).

In this case force field, see Eq. (19), has only two components,

[TABLE]

and can be used to propagate parameters $\beta$ and $\mu$,

[TABLE]

as soon as our system is approximately described by a Fermi-Dirac distribution. These equations are only valid to leading order in $\epsilon$.

Fig. 4 shows the time-evolution obtained from a numerical solution of the Boltzmann equation for $\epsilon=0.01$. We initialize our system with the non-equilibrium state $f_{e_{n}}=1$ for states with even and $f_{e_{n}}=0$ for states with odd $n$. Collisions lead to a rapid relaxation of this non-equilibrium state to a thermal state. We find that within our model this time scale is $\tau_{0}\sim 5$, independent of $\epsilon$ for small $\epsilon$. As the collision processes conserve energy and particle number, the approximate thermal state is determined by the initial values of energy and particle number, which leads to an initial relaxation towards a state with $\beta=0$ and $\beta\mu=0$ (in the thermodynamic limit). The subsequent dynamics is driven by the processes violating particle number and energy conservation (the gain and loss terms) and occurs on a time scale set by $1/\epsilon$. We find that these processes are quantitatively described by Eq. (65), as can be seen by comparing the dots and the solid line in the lower panel of Fig. 4.

Fig. 5 shows the force fields calculated according to Eq. (64) and the resulting dynamics in the space of Lagrange parameters.

Perturbation theory for the steady state — As a next step, we compare the exact steady-state solution of the Boltzmann equation for finite $\epsilon$ to the result of perturbation theory to linear order in $\epsilon$. We expand around the Fermi functions $f^{0}_{e}(\beta,\mu)$ with parameters $\beta$ and $\mu$ obtained from the ‘zeroth order’ described above. To linear order in $\epsilon$, we obtain using $f_{e}=f^{0}_{e}+\epsilon f^{1}_{e}$ to $\mathcal{O}(\epsilon)$

[TABLE]

which has precisely the same form as Eq. (6) in the Liouville case. For the discrete case matrices $M^{(0)},D^{(1)}$ are defined using a straightforward Taylor expansion

[TABLE]

Note that $D^{(1)}$ is a diagonal matrix in our example. To obtain $\bm{f}^{1}=\{f^{1}_{e_{1}},f^{1}_{e_{2}},\dots,f^{1}_{e_{L}}\}$ corrections we can use the perturbation theory developed in Sec. III.1. As before we need to treat the ($\parallel$) and ($\perp$) subspaces separately. Following Eqs. (30,31),

[TABLE]

where $P,Q$ are written as matrices after evaluating Eq. (63) for discretized energies and $\hat{Q}=\hat{1}-\hat{P}$. Fig. 6 shows steady state level occupation $f_{e_{n}}$ as a function of perturbation strength $\epsilon$. In the limit $\epsilon\to 0$ Fermi functions obtained from condition Eq. (IV.2) give the exact steady state occupation distributions. For finite $\epsilon\lesssim 0.05$ the corrections linear in $\epsilon$ describe well the exact result. In Fig. 7 we show similar results for the total energy and particle number of the system, which for the chosen model change only slightly.

It is also straightforward to calculate higher order corrections using perturbation theory applied to the Boltzmann equation. However, here one has to take into account that the Boltzmann equation is a non-linear equation in the occupation functions $f_{e}$ while the Liouville equation is a linear equation in the density matrix $\rho$. Therefore the equations used to obtain higher order corrections differ in the two cases but the overall structure of perturbation theory (separation into parallel and perpendicular sector) remains the same.

V Conclusions

The main observation, which is the basis of our study, is that many driven systems can be approximately but efficiently described by generalized Gibbs ensembles built from approximately conserved quantities. Such approximate conservation laws are important for many different systems. For example, in almost all solids, the coupling of electrons and phonons is weak due to the large mismatch of ionic and electronic masses. Therefore the difference of the electron and phonon Hamiltonians, $H_{e}-H_{\rm ph}$, is an approximate conservation law widely used in two-temperature models anisimov74 ; armstrong81 ; eesley86 ; allen87 ; singh10 which efficiently describe non-equilibrium states of solids. A more exotic example is the description of spin-chain materials using as the approximate conservation laws those of the underlying integrable Heisenberg model lange17 . For such situations, we developed in this paper a perturbation theory for the steady state by perturbing around generalized Gibbs ensembles. Here it was essential to treat perturbations perpendicular and parallel to the manifold of generalized Gibbs states separately.

To show the validity of our concepts we studied a small finite size quantum system with Lindblad dynamics. For the future inpreparation we plan to test numerically our approach also for larger quantum systems close to the thermodynamic limit. In this paper we considered as an example for a large system only a model described by the Boltzmann equation. From a more general point of view, one can consider the Boltzmann equation also as a rate equation for approximately conserved quantities, the quasi-particle occupations $c^{\dagger}_{k}c_{k}$, thus fixing the parameters of the appropriate generalized Gibbs ensemble $e^{-\sum_{k}g_{k}c^{\dagger}_{k}c_{k}}$ (as in our paper, we assume a translationally invariant system here). Therefore it would be interesting to use our formalism to calculate corrections to the non-equilibrium steady state not captured by the Boltzmann approach. A comparison of the two formalisms and a systematic calculation of beyond-Boltzmann corrections for weakly interacting quantum systems is left for further studies. In the classical context such a formalism has been developed by Gurarie gurarie95 . It would be in general interesting to test the validity of our approach in purely classical systems, for example, those described by the Fokker-Planck equation.

An obvious question is whether our perturbation theory in $\epsilon$, where $\epsilon$ is the prefactor of driving terms breaking the conservation laws, has a finite radius of convergence and/or whether all expansion coefficients are finite. Similarly to many other (highly successful) perturbative expansions in physics, the radius of convergence is formally expected to be zero. This is obvious in the case when the perturbation arises from a Lindblad operator, where negative $\epsilon$ correspond to negative friction and therefore an unstable situation for arbitrarily small $\epsilon$, which implies a vanishing radius of convergence. The expansion in $\epsilon$ should therefore be viewed as an asymptotic expansion similar to, e.g., the expansion in the fine-structure constant $\alpha$ in quantum-electrodynamics. A more difficult question is whether one can expect that all expansion coefficient in powers of $\epsilon$ are finite or not. This is an open question which requires further studies. There is, however, a possible mechanism which can induce non-analytic corrections: in the presence of exact symmetries, hydrodynamic modes exist which enforce long-ranged correlations and, for example, the presence of long-time tails after quenches lux14 ; bohrdt17 ; leviatan17 . These hydrodynamic modes obtain masses proportional to $\epsilon$ or $\epsilon^{2}$ when the relevant conservation laws are weakly violated. A perturbative expansion in these masses can lead to non-analytic corrections even when perturbing around a GGE.

For the future it will be interesting to generalize our perturbative expansion for the steady state to dynamical questions. For example, we have already shown for a simple example based on the Boltzmann equation that time-dependent generalized Gibbs ensembles are a good starting point to investigate the approach of the steady state. It will be interesting to develop a perturbative expansion for such situations which will also allow to describe the dynamical response of driven systems and/or situations where instead of a unique steady state a limit cycle (time crystal) is realized.

Acknowledgements.

We acknowledge useful discussions with M. Genske, J. Lux, M. Žnidarič and financial support of the German Science Foundation under CRC 1238 (project C04).

Appendix A Regularization of inverse

Since $\hat{\mathcal{L}}$ has a zero mode, $\hat{\mathcal{L}}^{-1}$ is singular therefore formally one has to regularize it as in Eq. (6), which at the same time also avoids solution $\delta\rho=-\rho_{0}$ that satisfies $\hat{\mathcal{L}}(\rho_{0}+\delta\rho)$ trivially. If we write $\rho_{0}=\sum_{\alpha}c_{\alpha}\rho^{\alpha}$ in terms of $\hat{\mathcal{L}}$ (right) eigenstates $\rho^{\alpha}$, $\hat{\mathcal{L}}\rho^{\alpha}=\lambda_{\alpha}\rho^{\alpha}$, then the regularization

[TABLE]

gives the correct result, $\rho_{0}+\delta\rho=\sum_{\alpha,\lambda_{\alpha}=0}\rho^{\alpha}$. The sign $\eta>0$ is obtained from the property that $\text{Re}\,\lambda_{\alpha}\leq 0$ guaranteeing the absence of exponentially growing solutions.

The regularization also guarantees that $\text{Tr}\,\delta\rho=0$. From $\text{Tr}\,\dot{\rho}=\text{Tr}[\hat{\mathcal{L}}\rho]=0$, it follows that $\lambda_{\alpha}\text{Tr}\,\rho^{\alpha}=0$ and therefore $\text{Tr}\,\rho^{\alpha}=0$ for all $\alpha$ with $\lambda_{\alpha}\neq 0$. Using (69), we therefore find $\text{Tr}\,\delta\rho=0$.

Appendix B Identities for super-operators

Below we check the properties of super-projector defined in Eq. (16).

Projection property

[TABLE]

*Inverse within $P$ subspace

*Since $(\hat{P}\hat{\mathcal{L}}_{1}\hat{P})^{-1}$ is essential for building up the perturbation theory we confirm the validity of Eq. (32)

[TABLE]

with $M_{pr}=-\text{Tr}[C_{p}\hat{\mathcal{L}}_{1}\frac{\partial\rho_{0}}{\partial\lambda_{r}}]=\langle\dot{C}_{p}C_{r}\rangle_{{}_{0,c}}$ by calculating

[TABLE]

Appendix C Properties of projectors and inversion of projected Lindblad operators

In this appendix we discuss the relation of the projectors $\hat{P}$ and $\tilde{P}$ and describe how projected Lindblad operators can be inverted.

We first note that any projected density matrix is traceless

[TABLE]

which follows from $\text{Tr}[\rho_{0}]=1$ and therefore $\partial\text{Tr}[\rho_{0}]/\partial\lambda_{i}=0$.

If we consider all conservation laws from the set $\mathcal{Q}=\{|n\rangle\langle m|\ {\rm with}\ E_{n}^{0}=E_{m}^{0}\}$, $H_{0}|n\rangle=E_{n}^{0}|n\rangle$, then we obtain

[TABLE]

where $\tilde{P}X$ is part of $X$ that can be written in terms of elements of $\mathcal{Q}$ (for non-degenerate case simply the diagonal part). Therefore $\hat{P}$ projects on conservation laws and subtracts the trace.

Assuming that $\hat{P}\hat{\mathcal{L}}_{1}\rho_{0}$, Eq. (20), has a unique solution with trace 1, then one can easily show that $\hat{P}\hat{\mathcal{L}}_{1}\hat{P}$ is invertible in $\hat{P}$ space. Assuming that $\hat{P}\hat{\mathcal{L}}_{1}\hat{P}\rho_{1}=0$, it follows immediately that $\hat{P}\hat{\mathcal{L}}_{1}(\rho_{0}+\hat{P}\rho_{1})=0$. As $\rho_{0}$ is by assumption a unique solution, we find $\hat{P}\rho_{1}=0$. Therefore $\rho_{1}$ has no component in $\hat{P}$ space.

As $\hat{P}\mathbb{1}=0$, the identity matrix is in the $\hat{Q}$ space, $\hat{Q}\mathbb{1}=\mathbb{1}$. This seems to be a problem as we have to invert $\hat{Q}\hat{\mathcal{L}}_{0}\hat{Q}$ and $\hat{\mathcal{L}}_{0}\mathbb{1}=0$ for all $\hat{\mathcal{L}}_{0}$ describing unitary evolution . In all of our formulas, however, the inverse of $\hat{Q}\hat{\mathcal{L}}_{0}\hat{Q}$ is only applied to traceless density matrices, therefore no singularity arises from this zero mode. For numerical implementations one can simply add a ’mass term’ $m_{0}$ to this zero mode by replacing $\hat{\mathcal{L}}_{0}\to\hat{\mathcal{L}}_{0}+m_{0}\hat{P}_{\mathbb{1}}$ where $\hat{P}_{\mathbb{1}}$ is the superoperator projecting onto the identity matrix defined by $\hat{P}_{\mathbb{1}}X=\mathbb{1}\,\text{Tr}[X]/\text{Tr}[\mathbb{1}]$.

Appendix D Correlation functions

Higher order corrections $\delta\rho$ to $\rho_{0}$, as derived from our perturbation theory, can also be expressed in terms of multi-time correlation functions. In the main text we have already introduced

[TABLE]

however, to express higher order $\mathcal{O}(\epsilon^{n})$ contributions one needs also more complicated correlation functions

[TABLE]

In order to write them in a compact form with an explicit time ordering we introduce a general notation of Lindblad super-operators in the interacting picture,

[TABLE]

where Lindblad operators $L_{\alpha}(t)=e^{iH_{0}t}\,L_{\alpha}\,e^{-iH_{0}t}$ are evolved with respect to $H_{0}$.

Using this notation the lowest order correlation functions can be, according to Eq. (69), written as regularized time integrals

[TABLE]

Note that the time ordering of (super-)operators is such that it is most naturally represented using Keldysh formalism.

*First order

*Using the correlation functions defined above and the equality \text{Tr}\big{[}C_{k}\hat{\mathcal{L}}_{1}\hat{Q}\hat{\mathcal{L}}_{0}^{-1}\hat{Q}\hat{\mathcal{L}}_{1}\rho_{0}\big{]}=\text{Tr}\big{[}C_{k}\hat{\mathcal{L}}_{1}\hat{\mathcal{L}}_{0}^{-1}\hat{\mathcal{L}}_{1}\rho_{0}\big{]}, following from $\hat{P}\hat{\mathcal{L}}_{1}\rho_{0}=0$ for properly chosen $\rho_{0}$, the linear contribution to the expectation values of conserved quantities, Eq. (34), gets the form

[TABLE]

*Second order

*Second order $\mathcal{O}(\epsilon^{2})$ in the tangential part of $\delta\rho$ (relevant for expectation values of conserved quantities) has two contributions, $\delta\rho_{2,\parallel}=\delta\rho^{(1)}_{2,\parallel}+\delta\rho^{(2)}_{2,\parallel}$

[TABLE]

which we obtain by inserting $\hat{Q}=\hat{1}-\hat{P}$ under condition $\hat{P}\hat{\mathcal{L}}_{1}\rho_{0}=0$ and using that $\hat{\mathcal{L}}_{0}^{-1}$ cannot change the subspace, i.e., $\hat{P}\hat{\mathcal{L}}_{0}^{-1}=\hat{\mathcal{L}}_{0}^{-1}\hat{P}=\hat{P}\hat{\mathcal{L}}_{0}^{-1}\hat{P}$. The two contributions then add up to

[TABLE]

Using this expression it is straightforward to write the $\mathcal{O}(\epsilon^{2})$ contribution to the expectation value of conserved quantities in the steady state in terms of correlation functions. One should notice that each $\hat{P}$ ’cuts’ the expression in the sense that content between to consequential $\hat{P}$ corresponds to one correlation function. The inverses $(\hat{P}\hat{\mathcal{L}}_{1}\hat{P})^{-1}$ are treated as in Eq. (71), contribute $M^{-1}$ components to the final expression. Using the definition of projector $\hat{P}$, Eq. (16), the definition of the inverse $(\hat{P}\hat{\mathcal{L}}_{1}\hat{P})^{-1}$, Eq. (71), and definitions of multi-time correlation functions, Eq. (75), one obtains the following expression

[TABLE]

*Lehmann representation

*All expressions obtained from our perturbation theory can be calculated in a straightforward way using eigenstates and eigenenergies of $H_{0}$. As Lehmann representations are perhaps less common in the case of a Lindblad dynamics we show, for example, the Lehmann representation of $N_{j}^{(1)}$, Eq. (77), needed to calculate $\langle C_{i}\rangle_{1}$. By inserting the identity operator $\sum_{m}|m\rangle\langle m|$ between all operators and super-operators one obtains

[TABLE]

where $m,n,p,r$ run over the eigenstates $|m\rangle$, $H_{0}|m\rangle=E_{m}^{0}|m\rangle$, while $\alpha,\beta$ run over the Lindblad operators defined in Eq. (4). We use notation $\langle m|L_{\alpha}|n\rangle=L^{(\alpha)}_{mn}$ and $\langle m|\rho_{0}|m\rangle=\rho_{m}^{0}$, $\langle m|C_{j}|m\rangle=C^{(j)}_{m}$ where we assumed that the conserved quantities $C_{i}$ are diagonal operators in the eigenbasis of $H_{0}$. As $\hat{\mathcal{L}}_{0}|n\rangle\langle m|=-i(E_{n}^{0}-E_{m}^{0})|n\rangle\langle m|$, the inverse is obtained as $\hat{\mathcal{L}}_{0}^{-1}|n\rangle\langle m|=\frac{i}{E_{n}^{0}-E_{m}^{0}-i\eta}|n\rangle\langle m|$ using the regularization of appendix A. Note that $\eta$ should always be chosen to be larger than the level spacing if one is interested in systems in the thermodynamic limit.

Appendix E Unitary driving

In the case of unitary driving, the projection of the perturbation onto the $P$ space vanishes, $\hat{P}\hat{\mathcal{L}}_{1}\hat{P}=0$, and therefore also $\hat{P}(\hat{\mathcal{L}}_{0}+\hat{\mathcal{L}}_{1})\hat{P}=0$. In this case, the formulas fixing the density matrix to zeroth order in the perturbations, Eq. (11), and the perturbation theory, presented in Sec. III.3 have to be modified compared to the Lindblad dynamics where the inverse of $\hat{P}\hat{\mathcal{L}}_{1}\hat{P}$ within $P$ space is well defined.

E.1 Projections of the inverse

Different projections of full inverse $(\hat{P}+\hat{Q})\hat{\mathcal{L}}^{-1}(\hat{P}+\hat{Q})$ can, for example, be obtained using the transformations

[TABLE]

where $\hat{U}$ and $\hat{V}$ are chosen in such a way that they transform $\hat{\mathcal{L}}$ into a block-diagonal form (in $\hat{P}$ and $\hat{Q}$ space),

[TABLE]

The inverse $\hat{\mathcal{L}}^{-1}$ is then obtained from

[TABLE]

where

[TABLE]

Using the definition of $\hat{U}$ and $\hat{V}$ we obtain for perturbations with $\hat{P}\hat{\mathcal{L}}_{1}\hat{P}=0$ the following expressions for the projections

[TABLE]

$(\hat{P}\hat{\mathfrak{L}}_{2}\hat{P})^{-1}$ and $(\hat{Q}\hat{\mathcal{L}}\hat{Q})^{-1}$ should be interpreted as inverses within $P$ and $Q$ subspace, respectively, while the left hand sides of Eqs. (91) correspond to the projections of the full (properly regularized) Liouville inverse onto $P$ and $Q$ subspace. Note that the last equation does not imply that $\hat{Q}(\hat{Q}\hat{\mathcal{L}}\hat{Q})^{-1}\hat{Q}$ vanishes (which is finite for $\epsilon\to 0$).

A direct consequence of the equations given above is that in the limit of small $\epsilon$, the inverse in the $P$ sector is given through $\hat{\mathcal{L}}_{2}$ defined in Eq. (50),

[TABLE]

Note that $\hat{\mathcal{L}}_{2}$ is obtained from $\hat{\mathfrak{L}}_{2}$ just by replacing $\hat{Q}\hat{\mathcal{L}}\hat{Q}$ by $\hat{Q}\hat{\mathcal{L}}_{0}\hat{Q}$ in the inverse. $\hat{P}\hat{\mathcal{L}}_{2}\hat{P}$ takes over the role of an effective Lindblad super-operator in $P$ space.

E.2 Zeroth order expansion point

Here we show that the condition (11),

[TABLE]

does give the correct reference point $\rho_{0}$ for situations where $\hat{P}\hat{\mathcal{L}}_{1}\hat{P}=0$. One way to show this is to use the perturbative analysis provided in Sec. III.3 where the power counting of diagrams worked only if the correct reference point was chosen. Below we give a more direct argument.

As described in the main text, Eq. (93) is obtained from the requirement that the dominant contribution to the time-averaged expectation value of conserved quantities must vanish,

[TABLE]

where $\hat{P}$ is used to extract the non-oscillatory component and $\text{Tr}[C_{i}\hat{P}\hat{\mathcal{L}}_{1}\rho_{0}]=0$ due to the cyclicity of trace. As discussed in the main text, the starting point is the exact formula for $\delta\rho$, Eq. (6) and the formula Eq. (10) which directly leads to

[TABLE]

If we use only $\delta\rho^{(I)}$ in Eq. (94), then Eq. (93) follows immediately. Equivalently, the condition Eq. (93) implies that the contribution from $\delta\rho^{(I)}$ vanishes in Eq. (94). In the following we will show that Eq. (93) also implies that the contribution from $\delta\rho^{(II)}$ to Eq. (94) vanishes, which is less obvious, and a useful consistency check.

Plugging $\delta\rho^{(II)}$ into Eq. (94) one finds

[TABLE]

where the third line differs from the first one by two extra $\hat{Q}$ super-operators enclosing $\hat{\mathcal{L}}^{-1}$. The first one can be inserted because we consider the case $\hat{P}\hat{\mathcal{L}}_{1}\hat{P}=0$ and therefore $\hat{P}\hat{\mathcal{L}}_{1}=\hat{P}\hat{\mathcal{L}}_{1}\hat{Q}$. The second one can be used as a consequence of Eq. (93), which states that the $\hat{P}$ projection of the operator to the right of $\hat{\mathcal{L}}^{-1}$ vanishes. Finally, we can use that $\hat{Q}\hat{\mathcal{L}}^{-1}\hat{Q}=0$, see Eq. (91), to prove that the whole expression vanishes.

To finish our argument, we still have to show that $\delta\rho$ is small for $\epsilon\to 0$ provided that Eq. (91) holds, which can be done using similar arguments as above. First, the combination $\hat{\mathcal{L}}_{0}^{-1}\hat{\mathcal{L}}_{1}\rho_{0}=\hat{\mathcal{L}}_{0}^{-1}\hat{Q}\hat{\mathcal{L}}_{1}\rho_{0}$ is non-singular for $\hat{P}\hat{\mathcal{L}}_{1}\hat{P}=0$ which implies that $\delta\rho^{(I)}\sim\mathcal{O}(\epsilon)$. Second, we used already above that $\delta\rho^{(II)}=\hat{\mathcal{L}}^{-1}\hat{\mathcal{L}}_{1}\hat{\mathcal{L}}_{0}^{-1}\hat{\mathcal{L}}_{1}\rho_{0}=\hat{\mathcal{L}}^{-1}\hat{Q}\hat{\mathcal{L}}_{1}\hat{\mathcal{L}}_{0}^{-1}\hat{\mathcal{L}}_{1}\rho_{0}$. As $\hat{Q}\hat{\mathcal{L}}^{-1}\hat{Q}=0$ and $\hat{P}\hat{\mathcal{L}}^{-1}\hat{Q}\sim\mathcal{O}(1/\epsilon)$ it follows immediately that also $\delta\rho^{(II)}\sim\mathcal{O}(\epsilon)$ which concludes the derivation of Eq. (93).

E.3 Monochromatic driving

Some of the diagrams depicted in Fig. 2 vanish when a perturbation which contains only oscillations with a single frequency $\omega$ is considered. In this case $\hat{P}(\hat{\mathcal{L}}_{1})^{3}\hat{P}=0$ since each application of $\hat{\mathcal{L}}_{1}$ changes the Floquet index $n$ of $\rho^{(n)}$, Eq. (12), by $\pm 1$. Consequently, for example, the $\mathcal{O}(\epsilon)$ diagram in the $(\parallel)$ sector vanishes and conservation laws are only changed by processes of order $\mathcal{O}(\epsilon^{2})$, represented by the diagram in the second line of Fig. 2 (note that the third one vanishes as well). The situation is different in the presence of driving with higher harmonics, e.g. $\cos(\omega t)$ and $\cos(2\omega t)$, when $\hat{P}(\hat{\mathcal{L}}_{1})^{3}\hat{P}\neq 0$. In this case the dominant stationary correction to the expectation value of $C_{i}$ is of $\mathcal{O}(\epsilon)$, as shown by the diagram in the first line in the $(\parallel)$ sector.