Closed hierarchy of correlations in Markovian open quantum systems

We study the Lindblad master equation in the space of operators and provide simple criteria for closeness of the hierarchy of equations for correlations. We separately consider the time evolution of closed and open systems and show that open systems satisfying the closeness conditions are not necessarily of Gaussian type. In addition, we show that dissipation can induce the closeness of the hierarchy of correlations in interacting quantum systems. As an example we study an interacting optomechanical model, the Fermi–Hubbard model, and the Rabi model, all coupled to a fine-tuned Markovian environment and obtain exact analytic expressions for the time evolution of two-point correlations.


Introduction
A complete description of nonequilibrium steady states is one of the main aims of the nonequilibrium statistical mechanics. In this regard exactly solvable models are of great importance. Their value lies in the opportunity to compare numerical simulations and an ability to recognize fundamental physical laws. In the classical setting steady states have been numerically and analytically investigated in stochastic lattice gas models [1] and in exactly solvable one-dimensional exclusion processes [2][3][4][5], where the steady state is framed in a matrix product form [6]. This ansatz has proven useful also in quantum setting, where it has been utilized to construct the nonequilibrium steady state of the boundary driven XXZ spin-1/2 chain [7][8][9]. Furthermore, the existing classes of Markovian many-body open quantum system models that can be exactly solved (by this we mean to analytically obtain the steady state) are strongly related to integrable closed quantum models, namely the quasi-free bosonic or fermionic models [10][11][12] and Yang-Baxter integrable models [13]. In addition, there are some solutions that do not evidently rely on the underlying integrability structure, e.g. the quantum symmetric exclusion process [14], the XX spin-1/2 chain with dephasing [15], and the noninteracting harmonic oscillator lattices in arbitrary dimensions with dephasing [16]. Exact long-time behaviour has been found also in time-dependent-driven quadratic systems [17,18]. Whereas in the mentioned cases the steady state is obtained as a fixed point of a nonunitary flow, generated by a Liouvillian of Lindblad form, a different approach to nonequilibrium statistical mechanics applicable to infinitely extended systems has been initiated by Ruelle [19][20][21] and has been actively developed since then [22]. A drawback of this approach is that the results are mainly limited to quasi-free systems. Thus, exact analytic results for steady states of interacting many-body models out of equilibrium are still scarce, and for non-integrable systems even scarcer.
In this paper we provide simple criteria (based on quantization in the space of operators) for a closed hierarchy of equations of motion for correlations of Markovian open quantum systems. The closeness property enables efficient calculation of the low-order steady-state correlations, although a complete description of the nonequilibrium steady state may still be difficult to assess. Interestingly, the closeness of correlation hierarchy can be exploited to analytically obtain exact steady-state correlations even in non-integrable interacting systems. We demonstrate this by calculating the two-point correlations in three paradigmatic interacting models, namely an interacting optomechanical model, the Fermi-Hubbard model, and the Rabi model. Solutions of these and similar out-of-equilibrium models can be important beyond the nonequilibrium setting, since they provide a new insight into the structure of interacting models.

Preliminaries
Before we proceed let us shortly introduce the method of quantization in the space of operators upon which the closeness criteria are based. We want to study the Lindblad master equation where {A, B} = AB + B A represents the anti-commutator, [A, B] = AB − B A the commutator, H is the Hamiltonian of the system and L k are the Lindblad operators. The first part on the righthand side of (1) generates the unitary evolution and the rest is responsible for dissipation. Our goal is to efficiently rewrite the Lindblad master equation (1) in the Liouville space by lifting the methods of second quantization to the space of operators. Since we want to describe fermions and bosons with the same formalism we have to introduce the space of trace class operators T (H) and the space of bounded 1 operators B(H) acting on the many-body Hilbert space H.
For each A ∈ B(H) we can define a functional (A| acting on |ρ ∈ T (H) as The functionals (A| form a dual space T (H) * [23]. We use the bra-ket notation and indicate the element of T (H) with |ρ and the element of T (H) * with (A|. The super-operators acting on the spaces T (H) from the left and T (H) * from the right shall be marked with a hat•. When expressing the action of the Liouvillian and dissipators on the ket vectors |ρ we shall omit the ket and write simply ρ, as in equation (1). In order to rewrite the Liouvillian, which determines the right-hand side of (1), in an efficient and clear way we introduce two fundamental superoperator maps, namely the leftL and the rightR multiplication maps defined aŝ Their action on the dual space can be deduced from the definition of the functionals (2). All bosonic (fermionic) Liouvillians of the form (1) can be rewritten with the aid of left L and rightR multiplication maps of creation a † j (c † j ) and annihilation a j (c j ) operators satisfying the canonical (anti-)commutation relations: . We use these sets of left and right multiplication maps to define creation and annihilation maps satisfying almost canonical commutation relations and almost canonical anti-commutation relations (almost-CAR) whereP =L (P) is the parity super-operator, is the parity operator and N is the number of bosonic or fermionic modes. The word almost indicates that the plus (+) inb + j orf + j indicates the bosonic or fermionic creation maps and not the adjoint maps (with respect to the inner product (2)) ofb j orf j , respectively. Equations (4) and (5) can be inverted and we can rewrite all Liouvillians of the Lindblad form in terms of creation and annihilation maps satisfying almost-CCR (almost-CAR) for bosons (fermions).
In the following, we shall use the super-operators in (4) and (5) to determine the closeness conditions for the hierarchy of equations for the correlation tensors in the bosonic, fermionic and mixed cases. In each case we separately consider the generators of unitary and nonunitary dynamics. We show that generators of the dissipative dynamics satisfying the closeness condition are not restricted to quadratic models, as in the unitary case. In the fermionic and bosonic cases we obtain general conditions for the closeness of the hierarchy for quadratic noise, whereas in the mixed case we show a simple example of quadratic noise with a closed hierarchy of equations. Finally, at the end of each section, we study three examples of interacting Hamiltonians, which by themselves do not exhibit a closed hierarchy, but in the presence of finetuned dissipation satisfy the closeness conditions. In particular, we find exact expressions for the time evolution of low-order correlations.

Bosonic Liouvillians
We begin with the bosonic case, which seems simpler due to canonical commutation relations.

Bosonic closeness condition
As discussed in the previous section we can always rewrite the bosonic LiouvillianB in terms of bosonic creation and annihilation mapŝ The non-negative integers m and n denote the number of creation and annihilation maps in the super-operatorsB (m,n) and the lengths of the vectors j and k, respectively. An important property of the creation mapsb + j defined in (4) is that they left-annihilate the identity operator, namely (1|b + j = 0. The trace preservation property (1|B = 0 and the condition (7) imply that in the sum (6) we always have m > 0. Let us now define a symmetric p-point correlation tensor (correlator) as The time evolution of correlator Z ( p) is determined by the equation Using the almost-CCR and the property (7) we observe that the expression where the minimum and maximum are taken over all paris (m, n) for whichB (m,n) in equation (6) is not zero. The property (7) implies that all correlators Z ( p−m+n) with p − m + n < 0 vanish. Therefore, the time evolution for the correlation tensors is closed iff M max 0 (bosonic closeness condition for the hierarchy of equations). For quadratic Hamiltonians and linear Lindblad operators we have m + n = 2 and since m > 0 the closeness condition is always satisfied. This has been exploited to find the properties of driven noninteracting bosonic systems [11] and simple harmonic chains coupled to heat baths [24]. Moreover, as we shall shortly show, the closeness condition for bosons can be satisfied also in a more general setting, e.g. harmonic chains in arbitrary dimension and in the presence of dephasing [16].

Unitary evolution
In this section we study as an example an interacting fourth-order Hamiltonian of the form The only term violating the bosonic closeness condition isB (1,3) . By using the symmetry of the tensor H with respect to permutation of the first and the last two indices, we find thatB (1,3) = 0 implies that H = 0. Not surprisingly, similar symmetry arguments can be applied for higherorder Hamiltonians showing that a hierarchy of equations for correlations cannot be closed for any non-quadratic (or linear) bosonic Hamiltonian (see appendix A).

Dissipation
In contrast to the unitary case, dissipators can contain a product of more than two creation and annihilation maps and nevertheless exhibit a closed hierarchy of equations for correlations. To illustrate this we consider an example of a purely dissipative evolution determined by a Lindblad operator L = j,k A jk d † j d k . The corresponding dissipative Liouvillian may be written aŝ The closeness condition M max = 0 implies that and is satisfied by quadratic Hermitian Lindblad operators. In fact, Hermitian A is the sole solution of equation (12), apart from an unimportant phase factor, showing that in the considered case quadratic Hermitian Lindblad operators are a necessary and sufficient condition for a closed hierarchy of correlation equations. In general we can discuss the closeness for a dissipator of the form where G is a positive Hermitian rate matrix, and µ, ν, µ , ν determine the lengths of vectors j, k, j , k , respectively (e.g. j = ( j 1 , j 2 , . . . , j µ )). In the bosonic case we may take the following coupling operators: In appendix B we show how the closeness condition imposes additional symmetry constraints on the rate matrix G. We find that the simplest additional requirements are satisfied by the rate matrix with the following symmetry: In the case of quadratic noise, i.e. restricting the set of coupling operators to L j,k ∈ {a j a k , a † j a † k , a † j a k }, the rate matrix G with the symmetries (15) has a closed hierarchy of equations. It remains an open question whether (15) is the only solution of the general conditions presented in appendix B.

Exact solution of a dissipative interacting Boson model
In this section we shall combine the results of the previous two sections and devote our attention to a combined unitary and dissipative evolution. The idea is to show that an interacting bosonic Hamiltonian can satisfy the closeness condition if coupled to fine-tuned Lindblad dissipators. For this purpose we study the simplest interacting bosonic Hamiltonian of the form which represents the interaction between one vibrational and one radiation mode and successfully describes most of the observed optomechanical phenomena [25]. The action of the adjoint mapâ d H VR can be expressed in terms of the bosonic creation and annihilation superoperators (4) aŝ Evidently, all but the last two interacting terms do not preserve the hierarchy of correlations. We wish to cancel them by adding dissipation of the form (13) with the coupling operator basis It can be shown that the rate matrix exactly cancels the unwanted terms in the adjoint map (18). For positive rates + 1 , + 2 , 1,2 2,1 , 1,1 > 0 and 4g 2 < min( + 1 1,2 , 1 2,1 , 2 1,1 ) the rate matrix (18) is positive and defines a valid Lindblad dissipatorD VR , which with the adjoint map (18) determines the complete LiouvillianB VR = −iâ d H VR +D VR . With the aid of creation and annihilation superoperators (4) we obtain Since m − n > 0 for all terms in the Liouvillian (19) the bosonic closeness condition is satisfied. This enables us to solve the differential equations (9) withB VR given by (19). In order to shorten the expressions we further restrict the rates + 1 = + 2 = 1,2 = 2,1 = + and 1 = 2 = 1,1 = . We assume that at initial time t = 0 the first and second modes are in a particle vacuum state, and find the following expressions for time evolution of non-vanishing first and second-order correlations: Since the remaining expressions for the time evolution of the occupations are too lengthy we write only their steady-state values All remaining one-and two-point correlations are zero. From the time evolution of the occupations (not shown) we see that the system is stable if > 3 + . Note that the hierarchy of equations for the considered model can be closed also by a different dissipator. In this regard it would be interesting to see how the properties of the systems change by changing the nonunitary part of the evolution.

Fermionic Liouvillians
In this section we define the fermionic closeness condition for unitary and dissipative dynamics and show by solving the dissipative Fermi-Hubbard model in one, two and three dimensions that fine-tuned dissipation can induce a closed hierarchy of equations also in the fermionic case.

Fermionic closeness conditions
As in the previous section we decompose the fermionic LiouvillianF in terms of creation and annihilation maps (5) The fermionic creation mapsf + j defined in (5) left-annihilate the identity operator, namely which again implies, by using the trace preservation property (1|F = 0, that in the sum (21) we always have m > 0. We define the anti-symmetric p-point correlation tensors as with α = 0, 1 determining the presence of the super-operatorP. An additional correlation tensor is necessary because of a parity super-operator in the fermionic Liouvillian (21). The time evolution of correlation tensors Due to the phase factor in the Liouvillian and two different types of correlations it is more difficult to find the closeness conditions than in the bosonic case. The structure of the derivative (24) is best explained by the following diagrams: The diagrams (25) show how the time derivatives of correlation tensors Z ( p,α) transform under the action of the super-operatorsF (m,n) andF (m,n)P appearing in the fermionic Liouvillian (21). The left and the right diagram have a different sign of m and n on the right-hand side of the diagrams. This is a consequence of an additional phase super-operatorP in Z p,1 , which satisfies the following relations: Hence, the super-operatorsF (m,n) , which lower the order of Z ( p,0) , increase the order of Z ( p, 1) . Observing this we find three fermionic closeness conditions for the time evolution of the correlation tensors Z ( p,α) : • m + n : even and n − m 0 for all non-vanishingF (m,n) in (21) ; α = 0; • m + n : even and m − n 0 for all non-vanishingF (m,n) in (21) ; α = 1; • (m − n = K , where K is an odd integer for all non-vanishingF (m,n)P in (21) with odd m + n) or (m − n = 0 for all non-vanishingF (m,n) in (21) with even m + n); α = 0, 1.
With α we mark which correlations Z ( p,α) exhibit a closed hierarchy. The first condition has been exploited to find steady states of a dissipative one-dimensional non-interacting spin system [10,[26][27][28], the XX model with dephasing [15], and the quantum analogue of a simple symmetric exclusion process [14].
Again we first independently consider the closeness conditions for unitary and non-unitary generators. In case of unitary evolution we find that the closeness condition cannot be satisfied with non-quadratic (or linear) Hamiltonians (see appendix A). In the case of a dissipative evolution defined with the dissipator (13) and the fermionic coupling operator basis the closeness condition implies additional symmetries of the matrix G. The procedure to attain these symmetries is outlined in appendix B. We explicitly provide only the 'lowest-order' constraints, which are satisfied by the rate matrix with the symmetries where we define the phase factor (µ, ν, µ ν ) = (−1) µ+ν 2 + µ +ν 2 +1+νν +µ ν+µ ν +µµ +µν+µν . In the simplest nontrivial example with the coupling operator basis L j,k ∈ {c † j c k , c † j c † k , c j c k } the symmetry (28) is also sufficient for the first fermionic closeness condition to be satisfied. Conditions (28) extend the result of [29], where it was shown that Hermitian Lindblad operators of the form L = j,k A jk c † j c k are sufficient for a closed hierarchy of correlations. In fact by using similar arguments as in the bosonic case for the dissipator (11) we can prove that in this simple case Hermitian A is also a necessary condition for the closeness of the hierarchy. However, it is not clear whether (28) is the only solution to the conditions for the general quadratic noise presented in appendix B.

Exact solution of a dissipative Fermi-Hubbard model
We argued that interacting fermionic models do not satisfy the fermionic closeness conditions whereas the 'interacting' fermionic dissipation does. Here we shall show that adding controlled dissipation can lead to a closed hierarchy of equations for the Fermi-Hubbard model with the Hamiltonian The Fermi-Hubbard model introduced in [30][31][32] is the simplest model describing basic concepts of condensed matter physics, such as superconductivity and magnetic and electronic properties of materials and can moreover be realized with cold Fermi gases in optical lattices [33]. Although the one-dimensional Fermi-Hubbard model is Bethe ansatz solvable [34] we shall abstain from using the quantum integrability insight into the structure of the model, but will instead employ a controlled dissipation to simplify the treatment of the model in arbitrary dimensions, i.e. for arbitrary hopping t j,k and on-site interaction U j ( j, k = 1, 2, . . . , N ). Since the only term in (29) violating the fermionic closeness conditions is the interacting term, we can without loss of generality focus on the Hamiltonian of the form For brevity we omit the sub-index j and define c ( †) 1,2 := c ( †) ↑,↓ . First, we write the adjoint map of (30) with aid of the creation and annihilation super-operatorŝ where the rest contains the hierarchy-preserving terms. Next, we wish to add dissipation of the form (13) with the coupling operator basis L j, which cancels the unwanted terms in (31). In appendix C we provide a positive Hermitian rate matrix for which the dissipator (13) exactly cancels the hierarchy violating terms such that the complete Liouvillian satisfies the fermionic closeness condition. For this reason it is possible to calculate the dynamics of low order correlations as well as their steady state values for a general Fermi-Hubbard model (29) with a fine-tuned dissipation (see appendix C). In particular, we find in one dimension with nearest-neighbour hopping that all steady-state twopoint correlations vanish aside from the occupation numbers c † j,σ c j,σ = G 55 G 11 −G 33 +G 55 . In two and three dimensions we were unable to find an analytic expression for the covariance matrix. However, we numerically find exponential decay of correlations with the distance from the diagonal in the two-dimensional (2D) case (see appendix C). In three dimensions we observe a revival of correlations at the edges of the cubic lattice, as can be seen in figure 1(a) where we show the correlations with respect to the site in the middle of the lattice. We also show in figure 1(a) the spatial distribution of the occupations and in figure 1(b) the non-vanishing part of the full two-point correlation matrix, which indicates exponential decay of correlations with the distance between sites. Interestingly, in the described cases the steady state expectation values of the two-point correlations do not depend on the interaction strength. However, their evolution before reaching the steady state and the higher order correlations still do. A different dissipation can induce interaction-dependent two-point steady-state correlations, and can perhaps reveal some interesting properties of the Fermi-Hubbard model in one or more dimensions.

Mixed Liouvillians
In this section we determine the mixed closeness conditions. We show that in contrast to the unitary evolution the purely dissipative evolution can exhibit a closed hierarchy. Finally, we show an example of dissipation-induced hierarchy of equations for the mixed correlations by providing exact analytic results for time evolution of second-order correlations in the dissipative Rabi model.

Mixed closeness conditions
We anew employ the same strategy as for bosons and fermions and rewrite the complete mixed Liouvillian withB (m ,n ) andF (m,n) as defined in (6) and (21) with α = 0, 1 determining the presence of the parity super-operator in the correlation matrix.
Resembling the fermionic case the two different correlators are necessary due to the presence of the fermionic parity super-operator in the mixed Liouvillian. The time evolution of the correlators (33) is governed by The action of the mixed Liouvillian (32) on the correlation tensors can be explained by the following diagrams: If the only non-vanishingF (m,n) in (32) has m + n = 0, above conditions reduce to simple bosonic closeness conditions and similarly if the only non-vanishingB (m ,n ) in (32) has m + n = 0 they simplify to fermionic closeness conditions. Once more we separately consider unitary and dissipative dynamics and observe that the mixed closeness conditions are not satisfied for any fermion-boson Hamiltonian (see appendix A). In contrast, we can find a fermion-boson dissipator satisfying one of the mixed closeness conditions. The simplest example considered in appendix B is a dissipative evolution determined by (13), with the coupling operator basis L ∈ {c, c † , a, a † } and the rate matrix The positivity of G implies |G 13 | < G 11 G 44 and G 11 , G 44 > 0, whereas the stability of the system implies δ > 0.

Exact solution of a dissipative Rabi model
In previous sections we showed that in contrast to unitary evolution purely dissipative evolution satisfies the closeness condition also for non-quasi-free higher-order Liouvillians, and can moreover induce a closed hierarchy in the case of interacting fermionic and bosonic Hamiltonians. Now we extend this result and show that controlled coupling to a Markovian environment can induce a closed hierarchy of equations also in the case of interacting mixed Hamiltonians. For this purpose, we study the Rabi model [35], which describes the interaction between a two-level system (spineless fermion) and a quantized harmonic oscillator (boson), and is given by the Hamiltonian Despite its simplicity it displays a rich behaviour and has a wide range of applicability in quantum optics [36], quantum information [37] and condensed matter [38]. Moreover, the time evolution of interesting observables is still restricted to numerical [39] and analytical [40] approximations. Here we show how the model can be simplified by adding Markovian dissipation which leads to a closed hierarchy of equations for the correlation tensors. We provide exact closed-form expressions for the time evolution of low-order correlations.
Let us start by writing the action of the adjoint map of the Hamiltonian in terms of the fermionic and bosonic creation and annihilation super-operatorŝ We observe that the closeness conditions are not satisfied due to the terms of the form f + jb kP , where j, k = 1, 2. The idea is to add noise (dissipation) which exactly cancels these terms. Again we employ the dissipator of the form (13) with the coupling operator basis L ∈ {c, c † , a, a † } and rate matrix This dissipator exactly cancels the unwanted terms in (39). Moreover, the so far unspecified rates G i, j in (40) can be chosen such that G is a positive Hermitian matrix and consequently defines a valid Lindblad dissipator. For simplicity we choose G 11 Rabi

expressed in terms of bosonic (4) and fermionic (5) creation and annihilation maps iŝ
Since the Liouvillian (41) does not contain any bosonic annihilation maps we have a closed hierarchy of equations for the correlation tensors, which are easily solvable for low order correlations. In particular, we find the following expressions for correlations of the first and second order: Since the remaining expressions for the time evolution of correlations are too lengthy we provide only their steady-state values.
2δ , The system is initially (at time t = 0) in bosonic vacuum and fermionic totally mixed state. The system is stable if δ > 0, i.e. when more bosonic modes are incoherently dissipated than pumped into the system. We may include also coherent pumping or more fermions or bosons and treat the model in a similar manner. In all three considered examples, i.e. the interacting bosonic Hamiltonian, Fermi-Hubbard model and Rabi model, the dissipators used to simplify the hierarchy of correlations are represented in the first standard form; hence it is difficult to comprehend which physical processes they represent. However, even if the rate matrix is diagonalized and the corresponding diagonal form is obtained the physical interpretation of all processes can be difficult. We shall illustrate this on the simplest of the considered examples, the Rabi model. Its rate matrix (40) can easily be diagonalized. The eigenvectors and the square roots of the corresponding eigenvalues determine up to an unimportant complex phase factor the following Lindblad operators: The first Lindblad operator describes a bit-phase flip channel with the rate flip = 2 . The second Lindblad operator describes an amplitude damping channel of the bosonic mode with the rate amp = + 2δ. Interpretation of the remaining channels is not so straightforward since they represent a collective decoherence of a fermionic and bosonic mode with rates c 1 = 1 2 + √ 2g and c 2 = 1 2 − √ 2g. Nevertheless, these channels could, in principle, be physically realized for example following ideas along the lines of [41] and using a stroboscopic unitary time evolution of the Rabi model coupled to a bosonic mode initially prepared in a vacuum state. This transforms the problem of channel engineering to the problem of constructing the appropriate coupling between the environment bosonic mode and the degrees of freedom of the Rabi model 2 . Similar collective channels are expected to appear in all cases of dissipation-induced closed hierarchy of correlations.
The rates of the bit-phase flip channel and the bosonic amplitude damping channel can in fact be larger without altering the closeness property of the model, namely flip > 2 and amp > + 2δ, respectively. Hence, only the rates corresponding to collective channels have to be fine-tuned, i.e. c 1 − c 2 = 2 √ 2g and c 1 , c 2 > 0.

Conclusions
We determined simple criteria for the closeness of the hierarchy of equations in Markovian open quantum systems. In particular, we discussed general quadratic fermionic and bosonic noise. We showed that dissipative systems satisfy the closeness property under more general conditions than closed systems. In the case of unitary evolution the closeness of the hierarchy is possible only in quadratic systems, in other words for an evolution which preserves Gaussian states. In the dissipative case, however, the closeness conditions can be satisfied by more general systems, of which steady states are not necessarily Gaussian states. Nevertheless, their low order correlations (as well as their dynamics before reaching the steady state) can be efficiently computed. Moreover, we considered three examples of dissipation-induced closeness of the hierarchy, where the unitary evolution, determined by interacting Hamiltonians, by itself does not satisfy the closeness condition, whereas the joint non-unitary evolution does. This enabled us to find exact expressions for the time evolution of the first-and second-order correlations of paradigmatic interacting quantum models, which are otherwise accessible only through numerical approximations. Although the solutions may seem rather artificial due to an additional fine-tuned noise they can perhaps provide new insight into the structure of interacting quantum models in one and more dimensions and can be extended by means of perturbation theory in the interaction strength.
Albeit the closeness conditions enable simple calculations of low-order correlations, a general procedure how to evaluate higher-order correlations, as Wick's theorem for Gaussian systems, remains an interesting open problem. Since the closeness condition also enables the calculation of some eigenvalues of the Liouvillian it would be interesting to see whether these eigenvalues include also the spectral gap, which is the eigenvalue with the largest real part and generically governs the convergence of the evolution to the steady state.
where the indices µ and ν determine the lengths of the vectors j and k, respectively. It is clear that due to the trace preservation and the property 1|b + = 0 linear and quadratic Hamiltonians automatically satisfy the closeness condition. In order to show that higher-order Hamiltonians do not respect this rule we look at terms of the adjoint mapâ d H ρ = [H, ρ] that contain only one bosonic creation mapb + j p orb + The rest contains all the remaining terms, which have two or more bosonic creation maps or in which the number of creation and annihilation maps is not equal to µ + ν. One can immediately see that due to the symmetry H P j,P k = H j,k , where P and P represent arbitrary permutations of elements of vectors j and k, respectively, the explicit sums in equation (A.2) cannot vanish. Therefore, the mapâ d H always contains terms with only one creation superoperator and ν + µ − 1 annihilation super-operators, which leads to the following closeness condition: M max = max(ν + µ − 2) 0, where the maximum is taken over all pairs µ, ν in (A.1) and (A.2) for which H (µ,ν) does not vanish. In other words, the closeness condition for bosonic Hamiltonians is satisfied only for quadratic and linear bosonic Hamiltonians.

A.2. Fermionic Hamiltonians
In this subsection we discuss the closeness conditions for general fermionic Hamiltonians Explicitly writing only terms with exactly one fermionic creation mapf + j p orf + N + j p the adjoint map of the fermionic Hamiltonian (A.3) readŝ The rest contains all terms which either have two or more creation maps or consist of less than µ + ν creation and annihilation maps. By changing the order of operators in the second and fourth sum in (A.4) we obtain (A.5) Using the symmetry H P j,P k = ε P j ε P k H j,k , with P and P being arbitrary permutations of elements of vectors j and k, respectively, and ε j a completely antisymmetric tensor, we observe that terms with one creation map and µ + ν − 1 annihilation maps do not vanish. Hence, for even µ + ν the first fermionic closeness condition reduces to ν + µ − 2 0. In order to show that the second and the third condition cannot be satisfied for µ + ν > 2, we examine other terms with a maximal number of creation and annihilation maps. Let us first write out only the terms consisting only of creation mapŝ (A.6) where the rest contains all nonvanishing terms with at least one annihilation map. For odd µ + ν explicit terms in (A.6) do not vanish. Due to (A.5) and (A.6) the fermionic adjoint map (21) with odd µ + ν always contains terms with n − m < 0 and terms with m > 1 and n = 0 and cannot satisfy the third fermionic closeness condition unless µ + ν = 1. Since for even µ + ν the explicit terms in (A.6) vanish we have to consider terms of the adjoint map containing only one annihilation map and µ + ν − 1 creation maps. After a tedious calculation we find the following simplified expression +rest, implying that the second fermionic condition can not be satisfied unless µ + ν = 2. We showed that similarly as in the bosonic case the only two fermionic Hamiltonians with the closed hierarchy of equations for correlations are quadratic and linear Hamiltonians. Since explicit terms in (A.6) vanish the Hamiltonian with quadratic and linear terms also satisfies the third fermionic closeness condition.

A.3. Mixed Hamiltonians
In this subsection we determine the closeness condition for the mixed Hamiltonian of the form with µ + ν 1 and µ + ν 1. The fermionic and the bosonic maps commute. Accordingly, the adjoint map of (A.7) contains the same fermionic terms as in (A.5), (A.6) and (A.7) multiplied by µ + ν bosonic annihilation maps. Therefore, Hamiltonians containing fermion-boson interaction cannot satisfy the closeness conditions.

Appendix B. Closed hierarchy for Lindblad dissipators
In this appendix we discuss Lindblad dissipators of the form Dρ = µ,ν,µ ,ν j,k; j ,k G (µ,ν;µ ,ν ) where G is a positive Hermitian matrix rate-matrix, µ, ν, µ ν designate the lengths of multiindices j, k, j , k (e.g. j = ( j 1 , j 2 . . . , j µ ) and L jk designates the coupling operators, which shall be defined in the following sections. Our aim is to find necessary and sufficient conditions for the closed hierarchy of equations for correlations (8), (23) and (34). We separately study the bosonic, fermionic and mixed cases. For quadratic bosonic and fermionic Lindblad operators we determine a set of symmetry conditions for G which is equivalent to the closeness condition.
In the mixed case we study the simplest nontrivial example and show that in contrast to the unitary case the closeness conditions can be satisfied.
These equations have to be satisfied for all values of j 1 , j 2 , t 1 , t 2 , k 1 , k 2 and determine the most general set of quadratic Markovian bosonic noise which satisfies the bosonic closeness condition.
The 'higher-order' requirements can be obtained in a similar manner by writing the vanishing conditions for unwanted terms in the dissipators.

B.3. Mixed dissipator -Example
In this appendix we show an example of a mixed dissipator satisfying the third mixed closeness condition. We consider a dissipative evolution determined by the dissipator (13) with the coupling operator basis L j ∈ {c, c † , a, a † }. Since our purpose is only to give an example of a Lindblad dissipator satisfying the mixed closeness conditions, we restrict ourselves to the following rate matrix: and satisfies the third mixed closeness condition. For the system to be stable δ should be positive; in other words the bosonic mode should be more dissipated than incoherently pumped into the system.   This form is suitable for the efficient solutions of (24). We calculate the steady-state twopoint correlations in one, two, and three dimensions. We assume nearest-neighbour hopping.
In one dimension all two-point correlations vanish except for the occupations which are given by c † j,σ c j,σ = G 55 G 11 −G 33 +G 55 and are independent of the length of the chain and disorder in the hopping. In 2D the two-point correlations c † j,σ c k,σ do not vanish and decay exponentially with the distance from the diagonal | j − k| (see figure C.1); all other two-point steady state correlations still vanish. Similar indications of exponential decay of correlations are observed in three dimensions, although we also observe small correlation revivals at the edges of the chain (see main text figure 1).