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 necessary 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 opportunity to compare numerical simulations and 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 non-interacting harmonic oscillator lattices in arbitrary dimensions with dephasing [16]. Exact long-time behavior 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 non-unitary 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 [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 manybody 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 loworder 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 twopoint 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 + BA represents the anti-commutator, [A, B] = AB − BA the commutator, H is the Hamiltonian of the system, and L k are the Lindblad operators. The first part on the right hand 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 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 superoperators 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 super-operator 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 leftL and rightR multiplication maps of creation a † j (c † j ) and annihilation a j (c j ) operators satisfying the canonical (anti-)commutation relations; [a j , a . We use these sets of left and right multiplication maps to define creation and annihilation maps satisfying almost canonical commutation relations (almost-CCR) 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. The 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 case. In each case we separately consider the generators of unitary and non-unitary 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 case 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 fine-tuned 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 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 (1|b j 1b j 2 . . .b jpB |ρ(t) is a function of correlators Z (p+M min ) , Z (p+M min +1) , . . . Z (p+Mmax) only, with M min = min(n − m) and M max = max(n − m), 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 H = N i,j,k,l=1 H ijkl a † i a † j a k a l . Due to Hermicity of the Hamiltonian we have H i,j,k,l =H l,k,j,i , where the bar• denotes complex conjugation. By using (4) we write the adjoint map of the Hamiltonianâd H := [H, •] in the form (6) 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 H = 0. Not surprisingly, similar symmetry arguments can be applied for a higher order Hamiltonians showing that a hierarchy of equations for correlations cannot be closed for any non-quadratic (or linear) bosonic Hamiltonian (see the 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ŝ 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 necessary and sufficient condition for closed hierarchy of correlation equations. In in general we can discuss the closeness for a dissipator of the form where G is a positive Hermitian rate matrix, µ, ν, µ , ν 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 the Appendix B we show how the closeness condition imposes additional symmetry constraints to the rate matrix G. We find that the simplest additional requirements are satisfied by rate matrix with the following symmetry In 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 if (15) is the only solution of the general conditions presented in the 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 finetuned 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 super-operators (4) aŝ 3 . 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 exactly cancels the unwanted terms in the adjoint map (17). For positive rates (18) is positive and defines a valid Lindblad dissipatorD VR , which with the adjoint map (17) determines the complete LiouvillianB VR = −iâd H VR +D VR . With the aid of creation and annihilation super-operators (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 mode 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Γ + . Notice 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 non-unitary 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 relationŝ 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 dissipative one-dimensional noninteracting spin system [10,26,27,28], the XX model with dephasing [15], and the quantum analog of simple symmetric exclusion process [14].
Again we first independently consider the closeness conditions for unitary and nonunitary generators. In case of unitary evolution we find that the closeness condition cannot be satisfied with non-quadratic (or linear) Hamiltonians (see the Appendix A). In 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 the 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) +1+νν +µ ν+µ ν +µµ +µν+µν . In the simplest nontrivial example with the coupling operator basis (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 a sufficient for 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 if (28) is the only solution to the conditions for the general quadratic noise presented in the 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 closed hierarchy of equation 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, as superconductivity and magnetic end 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 dimension, i.e. for arbitrary hopping t j,k and onsite 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 ( †) ↑,↓ . 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 which cancels the unwanted terms in (31). In the 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 the Appendix C). In particular we find in one-dimension with nearest neighbour hopping that all steady state two-point 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 case (see the 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) 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 the Appendix A). On the contrary, we can find a fermion-boson dissipator satisfying one of the mixed closeness conditions. The simplest example considered in the 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 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 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 behavior 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 first 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 + jbkP , 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 = G 22 = G 44 = Γ > 2 √ 2|g|, G 33 = Γ + 2δ > 0, and G 2,1 = G 1,2 = G 3,4 = G 4,3 = 0. In this case the eigenvalues of G are positive (Γ, Γ + 2δ, Γ − 2 √ 2 g, Γ + 2 √ 2 g). The full LiouvillianL Rabi = −iâd H Rabi +D 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. .
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.

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 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 necessary 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 the 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 if this 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.
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 + jp orb + N +jp The rest contains all 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 super-operator 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.

Appendix 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 + jp orf + N +jp 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 an fourth sum in (A.4) we obtain 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 can not be satisfied for µ + ν > 2, we examine other terms with maximal number of creation and annihilation maps. Let us first write out only the terms consisting only of creation mapŝ where the rest contains all non-vanishing 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 can not 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 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.

Appendix B. Closed hierarchy for Lindblad dissipators
In this Appendix we discuss Lindblad dissipators of the form where G is a positive Hermitian matrix rate-matrix, µ, ν, µ ν designate the lengths of multi-indices 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 bosonic, fermionic and mixed case. 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.

Appendix B.1. Bosons
In general the bosonic coupling operators can be written as We take the same approach as in the unitary case. We first write the dissipator by using creation and annihilation maps and then consider terms with different number of creation maps. From this expansion we extract independent, necessary requirements for the closeness condition. Let us first focus on terms containing only one creation map and m + n + m + n − 1 annihilation mapŝ Since the tensor G is symmetric with respect to permutations of elements of each multi index j, k, t, v we obtain the conditions Now we introduce the vectors x = (j, v) and y = (t 2 , t 3 . . . , t m , k 1 , k 2 , . . . k n ) and extract from the equation (B.4) the following symmetry conditions for the rate matrix G which have to be satisfied for all variations of x, y and t 1 . The simplest solution of the equations (B.5) is if all terms in the sum vanish separately, namely if G . The question remains if this is in general the only solution which satisfies also the positivity constraint. Equations (B.5) are necessary condition for the closeness of the hierarchy. However, they are sufficient only int the case where µ + ν + µ + ν = 4. In general we obtain in a similar manner also "higher-order" necessary conditions for the closeness of the hierarchy and it is not clear if they can always be satisfied. Nevertheless, for a general quadratic noise determined by coupling operators from the set {a j a k , a † j a † k , a † j a k } equations (B.5) are sufficient for the closeness of the hierarchy and simplify to 0 = 4(G This equations have to bee 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 obtain in a similar manner, by writing the vanishing conditions for unwanted terms in the dissipators.

Appendix 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 which for |G 24 | 2 < |G 22 G 44 | and |G 13 | 2 < |G 11 G 33 | defines a valid dissipator. Rewriting the dissipator in terms of the fermionic and bosonic creation and annihilation maps it is not difficult to see that the mixed closeness conditions are satisfied if G 11 = G 22 and G 13 = G 24 . For simplicity we write G 33 = G 44 + 2δ. In this case the dissipator may be written asD This form is suitable for the efficient solutions of (24). We calculate the steady-state two-point 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 two dimensions the two-point correlations c † j,σ c k,σ do not vanish and decay exponentially with the distance from the diagonal |j − k| (see Figure C1); 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