Dissipative flow equations

We generalize the theory of flow equations to open quantum systems focusing on Lindblad master equations. We introduce and discuss three different generators of the flow that transform a linear non-Hermitian operator into a diagonal one. We first test our dissipative flow equations on a generic matrix and on a physical problem with a driven-dissipative single fermionic mode. We then move to problems with many fermionic modes and discuss the interplay between coherent (disordered) dynamics and localized losses. Our method can also be applied to non-Hermitian Hamiltonians.

The case of a one-dimensional system with a linear spectrum 31 References 33

Introduction
The study of many-body quantum physics has been recently challenged by the appearance of an increasing number of experimental platforms where genuine quantum phenomena take place notwithstanding the presence of an environment and of dissipation. Exciton polaritons [1,2], lossy atomic and molecular gases [3], cavity and circuit QED arrays [4,5], arrays of trapped ions [6] and Rydberg atoms [7], are only few prominent examples of a long list. Whereas the notion of equilibrium has been a fruitful guide to the development of standard manybody physics, these setups are inherently out of equilibrium and their modelisation requires the introduction of a Lindblad master equation that describes in an effective way the weak coupling to a bath under the Markov approximation [8].
Several methods for addressing this dissipative out-of-equilibrium dynamics have been proposed, based for instance on quantum trajectories [9], tensor networks [10,11], extensions to mean field theories [12,13], and machine learning [14][15][16][17]; yet, the solution of many-body physics for open quantum systems remains a formidable task. Clearly, techniques developed in the framework of Hamiltonian closed systems are a continuous source of inspiration for novel developments, and in this article we present the generalization of one such technique, the so-called flow equations [18], to the dissipative framework.
The method of flow equations has been independently developed by Wegner in the context of condensed-matter systems [19], and by G lazek and Wilson in the context of high-energy physics [20,21]. The main idea is the search for a parameter-dependent unitary transformation that transforms the Hamiltonian into a diagonal operator where eigenvalues can be easily read out. The approach has been successfully applied to several problems; within condensed-matter physics we can briefly mention Kondo and impurity problems [18] or quantum quenches in the Fermi-Hubbard model [22], quantum chemistry [23,24], and more recently many-body localisation [25][26][27][28][29][30][31][32][33][34].
In this article we develop the theory of dissipative flow equations. Technically, this requires to work with generators of the real-time dynamics that are not Hermitian, whereas the established theory relied significantly on the fact that Hamiltonians are Hermitian. We propose three different kinds of flow equations that are inspired by the original contributions by F. J. Wegner [19] and S. R. White [23]. After elaborating on the links with the dissipative Schrieffer-Wolff transformation [35] (for which we present a new derivation), we show how to employ the method to infer stationary and time-dependent properties of the dynamics. We discuss several examples, three of which deal with fermionic systems, where the study of the eigenvalues of the generator of the dissipative dynamics is particularly interesting. The formulation of the flow equation for fermions requires the use of the superoperator fermionic formalism [36,37], which is briefly reviewed in an appendix.
It is interesting to stress that attempts to using the flow equations for studying dissipative quantum systems have already appeared in the literature [38][39][40][41][42][43][44][45][46]. However, these approaches have typically described the global unitary dynamics of the coupled system and environment, rather than only focusing on the system, as we are proposing here. Our goal is not to follow a microscopic path, but rather to start from the beginning with a dynamics that focuses only on the system and takes into account the bath in an effective way. For this reason, our work will mainly focus on the Lindblad master equation, which is the most generic way of describing the dynamics of a system coupled to a Markovian environment. However, the method can also be used for non-Hermitian Hamiltonians.
Before concluding this introduction, it is important to stress the long-term motivation of this study. In this article we apply the dissipative flow equations only to quadratic fermionic systems, for which well-developed techniques already exist for solving the dynamics. While they perfectly serve as a benchmark for our novel method, it is also clear that our method cannot compete with them in any respect. We see this article as a first study in the direction of applying the dissipative flow equations to interacting systems which cannot be solved exactly and where approximations are necessary. Our perspective is the study of the dissipative flow equations in this context, where they could generate a new set of approximations and lead to novel solutions in a renormalization-group-like spirit (see Ref. [18] for a similar discussion in the Hamiltonian case).
The article is organized in two parts; in the former we present our theory of dissipative flow equations. In particular, in Sec. 2 we introduce the main general framework, whereas in Sec. 3 we present the details of three generators of the flow that accomplish the task of diagonalizing the Lindblad master equation in the long-flow limit. The second part is devoted to the discussion of several examples where we compare our approach with results obtained using more established techniques. In Sec. 4 we test our method on the diagonalization of a generic non-Hermitian matrix. We then move to physically-motivated problems with fermions and in Sec. 5 we discuss the problem of a single fermionic mode coupled to an environment inducing losses and gain. We then consider the problem of many fermionic modes in the presence of a localised source of losses, without disorder (Sec. 6) and with disorder (Sec. 7). This is also the occasion to discuss the dissipative flow equations in momentum space (Sec. 6) and in real space (Sec. 7). Our conclusions are presented in Sec. 8. The article is concluded by three appendices on the dissipative Schrieffer-Wolff transformation (Appendix A), on the superoperator formalism for fermions (Appendix B) and on the dissipative scattering model (Appendix C).

Definitions
We study the dynamics of an open quantum system in contact with a reservoir within the framework of the Markovian Lindblad master equation: Here, ρ(t) is the density matrix of the system at time t, H is the Hamiltonian of the system and L α are the quantum-jump operators effectively describing the coupling to the environment. The superoperator L is linear and not Hermitian: its eigenvalues λ are complex and such that [λ] ≤ 0; due to its specific form, if λ is eigenvalue then λ * is also eigenvalue. The eigenvalues determine the normal decaying modes of the dynamics: indeed, if L[τ ] = λτ , then τ (t) = τ exp[λt]. The zero eigenvalue λ = 0 is particularly important because its eigenvectors represent stationary states of the dissipative dynamics: if L[ρ 0 ] = 0 then ρ(t) = ρ 0 . Given a Lindbladian L, we look for a parameter-dependent invertible transformation S( ) with ∈ R + such that S(0) = I and such that becomes diagonal in the limit → +∞. We parametrize the invertible transformation introducing a generator η( ) which is a generic matrix, so that: Accordingly, it follows that where the second equality has been reported for later convenience. In Sec. 3 we will present three generators that diagonalize L in the infinite-flow limit.

Non-Hermitian Hamiltonians
This approach can be extended to non-Hermitian Hamiltonians, which constitute the other main theoretical tool for describing open quantum systems. In this case the dynamics is described by: whereH is a non-Hermitian operator; its eigenvalues have a clear physical importance: real parts are related to energies and imaginary parts to gain and loss rates. The discussion proposed in the previous paragraph can be easily adapted to this situation by applying the parameter-dependent invertible transformation S( ) toH in order to define a non-Hermitian operatorH( ) that is diagonal in the infinite-flow limit. Although in the article we will explicitly consider only Lindblad master equations, all results can be easily remapped to the framework of non-Hermitian Hamiltonians; indeed, the master equations discussed in Secs. 6 and 7 require the application of the flow-equation technique to two non-Hermitian Hamiltonians.

Invariants of the flow
We now identify quantities that do not change during the flow L( ); the simplest example is its characteristic polynomial: for this reason, it is an invariant of the flow. As a consequence, every eigenvalue of the matrix L( ) is an invariant of the flow. In the following, it will be useful to use a set of quantities I n (n ∈ N + ) that do not depend on : It is worth mentioning that differently from the Hamiltonian setting for a generic Lindbladian L, its 2-norm L( ) 2 2 = tr L( ) † L( ) is not an invariant of the flow.

Observables
Stationary-state properties The theory presented above gives direct access to the spectrum of L, but in addition to that, one might be interested in computing the average of an observable O over the stationary-state density matrix ρ 0 (we are here assuming that it is unique, but the degenerate case is treated in the same way): We can rewrite the latter expression as where This equation, valid at any , becomes particularly useful for → ∞ where the Lindbladian is diagonal and the stationary state ρ 0,∞ can be readily obtained: However, the use of Eq. (11) requires the parallel evolution of the operator from = 0 up to → ∞ using the same flow: Roughly speaking, the operator O has to be transformed in the basis where the Lindbladian is diagonal. The alternative interpretation is of course equally valid, namely, in order to evaluate Eq. (8) one needs to obtain ρ 0 from the knowledge of ρ 0,∞ , which amounts to evolve the latter backwards in flow: ρ 0 = lim →∞ S −1 ( ) ρ 0.∞ S( ).

Time evolution
Another interesting question concerns the real-time dynamics of the expectation value of a given operator: with ρ(t) = e Lt [ρ(0)] where e Lt is the exponential of the superoperator tL and ρ(0) is the density matrix at the initial time. We can plug in the identity S( ) −1 S( ) and write the time evolution in the running basis where we have used the following property of the operator exponential: S( )e Lt S( ) −1 = e L( )t . Once more, the expression simplifies for → ∞, since in this basis the Lindbladian is diagonal and the solution of the time evolution becomes trivial: As in Eq. (11), the solution of this problem requires the knowledge of the operator O( ) in the → ∞ basis, to be obtained with Eq. (12). In addition to that, also the initial density matrix has to be written in such basis by solving the differential equation: With respect to the Hamiltonian dynamics [18,47], which requires forward evolution, realtime propagation and then backward flow evolution, in the dissipative case one can obtain the real time dynamics with a forward flow evolution of the initial density matrix and of the operator.

Generators of the flow
In this section we present three generators of the flow that accomplish the goal described in the previous section, namely the diagonalisation of the Lindbladian L in the infinite-flow limit → ∞. It will be useful to address separately the diagonal and the off-diagonal parts of the Lindbladian, that we dub respectively D( ) and V( ). Obviously, this choice is basis dependent and L( ) = D( ) + V( ). In the following, in order to show that the generators diagonalize the Lindbladian, we will in particular focus on V 2 2 and on I 2 .

First choice of the generator
Inspired by the original work by Wegner [18,19], we propose the following generator of the flow: In order to show that it diagonalize L, we now show that it induces a flow such that V 2 2 ≥ 0 cannot increase as a function of .

Proof
We first compute the derivative with respect to the flow parameter as follows: Let us now discuss two interesting identities, which follow from the fact that V( ) multiplied by any diagonal operator is an off-diagonal operator: Concerning the second identity, it is important to stress that dD( ) At this point, we can use the second of the identities (19) and write: If we take the definition in Eq. (17) we obtain:

Possible issues
Similarly to the Hamiltonian case, Eq. (21) does not rule out the possibility of a flow that does not start. This happens for instance when D(0) is the zero matrix. This issue can be circumvented numerically by applying at the beginning of the dynamics a random invertible operator R 0 to the Lindbldian: L(0) → R 0 L(0)R −1 0 .

Second choice of the generator
As a second generator, we propose the following one, which is a slight modification of the previous one: For later convenience, we have reported the matrix elements η nk of η in a basis of choice, using the notation d nn for the diagonal matrix elements of D and V nk for the off-diagonal matrix elements of V. We do not have a proof that this generator does the desired job apart from the numerical evidence reported in the next sections, which also shows that it is more efficient than the previous generator. We can however present some considerations on its convergence based on perturbative arguments.

Perturbative solution
Let us begin by writing the flow equations for all matrix elements: The solution of this set of equation looks rather demanding, but it can be simplified if we assume that at the initial time the following condition holds: We introduce the small parameter ξ proportional to the off-diagonal part of the Lindbladian. We now expand each matrix element as: with initial conditions: The flow equations for each term of the expansion are obtained by comparing terms of the same power of ξ: At the lowest non-trivial order for each term, the equations are solved by: We obtain a correction that is thus completely consistent with perturbation theory (although we remark that flow equations are a technique that is not perturbative in spirit). It is interesting to observe that in this limit the off-diagonal matrix elements V nk decay to zero with a typical flow scale that is proportional to their eigenvalue difference |d nn (0) − d kk (0)| 2 ∼ |d nn ( ) − d kk ( )| 2 . This means that an interpretation of the flow as a generalized renormalization group where the decimation suppresses terms that are distant in energies might be possible as in the Hamiltonian case [18].

Third choice of the generator
Let us now consider the following generator: which resembles the one proposed by S. White in the Hamiltonian setting [23].

Perturbative limit and dissipative Schrieffer-Wolff transformation
Even if we are going to present a proof that shows that this generator works also in nonperturbative regimes, we motivate its introduction considering the case of an off-diagonal part of L( ) that is a small perturbation; Its strength is controled by a small parameter ξ.
We thus propose an expansion of the generator η( ) in powers of ξ: No zero-th order term is introduced because for ξ = 0 the matrix is diagonal and the generator η( ) must be equal to zero, so that S( ) = I. At first order in ξ, the flow equation reads: In order to kill the off-diagonal part of L( ), it is reasonable to ask: This latter equation defines the third generator of the dynamics, whose matrix elements were presented in Eq. (29). It is important to observe that Eq. (32) also allows us to establish a link with the theory of dissipative Schrieffer-Wolff transformation, whose form was first derived in Ref. [35], and for which we present another derivation in Appendix A. In this approach, one looks for an invertible transformation that rotates the Hilbert space and decouples subspaces related to different eigenvalues of the unperturbed system. It is customary to introduce a generator also for the Schrieffer-Wolff transformation, and at first order in the strength of the perturbation one obtains exactly Eq. (32) (see Eq. (94) in Appendix A). Thus, the third generator implements a flow that reproduces the action of the Schrieffer-Wolff transformation.

Proof
We consider the second invariant of the flow: We now study the derivative with respect to of I off 2 : using the fact that d d I 2 = 0 we obtain d d This differential equation is solved by: and shows a very efficient decay to 0 with a typical flow-scale that does not depend on any system parameter. This scaling is expected to be particularly useful in numerical applications.

Possible issues
The fact that I off 2 is equal to zero does not mathematically guarantee that the off-diagonal part of the Lindbladian matrix is equal to zero; one possible example is constituted by an initial matrix that is triangular: in this case I off 2 (0) = 0. This issue, that is not present in the Hamiltonian case, can be circumvented numerically by applying at the beginning of the dynamics a random invertible operator R 0 to the Lindbladian:

First example: Numerical tests on a generic matrix
As a first example, we apply the flow equations to a generic non-Hermitian complex matrix A with size 15 × 15; the real and imaginary parts of the matrix elements are randomly taken from the uniform distribution in the range [−1, 1]. We have verified that the main qualitative features of the data that we present in this Section do not depend on the specific choice of A.
We implement the numerical solution by means of a 5-th order Runge-Kutta algorithm (Butcher's scheme); the flow step is d = 10 −3 and N step = 15000, so that max = 15. During the flow, we set to zero the off-diagonal matrix elements whenever their absolute value is less than the 1% of the initial value, as routinely done in Hamiltonian setting [18]. We verify that at the end of the simulation the matrix A( max ) is well approximated by a diagonal matrix; the real and imaginary parts of the diagonal elements are compared with the values obtained by standard matrix diagonalization and we benchmark the convergence using the discrepancy We also study the invariants I n , which should be conserved by the flow, and in particular focus on the relative error δI n = |I n ( max ) − I n ( = 0)|/|I n ( = 0)|.
We study the dissipative flow equations using the first, the second and the third generator, for which the real and imaginary parts of the diagonal matrix elements are plotted in Figs. 1, 2 and 3, respectively. We plot the real and imaginary parts of the diagonal elements; a comparison with the exact results obtained with standard diagonalization routines is also presented.
An important property of the third generator is the fact that the quantity I off 2 ( ) satisfies the differential equation (35); as a consequence, the convergence to the diagonal form is very  17); dashed red lines represent the eigenvalues obtained by standard diagonalization routines. For this specific case, the discrepancy is ∆ 8.5 · 10 −3 ; we have verified that the relative error of the invariants δI n with 1 ≤ n ≤ 15 never exceeds 10 −5 , the worst case being I 15 .  22); dashed red lines represent the eigenvalues obtained by standard diagonalization routines. For this specific case, the discrepancy is ∆ 7.3 · 10 −3 ; we have verified that the relative error of the invariants δI n with 1 ≤ n ≤ 15 never exceeds 10 −8 , the worst case being I 15 .  29); dashed red lines represent the eigenvalues obtained by standard diagonalization routines. For this specific case, the discrepancy is ∆ 1.9 · 10 −7 ; we have verified that the relative error of the invariants δI n with 1 ≤ n ≤ 15 never exceeds 10 −7 , the worst case being I 15 . Figure 4: Decay during the flow of |I off 2 ( )| for the first (blue), second (orange) and third (green) generator. In this latter case, we recover the behaviour predicted in formula (35). efficient. The exponential decay I off 2 ( ) = I off 2 (0)e −2 has been numerically verified, see Fig. 4. A similar numerical calculation with the first and second generator shows a significantly slower decay, although in both cases we eventually recover an exponential law (see Fig. 4). We thus conclude that for numerical applications the third generator is the most efficient choice.

Second example: Fermionic mode with losses and gain
In this and following Sections we test the the dissipative flow equations with several physical examples for which exact solutions can be easily obtained analytically and numerically. We will show that the theory is correct and discuss the specific properties of each generator, highlighting advantages and disadvantages. We begin in this Section with the study of a single fermionic mode coupled to a bath inducing incoherent losses and gain. The interest of this simple example relies on the fact that the flow equations can be solved analytically in several limits.

The problem
We study a single fermionic mode with energy ε coupled to a bath; fermions are lost at a rate Γ 1 and gained at a rate Γ 2 . We introduce the canonical fermionic operatorsĉ andĉ † , and model the dynamics of the site with the following master equation: We propose the study of this model using the formalism of superfermion representation, that is presented in detail in Refs. [36,37] and that is also reviewed in Appendix B. Roughly speaking, the idea is to treat the density matrix as a vector of an appropriate Hilbert space isomorphic to H ⊗ H, where H is the two-dimensional Hilbert space of a single fermionic mode: ρ(t) → |ρ(t) . We subsequently need to introduce superoperators c andc that describe the action of fermionic operators on the left or on the right of the density matrix; they square to zero (c 2 = 0 andc †2 = 0) and satisfy mutual canonical anticommutation relations {c ( †) ,c ( †) } = 0. To all effects, this formalism describes the dynamics of a single fermionic mode coupled to a bath as a fermionic two-mode problem. This approach has the great advantage to allow to represent a quadratic fermionic master equation as the one in Eq. (36) as a matrix, and to link the normal decaying modes to its eigenvalues.
A detailed discussion of model (36) using superoperators is reported in Ref. [36]; here we briefly mention some relevant aspects. According to this formalism, (36) can be cast in the following form: where L is an operator that is quadratic in the fermionic superoperators: Note that roughly this expression can be obtained from (36) by using ac operator each time the fermions act to the right of the density matrix. Since it is a quadratic operator, we can put it into a 2 × 2 matrix form by exploiting the aforementioned anticommutation relations. We obtain: It is possible to diagonalize the matrix with an invertible and non-unitary transformation (not specified for the moment) and introducing the operators d, D † ,D andd † , so that: The steady state is annihilated by both d andD operators: and is defined by this relation. The normal decaying modes are obtained by applying the D † andd † operators onto the steady state; the corresponding Lindbladian eigenvalue determines their time evolution. In this case, for instance, by looking at the eigenvalues of L, which are λ ± = ε ± i 2 (Γ 1 + Γ 2 ), we can infer that decays take place according to the typical time τ = (Γ 1 + Γ 2 ) −1 .

Three approaches with dissipative flow equations
We now apply the techniques of the flow equations to the matrix representation of the superoperator L in Eq. (39) in order to put it in diagonal form and extract the two eigenvalues λ ± .
We parametrize the Lindbladian matrix L( ) in the following way:

Invariants of the flow
Since the matrix is 2 × 2, there are only two independent invariants of the flow: From the invariance of these expressions, we obtain that: First generator The first generator of the dynamics reads: and its commutator with the Lindbladian matrix: This leads to a set of coupled non-linear differential equations: By using the invariants listed above, we can reduce the three equations to a single one: We do not give an explicit (and useless) analytical solution here; it is however clear that the stationary values of α( ) areᾱ 1,2,3 = {±(Γ 1 + Γ 2 )/2, 0}, among which we find the correct value. The specific value is determined by the initial conditions. Using the second invariant, we can conclude that the stationary valuesᾱ 1,2 = ±(Γ 2 + Γ 2 )/2 are accompanied by the stationary valueμ 1μ2 = 0; Eqs. (47) additionally say that eachμ 1 andμ 2 should be equal to zero. Vice-versa, the stationary valueᾱ 3 = 0 is accompanied byμ 1μ2 = (Γ 1 + Γ 2 ) 2 /4. In order to test the stability of the three stationary solutions of the flow, we consider Eq. (48), which we write in the form d d α( ) = f (α). We then evaluate f (α 1,2,3 ) for the three stationary values and obtain: Given the presence of the second invariant, that links the flow of α( ) to that of µ 1 ( )µ 2 ( ), we can conclude that ifᾱ 1,2 is a stable stationary values, then alsoμ 1 = 0 andμ 2 = 0 are stable stationary values. This simple analysis reveals that the stationary points with vanishining off-diagonal part are locally stable, and depending on the initial conditions, one of the two specific valueᾱ 1,2 is selected.

Second generator
The second generator of the dynamics reads: and its commutator with the Lindbladian matrix: This leads to the following set of coupled non-linear differential equations: These equations are very similar to those obtained with the first generator and reported in Eq. (47). The stationary values can be obtained with the invariants of the flow. This approach leads us again to Eq. (48), for which the analysis discussed above can be repeated.

Third generator
The third generator of the dynamics reads: and its commutator with the Lindbladian matrix: This leads to the following set of coupled non-linear differential equations: It is very easy to obtain µ i ( ) = µ i (0)e − and accordingly to verify property (35), concerning the flow-evolution of I off 2 : Using the second invariant, we can also obtain α( ): A simple analytical case: Γ 2 = 0 In order to shed more analytical light on the technique, we consider here the special case Γ 2 = 0. In order to study this special case, we re-parametrize L( ) setting µ 2 ( ) = 0; the matrix is now triangular, its eigenvalues can be directly read from the diagonal, which is not supposed to evolve during the flow. Indeed, from the second invariant we obtain α( ) = − Γ 1 2 . We notably simplifies the differential equations satisfied by µ 1 ( ); we list them here below for the three generators: Whereas the second and third equations are trivial, and are solved by two exponentials, in this limit it is possible to give a simple and analytical solution also to the former: This results, together with those presented in Sec. 4 highlight the fact that the third generator of the flow is the most efficient from a numerical viewpoint. On the other hand, the first and second generators share several similarities, and given that between the two the second is the most efficient and simplest, we disregard from now on the first generator of the flow.

Third example: Dissipative scattering model
We now move to the application of dissipative flow equations to a problem involving many fermionic modes.

The problem
We now discuss a gas of spinless fermions in a d-dimensional square box of volume L d in the presence of a loss mechanism that acts locally; the Lindblad master equation reads: (60) We consider in particular the case of a loss mechanism that is localized around x = 0, and choose the form: Γ(x) = γ δ(x). Note that γ has the dimension of a volume divided by a time. We furthermore assume that the Hamiltonian is single-particle H = k ε(k)ĉ † kĉ k and the full dynamics can be written in momentum space as follows: This master equation is quadratic in the fermionic operators and amenable to the treatment with fermionic superoperators recalled in Sec. 5 and detailed in Appendix B. We now Figure 5: Plot of the imaginary part of the eigenvalues of the dissipative scattering model for L = 201 considering 601 states (the parameter j Λ defined in the appendix is set to 300). The plot highlights the emergence of a strongly dissipative state for γ > 4v which is well described by the formula in (64). For γ < 4v the formula in 135 in Appendix C describes the imaginary part of the eigenvalue with real part equal to zero that evloves into the strongly dissipative state.
pass to the superoperator representation: i d dt |ρ(t) = L|ρ(t) with In matrix form: where H is a diagonal matrix with entries ε(k) and Λ 1 is a matrix with all matrix elements equal to γ/L d . Since the matrix is block triangular, in order to study its eigenvalues it is sufficient to look for the eigenvalues of the matrix M = H − i Λ 1 /2; in the following we will only concentrate on it. This of course prevents a correct reconstruction of the observables, as detailed in Sec. 2.2, but as long as the focus is on the spectrum it is a legitimate restriction. The problem has been discussed in several articles, see for instance Refs. [48][49][50][51][52]. Here we focus on a simplified situation, that of a one-dimensional setup (d = 1) with a linear dispersion relation: ε(k) = v 2π L j, where v is a velocity and j ∈ Z. In particular, we are interested in an important spectral feature of the model: for γ ≥ 4v, we observe the appearance of a strongly dissipative state, namely of an eigenvalue with real part equal to 0 and imaginary part much larger then that of the other eigenvalues: where Λ is an appropriate energy cutoff. This eigenvalue marks the appearance of a single state where almost all dissipation is concentrated, whereas all other ones are weakly dissipative (see Fig. 5). This separation of time-scales is typical of the quantum Zeno effect. In appendix C we present the analytical solution of this model and the necessary details.

Dissipative flow equations
We will use this example to discuss the dissipative flow equations in momentum space. Following the expression in Eq. (62), we propose the following parametrization for the flow equations: with the following initial conditions: We now divide the super-operator M ( ) into a diagonal and an off-diagonal part: and apply the theory that we have developed for the dissipative flow equations. The flow is characterized by several invariants I n , the first and the second read:

Second generator
We proceed to the solution of the problem using the second generator of the dynamics: whose explicit expression reads (we omit for brevity the dependence on ): where we have used the important identity: We now compute the commutator [η( ), M ( )] by splitting it into two parts; first the commutator with D( ) and then the commutator with V ( ): With this information we can write the flow equations for the coupling constants:

Third generator
According to the general definition, we propose the following generator of the dynamics: In order to compute the commutator [η( ), L( )], we can reemploy some of the results obtained for the second commutator, and obtain: We are now ready to write the flow equations for the third generator:

Numerical solutions
We have solved numerically the flow equations in (74) and in (77). In the following we present the numerical results obtained by solving the flow equations for the dissipative scattering model both with the second and third generator. We set γ = 5v and j Λ = 15 (so that we consider 31 states in total). The numerical algorithm that has been used to solve the system of coupled ordinary differential equations is a 5 th order Runge-Kutta (Butcher's scheme) with a flow step d = 10 −4 and N step = 10 5 for the second generator, whereas N step = 8.5 · 10 4 for the third generator.  for the second (light blue) and third (red) generator. In this latter case, we recover the behaviour predicted in formula (35). The second generator instead has a less uniform behaviour, although at long times it also show an exponential decay. 20 The results are summarized in Fig. 6. We first observe that both generators let the system converge to a diagonal form. In order to be more quantitative, the discrepancy ∆( max ) for the second generator equals 6.0 · 10 −3 , whereas for the third generator is 8.11 · 10 −11 . Concerning the conserved quantities, we observe that δI 31 equals 1.8·10 −2 for the second flow and 3.8·10 −6 for the third one. In Fig. 7 we show the flow of V( ) 2 2 and observe that the third generator produces a more effective decay to zero.
This study corroborates the previous conclusions on the fact that for numerical applications the third generator is more effective than the second one. On the other hand, the dynamics generated by the second generator is more uniform than that of the third generator. In this respect, we envision that the use of the second generator would be more useful in situations where an approximated treatment is necessary (e.g. for interacting non-quadratic systems).

Fourth example: Disordered dissipative scattering model
In this Section we continue our analysis of dissipative flow equations applied to a many-fermion problem and consider a one-dimensional lattice with a fermionic disordered tight-binding model: Here, h j takes random values and is uniformly distributed in the range [−W, W ]. We consider a localized loss source at the center of the lattice (j = 0) with loss rate γ, so that the master equation reads: d dt ρ(t) = − i [H, ρ(t)] + γ ĉ 0 ρ(t)ĉ † 0 − 1 2 {n 0 , ρ(t)} . In this example we aim at investigating the interplay between disorder and losses, and at discussing the emerging behaviour of the system as a dissipative insulator. This analysis can be performed by focusing on the size-scaling of the asymptotic decay rate, namely the longest typical decay time of the system, which is related to the spectrum of the Lindbladian. Similarly to what done in Sec. 6, we rewrite the master equation using the fermionic superoperators (see Appendix B) and link the spectrum of the master equation to the eigenvalues λ α of the matrix M = H − i Λ 1 /2, which represents the master equation in this formalism. The asymptotic decay rate reads: In Fig. 8 we show the asymptotic decay rate of the problem, obtained with exact diagonalization, for several values of the lattice length L and of the decay rate γ. Its scaling is exponential in the system size: Γ adr ∼ exp[−L] and this behaviour is solely dictated by the presence of disorder. Indeed, we verified that for a clean system the scaling is always algebraic L −α . Similarly to the model discussed in Sec. 6, the system also displays a strongly dissipative state for γ > 4J, both in the clean [49,51] and in the disordered case (see Fig. 8). The latter result, related to the disordered model, has not been thoroughly discussed yet and deserves further investigation.

Dissipative flow equations
We propose to study this model using flow equations that are formulated in real space: with initial conditions: The flow equations are not different from those presented in Sec. 6.2 and we limit ourselves here to writing the final results. For the second generator of the flow: For the third generator of the flow: We present a numerical solution of the flow dynamics; we use a 5-th order Runge-Kutta algorithm (Butcher's scheme) with adaptive flow steps keeping an estimated error below the threshold of 10 −16 according to Butcher's tableau.
In the top panels of Fig. 9 we show the flow of the imaginary part of the diagonal matrix elements g jj ( ) averaged over 10 4 disorder realizations obtained with the second and the third flow generators. In both cases we see a good convergence to the expected values, computed with standard linear-algebra numerical packages. This convergence is associated to the decay with the flow of the off-diagonal part of the matrix. The bottom panels of Fig. 9 show that with the second generator the decay is algebraic, whereas with the third generator the decay is exponential in .

Conclusion
In this paper we have generalized to open quantum systems the flow equations that have originally been developed by Wegner, G lazek and Wilson. Specifically, our work shows how to generalize the flow equations to operators that are not Hermitian, focusing in particular on fermionic Lindblad master equations. Altough we did not discuss it explicitly, our results can also be employed for non-Hermitian Hamiltonians. We have described three generators of the flow and have highlighted their peculiarities and strong points. In particular, we have argued that the third generator is the most suited for numerical applications.
The main perspective of this work is related to the believe that the second generator of the flow could find a fruitful application in novel approximations schemes aiming at the development of renormalization group-like approaches. We have indeed shown that, in the long flow limit, the off-diagonal matrix elements decay with a flow scale that depends on the difference of the two eigenvalues that they connect. This behaviour is reminiscent of decimation schemes proposed for renormalization groups of Hamiltonians, as discussed for instance in Ref. [18]. Whereas we have benchmarked our new method with four different examples, in all these cases efficient techniques for the solution of the dynamics exist. On the other hand, the interest of dissipative flow equations might reside in the development of novel truncation and approximation schemes for the treatment of interacting problems in the presence of dissipation, where instead we lack a general and well-established method to use.

A Dissipative Schrieffer-Wolff transformation
We present a new derivation of the Schrieffer-Wolff transformation for Lindbladian operators (and in general for non-Hermitian linear operators) based on the considerations presented in Ref. [53] for Hamiltonian perturbation theory. The same results have already been derived in Ref. [35] using the resolvant method detailed in Ref. [54].
We consider a quantum system whose dynamics is described by a Lindblad master equation and assume that the superoperator L can be written as the sum of two parts where ξ is a dimensionless quantity which plays the role of a perturbative parameter. We assume that L 0 can be easily diagonalized and that its eigenvalues λ αi can be grouped into well-separated sets labeled by α. Eigenvalues relative to different values of α are very different, but within each set α there is not exact degeneracy, so that an additional index i is necessary.
The right eigenvectors {|α, v i } i span the subspaces E (0) α , the left eigenvectors are noted |α, u j and the following relations hold: where P α is the projector onto the subspace E (0) α thanks to the relation α, u j |α, v i = δ ij .
The key idea of a Schrieffer-Wolff transformation is to obtain an effective operator L that is block-diagonal with respect to the subspaces E but that at the same time takes into account presence of ξL 1 , which in principle is not block-diagonal (if it were, there would be no need for this theory). This process will be implemented by an invertible transformation Q (remember that invertible transformation preserve the spectrum), which we write as: Using the Taylor series of the matrix exponential, it is possible to give another expression to this latter formula using nested commutators: In order to determine the matrix elements of η, we will employ a perturbative approach in ξ and expand η in powers of ξ: We do not include any zero-th order term because in the case ξ = 0 the operator is already block diagonal and the invertible transformation should be Q = I, which is what one obtains with η = 0. At this point we have to consider Eq (89) and compare terms of the same order in ξ. For the left-hand-side of the equation we introduce the notation: For terms at zero-th order in ξ we obtain L (0) = L 0 ; for those proportional to ξ we obtain the interesting equality: Exploiting the fact that the matrix elements of L (1) between two different manifolds must be zero, see Eq. (87), one obtains the equation: The latter equation determines the matrix elements of η (1) between two subspaces with different label α: For α = β, the matrix element is not unambiguously determined. This ambiguity follows from the fact that once the matrix η has been found, it is possible to construct an infinite number of other solutions by simply applying an arbitrary invertible transformation to Q that does not mix different subspaces. One possibility to remove this uncertainty is to impose that η has no matrix elements inside each manifolds: P α ηP α = 0, ∀α. For this reason, we set α, u j |η (1) |β, v i = 0.
We can now determine the matrix elements of L eff,α up to second order in ξ. For the zero-th order, we simply obtain P α L 0 P α . For the first order, instead, unsing Eq. (91) and observing that η (1) , L 0 has only matrix elements between states with different values of α, we obtain: P α VP α .

B Superoperator formalism for fermions
We briefly review the superoperator formalism for fermions (see Refs. [36,37] for more extensive discussions). It is also worth to mention that this formalism is also connected to the third-quantization formalism presented in Ref. [55], where a real fermion representation is preferred to a complex one. Section B.2 presents some remarks that we did not find explicitly written elsewhere.

B.1 Generalities
The spirit of the superoperator formalism is to represent density matrices ρ in a space H ⊗H, where H is the Hilbert space associated to the physical system (it could be a Fock space), whereasH is a space that is isomorphic to it. We introduce the basis {|n } n for H and the basis {|ñ } n forH and the isomorphism U : H →H that maps |n → |ñ . With this notation we can introduce the left vacuum vector : and the representation of the density matrix ρ = nm ρ nm |n m| as a vector of this space: The normalization of the density matrix tr[ρ] = 1 reads I|ρ = 1, whereas the expectation value of an observableÂ, defined as A = tr[ρÂ], reads A = I|Â ⊗Î|ρ .
We now consider the case where H is a fermionic Fock space, with L associated fermionic operatorsĉ m (without loss of generality, we do not consider explicitly spin) satisfying canonic anticommutation relations: In this case it is customary to define the left vacuum state in a slightly different way with respect to Eq. (99), and namely: |I = n 1 ,n 2 ...n L (i) n 1 +n 2 +...+n L |n 1 , n 2 , . . . n L | n 1 , n 2 , . . . n L The fermionic superoperators for the space H ⊗H are defined as follows: As a consequence, we can think of H ⊗H as an enlarged Fock space with 2L anticommuting modes. Once applied onto the left vacuum state (102), the c m andc m operators satisfy the fundamental tilde conjugation rules, that are crucial in all subsequent calculations: The definitions of thec m operators in Eq. (103) are crucial for ensuring the anticommuting properties in Eq. (104c). In principle they are a choice, but they are highly recommended because they allow to treat c m andc m on the same footing.

Matrix representation based on superoperators
We now investigate the writing of a Lindblad master equation (1) in superoperator representation assuming that the Hamiltonian is quadratic in the fermionic operators and that the jump operators are linear: Before continuing, we make the following physical assumption: we are only interested in the study of jump operators that either inject particles in the system or take them out. This means that for a fixed α, either the l 1αm or the l 2αm are all zeros. As a consequence, it will always be true that l 1αm l 2αn = 0. For later convenience, let us write: In the last line we have defined the Hermitian matrices Λ 1 and Λ 2 with matrix elements Λ 1mn = α l * αm l 1αn and Λ 2mn = α l αm l * 1αn . It is customary to represent the operator i L instead of L for its formal similarity with the Schrödinger equation for pure states. Since it is a linear operator, in the superoperator representation it will have quadratic matrix form i L[ρ] → L|ρ : In order to obtain this expression, we have directly promoted every original operatorĉ ( †) m to an operator c ( †) m acting on H ⊗H. Subsequently, we have used the fact that the density matrix ρ commutes with the parity operatorP = (−1) ĉ † mĉm (we also say that it is an even operator) and thus that it commutes with every operatorc ( †) m . With these simple rules, every term can be readily obtained. For instance, for what concerns the Hamiltonian dynamicsĤρ − ρĤ, the first term is easily recast in the superoperator language:Ĥρ|I = mn h mn c † m c n |ρ ; the second instead requires some manipulation: (109) Let us now attempt to put the operator L in Eq. (108) in matrix form: where M is a 2L × 2L complex matrix and K is a complex constant. We need to rewrite Eq. (108) as follows: so that the matrix M has the following block-diagonal form: and the matrices H, Λ 1 and Λ 2 are Hermitian; moreover We observe that the expression for L derived in the previous equations is a generalisation of Eq. (39) presented in the main text for a dissipative fermionic mode.

Diagonalization of a quadratic Lindblad master equation
The matrix M in Eq. (112) satisfies a strong symmetry: where I is the identity; this is a generic result that is true for any matrix of the form provided B and C are hermitian. As a consequence, M and M † have the same spectrum; indeed, if we look at the characteristic polynomial: Since the eigenvalues of M † are the complex conjugates of those of M , we obtain that if λ is an eigenvalue of M , then this is true also for λ * . This means in particular that either the eigenvalues are complex and come in pairs, or they are reals. We now make the assumption that the Jordan canonical form of the matrix M does not contain any nilpotent part; we are not aware of any physical problem in quantum physics where this matematical object played a role. For this reason, we assume that there is an invertible transformation S that puts M in diagonal form: Thanks to this, we can write: where the operators are defined through the matrix elements of S and S −1 : It is interesting to observe that when Λ 2 is a real matrix, in order to extract the spectrum {λ} it is not necessary to diagonalize the full matrix M . Indeed, we can rewrite the matrix M in Eq. (112) in the following compact way: We now propose the following unitary transformation: σ x → −σ y , σ y → −σ z , σ z → σ x and obtain the following matrix representation: If Λ 2 is a real matrix, M becomes a block-diagonal matrix and thus, in order to study the spectrum of M , it is sufficient to study the spectrum of H ± i 2 (Λ 1 + Λ 2 ). We exploit this possibility in Secs. 6 and 7, where Λ 2 = 0.

Back to a fermionic master equation
We conclude this section with a discussion of the physical meaning of Eq. (118). First of all, by exploiting the canonical anticommutation relations in (120), we rewrite it as: and for brevity we also introduce 0 = v 2π L . Let us begin by considering Eq. (128) and let us show that when λ is purely imaginary (we thus take [λ] = 0 and parametrize it as λ = iλ I ) it satisfies the first constraint. We fix a high energy cutoff Λ = 0 j Λ with j Λ 1 and such that −j Λ < j < j Λ and write that: It is important to consider the cutoff otherwise one would obtain that every λ of the form λ = ε(k) + iλ I would satisfy the first constraint. We continue with the second equation: The series can be analytically evaluated: where ψ(z) is the Digamma function. We now consider the large band-width limit, with |j Λ ± iλ I / 0 | 1, so that the following asymptotic expansion can be used: Note that by definition j Λ ± iλ I / 0 cannot lie on the negative real axis, where the expansion would be problematic. Concluding, we obtain the following equation for the eigenvalue: The equation can be simplified by considering that j Λ is larger than Λ I / 0 , so that we can approximate arctan x ∼ x. By introducing the variable x = πλ I / 0 , the equation reads: Although the equations has formally two solutions, for physical reasons we only retain the negative one. The peculiar property of this equation results from the fact that coth(x) is always smaller than −1 for x < 0. Thus, we can identify three regimes depending on whether γ/v is smaller than 4, larger than 4 or approximately 4. We discuss analytically two of them, whose results can be compared with exact numerics in Fig. 5.
The case γ/v 4 In this case 4v/γ 1 and the solution must satisfy |x| 1. The correction due to the arctan(x/(πj Λ )), that in this case can be approximated by x/(πj Λ ), is negligible and can be safely disregarded. The eigenvalue reads: The formula is well-defined only for γ < 4v and displays a divergence for γ → 4v − that we do not consider as physical because it is not in the regime of validity of the approximations. In the deep perturbative limit γ v the result gives: The case γ/v 4 In this case 4v/γ 1 and the solution must satisfy |x| 1. In this region, we can approximate coth(x) ∼ −1 and obtain: We thus obtain that the result is proportional to the band edge, with limiting value for γ/v → ∞ equal to −∞.