Closed hierarchy of correlations in Markovian open quantum systems

As discussed in the previous section we can always rewrite the bosonic Liouvillian $\skew3\hat{\mathcal {B}}$ in terms of bosonic creation and annihilation maps

$$\begin{eqnarray} &&\begin{array}{ll} \displaystyle\skew3\hat{\mathcal{B}}=\sum_{m,n}\skew3\hat{\mathcal{B}}^{(m,n)},\\ \displaystyle\skew3\hat{\mathcal{B}}^{(m,n)}=\sum_{j_1\ldots j_m,k_1\ldots k_n}B^{(m,n)}_{j_1\ldots j_m,k_1\ldots k_n}\hat{b}_{j_1}^+\ldots\hat{b}_{j_m}^+\hat{b}_{k_1}\ldots\hat{b}_{k_n}. \end{array} \end{eqnarray} \tag{ 6 }$$

The non-negative integers m and n denote the number of creation and annihilation maps in the super-operators $\skew3\hat{\mathcal {B}}^{(m,n)}$ and the lengths of the vectors $\underline {j}$ and $\underline {k}$, respectively. An important property of the creation maps $\hat{b}_{j}^+$ defined in (4) is that they left-annihilate the identity operator, namely

$$\begin{equation} \left (\mathbbm{1}\right | \hat{b}^+_j=0. \end{equation} \tag{ 7 }$$

The trace preservation property $\left ( \mathbbm{1}\right | \skew3\hat{\mathcal {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

$$\begin{equation} Z^{(p)}_{j_1,j_2,\ldots,j_p}(t)=\left ( 1\right | \hat{b}_{j_1}\hat{b}_{j_2}\ldots\hat{b}_{j_p}\left | \rho(t)\right \rangle. \end{equation} \tag{ 8 }$$

The time evolution of correlator Z^(p) is determined by the equation

$$\begin{equation} \frac{{\rm d}}{{\rm d} t} Z^{(p)}_{j_1,j_2,\ldots,j_p}(t)=\left ( 1\right | \hat{b}_{j_1}\hat{b}_{j_2}\ldots\hat{b}_{j_p}\skew3\hat{\mathcal{B}}\left | \rho(t)\right \rangle. \end{equation} \tag{ 9 }$$

Using the almost-CCR and the property (7) we observe that the expression $\left ( 1\right | \hat{b}_{j_1}\hat{b}_{j_2}\ldots \hat{b}_{j_p}\skew3\hat{\mathcal {B}}\left | \rho (t)\right \rangle$ is a function of correlators ${\bf{ Z}}^{(p+M_{\min })},{\bf{ Z}}^{(p+M_{\min }+1)},\ldots {\bf{ Z}}^{(p+M_{\max })}$ 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 which $\skew3\hat{\mathcal {B}}^{(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 }\leqslant 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].

In this section we study as an example an interacting fourth-order Hamiltonian of the form $H=\sum _{i,j,k,l=1}^NH_{ijkl}a_i^\dagger a_j^\dagger a_ka_l$. Due to hermicity of the Hamiltonian we have $H_{i,j,k,l}=\bar {H}_{l,k,j,i}$, where the bar $\bar {\bullet }$ denotes complex conjugation. By using (4) we write the adjoint map of the Hamiltonian $\skew3\hat{{\rm ad}}_{H}:=[H,\bullet ]$ in the form (6)

$$\begin{eqnarray} &&\begin{array}{ll} \displaystyle\skew3\hat{\mathcal{B}}_{\rm int}=-{\rm i}\,\skew3\hat{\rm ad}_{H}=\skew3\hat{\mathcal{B}}_{\rm int}^{(1,3)}+\sum_{m\geqslant n}\skew3\hat{\mathcal{B}}_{\rm int}^{(m,n)}, \\ \displaystyle\skew3\hat{\mathcal{B}}_{\rm int}^{(1,3)}=-{\rm i}\sum_{i,j,k,l}\left(H_{ijkl}\, \hat{b}_{i}^+ \hat{b}_{N+j} \hat{b}_{k} \hat{b}_{l} + \,\hat{b}_{N+k}^+ \hat{b}_{N+i} \hat{b}_{N+j} \hat{b}_{l}\right). \end{array} \end{eqnarray} \tag{ 10 }$$

The only term violating the bosonic closeness condition is $\skew3\hat{\mathcal {B}}^{(1,3)}$. By using the symmetry of the tensor H with respect to permutation of the first and the last two indices, we find that $\skew3\hat{\mathcal {B}}^{(1,3)}=0$ implies that H = 0. Not surprisingly, similar symmetry arguments can be applied for higher-order Hamiltonians showing that a hierarchy of equations for correlations cannot be closed for any non-quadratic (or linear) bosonic Hamiltonian (see appendix A).

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=\sum _{j,k}A_{jk}d_j^\dagger d_k$. The corresponding dissipative Liouvillian may be written as

$$\begin{eqnarray} &&\fl \skew3\hat{\mathcal{B}}_{\rm dis}=\skew3\hat{\mathcal{B}}_{\rm dis}^{(1,3)}+\sum_{m\geqslant n}\skew3\hat{\mathcal{B}}_{\rm dis}^{(m,n)}, \nonumber\\ &&\fl \skew3\hat{\mathcal{B}}_{\rm dis}^{(1,3)}=\frac{1}{2}\sum_{i,j,k,l}((A_{il} \bar{A}_{kj}-A_{j,l} \bar{A}_{k,i}) \hat{b}_{i}^\dagger \hat{b}_{N+j} \hat{b}_{k} \hat{b}_{l}+ (A_{ki} \bar{A}_{lj}-A_{kl} \bar{A}_{ij}) \hat{b}_{N+i}^\dagger \hat{b}_{N+j} \hat{b}_{N+k} \hat{b}_{l}). \end{eqnarray} \tag{ 11 }$$

The closeness condition $M_{\max }=0$ implies that

$$\begin{equation} A_{il} \bar{A}_{kj}-A_{j,l} \bar{A}_{k,i}+A_{ik} \bar{A}_{lj}-A_{j,k} \bar{A}_{l,i}=0, \end{equation} \tag{ 12 }$$

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

$$\begin{eqnarray} &&\begin{array}{ll} \displaystyle \skew3\hat{\mathcal{D}}=\sum_{\mu,\nu,\mu',\nu'}\sum_{\underline{j},\underline{k},\underline{j}',\underline{k}'}G^{(\mu,\nu;\mu',\nu')}_{\underline{j},\underline{k};\underline{j}',\underline{k}'}\skew3\hat{\mathcal{D}}(L_{\underline{j},\underline{k}}, L_{\underline{j}',\underline{k}'}),\\ \displaystyle \skew3\hat{\mathcal{D}}(L_{\underline{j},\underline{k}}, L_{\underline{j}',\underline{k}'})\rho=L_{\underline{j},\underline{k}}\rho L_{\underline{j}',\underline{k}'}^\dagger-\frac{1}{2}\{L_{\underline{j}',\underline{k}'}^{\dagger} L_{\underline{j},\underline{k}},\rho\}, \end{array} \end{eqnarray} \tag{ 13 }$$

where G is a positive Hermitian rate matrix, and μ, ν, μ', ν' determine the lengths of vectors $\underline {j},\,\underline {k},\,\underline {j}',\,\underline {k}'$, respectively (e.g. $\underline {j}=(j_1,j_2,\ldots ,j_\mu )$). In the bosonic case we may take the following coupling operators:

$$\begin{equation} L_{\underline{j},\underline{k}}=a^\dagger_{j_1}a^\dagger_{j_2}\ldots a^\dagger_{j_\mu} a_{k_1}a_{k_2}\ldots a_{k_\nu}. \end{equation} \tag{ 14 }$$

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:

$$\begin{equation} G^{(\mu,\nu;\mu',\nu')}_{\underline{j},\underline{k};\underline{j}',\underline{k'}}-\bar{G}^{(\nu,\mu;\nu',\mu')}_{\underline{k},\underline{j};\underline{k}',\underline{j}'}=0. \end{equation} \tag{ 15 }$$

In the case of quadratic noise, i.e. restricting the set of coupling operators to L_j,k∈{a_ja_k,a^†_ja^†_k,a^†_ja_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.

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

$$\begin{equation} H_{\rm VR}=\omega \,a_1^\dagger a_1+\Omega \,a_2^\dagger a_2 + g \,a_1^\dagger a_1(a_2+a_2^\dagger), \end{equation} \tag{ 16 }$$

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 $\skew3\hat{\mathrm { ad}}_{H_{\mathrm { VR}}}$ can be expressed in terms of the bosonic creation and annihilation super-operators (4) as

$$\begin{eqnarray} &&\fl \skew3\hat{\mathrm{ad}}_{{H}_{\rm VR}}=\omega(\hat{b}_1^+\hat{b}_1-\hat{b}_3^+\hat{b}_3)+\Omega(\hat{b}_2^+\hat{b}_2-\hat{b}_4^+\hat{b}_4)+g\Big(\hat{b}^+_1 \hat{b}_2 \hat{b}_1 + \hat{b}^+_1 \hat{b}_4 \hat{b}_1 + \hat{b}^+_2 \hat{b}_3 \hat{b}_1 - \hat{b}^+_3 \hat{b}_3 \hat{b}_2 \nonumber\\ &&- \hat{b}^+_3 \hat{b}_4 \hat{b}_3 - \hat{b}^+_4 \hat{b}_3 \hat{b}_1 + \hat{b}^+_2 \hat{b}^+_1 \hat{b}_1 - \hat{b}^+_4 \hat{b}^+_3 \hat{b}_3\Big). \end{eqnarray} \tag{ 17 }$$

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 $L_{\underline {j},\underline {k}}\in \{ a_1^\dagger ,a_2^\dagger ,a_1,a_2,a^\dagger _1a_2^\dagger ,a^\dagger _2a_1,a^\dagger _1a_1,a_2a_1,a^\dagger _1a_2 \}$. It can be shown that the rate matrix

$$\begin{eqnarray} {\bf G}=\left(\begin{array}{@{}ccccccccc@{}} \Gamma_1^+ & 0 & 0 & 0 & -2 {\rm i} g & 0 & 0 & 0 & 0\\ 0 & \Gamma_2^+ & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & \Gamma_1 & 0 & 0 & -2 {\rm i} g & 0 & 0 & 0\\ 0 & 0 & 0 & \Gamma_2 & 0 & 0 & -2 {\rm i} g & 0 & 0\\ 2{\rm i} g & 0 & 0 & 0 & \Gamma_{1,2} & 0 & 0 & 0 & 0\\ 0 & 0 & 2{\rm i} g & 0 & 0 & \Gamma_{2,1} & 0 & 0 & 0\\ 0 & 0 & 0 & 2{\rm i} g & 0 & 0 & \Gamma_{1,1} & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \Gamma_{1,2} & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \Gamma_{2,1} \end{array}\right) \end{eqnarray} \tag{ 18 }$$

exactly cancels the unwanted terms in the adjoint map (18). For positive rates Γ⁺₁,Γ⁺₂,Γ_1,2Γ_2,1,Γ_1,1 > 0 and $4g^2<\min (\Gamma _1^+\Gamma _{1,2},\,\Gamma _1\Gamma _{2,1},\,\Gamma _2\Gamma _{1,1})$ the rate matrix (18) is positive and defines a valid Lindblad dissipator $\skew3\hat{\mathcal {D}}_{\mathrm { VR}}$, which with the adjoint map (18) determines the complete Liouvillian $\skew3\hat{\mathcal {B}}_{\mathrm { VR}}=-{\mathrm { i}}\,\skew3\hat{\mathrm { ad}}_{H_{\mathrm { VR}}}+\skew3\hat{\mathcal {D}}_{\mathrm { VR}}$. With the aid of creation and annihilation super-operators (4) we obtain

$$\begin{eqnarray} &&\fl \skew3\hat{\mathcal{B}}_{\rm VR}=-{\rm i} \omega(\hat{b}_1^+\hat{b}_1-\hat{b}_3^+\hat{b}_3)-{\rm i} \Omega(\hat{b}_2^+\hat{b}_2-\hat{b}_4^+\hat{b}_4)+\,{\rm i}\, g (\hat{b}^+_2 - \hat{b}^+_4 - 2 (\hat{b}^+_2 \hat{b}^+_1 \hat{b}_1 + \hat{b}^+_3 \hat{b}^+_1 \hat{b}_2 - \hat{b}^+_3 \hat{b}^+_1 \hat{b}_4 \nonumber\\ &&- \hat{b}^+_3 \hat{b}^+_2 \hat{b}_3+ \hat{b}^+_4 \hat{b}^+_1 \hat{b}_1 - \hat{b}^+_4 \hat{b}^+_3 \hat{b}_3 +\hat{b}^+_4 \hat{b}^+_3 \hat{b}^+_1- \hat{b}^+_3 \hat{b}^+_2 \hat{b}^+_1 )) \nonumber\\ &&+\frac{1}{2}\left( ( \Gamma_{1,2} + \Gamma_1^+- \Gamma_{2,1} - \Gamma_1 - \Gamma_{1,1}) \left(\hat{b}_1^+ \hat{b}_1+\hat{b}_3^+ \hat{b}_3\right)\right.\nonumber\\ &&\left. + (\Gamma_{1,2} + \Gamma_2^+- \Gamma_{2,1}- \Gamma_2 ) \left(\hat{b}_2^+ \hat{b}_2+\hat{b}_4^+ \hat{b}_4\right)\right. \nonumber\\ &&\left.+ \Gamma_{1,1} \left( 2 \hat{b}_3^+ \hat{b}_1^+ \hat{b}_3 \hat{b}_1 -\hat{b}_1^+ \hat{b}_1^+ \hat{b}_1 \hat{b}_1 - \hat{b}_3^+ \hat{b}_3^+ \hat{b}_3 \hat{b}_3\right)\right)\nonumber\\ &&+ (\Gamma_{1,2} + \Gamma_{2,1})\left( \hat{b}_3^+ \hat{b}_1^+ \hat{b}_4 \hat{b}_2 +\hat{b}_4^+ \hat{b}_2^+ \hat{b}_3 \hat{b}_1\right) +\Gamma_{1,2}\left( \hat{b}_3^+ \hat{b}_2^+ \hat{b}_3 \hat{b}_2 + \hat{b}_4^+ \hat{b}_1^+ \hat{b}_4 \hat{b}_1\right) \nonumber\\ &&-\Gamma_{2,1}\left( \hat{b}_2^+ \hat{b}_1^+ \hat{b}_2 \hat{b}_1 + \hat{b}_4^+ \hat{b}_3^+ \hat{b}_4 \hat{b}_3\right) + (\Gamma_{1,2} +\Gamma_1^+) \hat{b}_3^+ \hat{b}_1^+ + (\Gamma_{1,2} + \Gamma_{2}^+)\hat{b}_4^+ \hat{b}_2^+ \nonumber\\ &&+ \Gamma_{1,2} \left(\hat{b}_3^+ \hat{b}_2^+ \hat{b}_1^+ \hat{b}_2 + \hat{b}_4^+ \hat{b}_2^+ \hat{b}_1^+ \hat{b}_1 + \hat{b}_4^+ \hat{b}_3^+ \hat{b}_1^+ \hat{b}_4 + \hat{b}_4^+ \hat{b}_3^+ \hat{b}_2^+ \hat{b}_3 + \hat{b}_4^+ \hat{b}_3^+ \hat{b}_2^+ \hat{b}_1^+\right). \end{eqnarray} \tag{ 19 }$$

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) with $\skew3\hat{\mathcal {B}}_{\mathrm { VR}}$ given by (19). In order to shorten the expressions we further restrict the rates Γ⁺₁ = Γ⁺₂ = Γ_1,2 = Γ_2,1 = Γ⁺ and Γ₁ = Γ₂ = Γ_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:

$$\begin{eqnarray} &&\begin{array}{ll} \displaystyle {\langle a_2\rangle}=-\frac{2\mathrm{ i} g \left(-1+\mathrm{e}^{-\frac{1}{2} t \left(-\Gamma ^++\Gamma +2 \mathrm{i} \Omega \right)}\right)}{-\Gamma ^++\Gamma +2\mathrm{ i} \Omega },\\ \displaystyle {\langle a_2a_2\rangle}=-\frac{4\,g^2 \left(-1+\mathrm{e}^{-\frac{1}{2} t \left(-\Gamma ^++\Gamma +2 \mathrm{i} \Omega \right)}\right)^2}{\left(-\Gamma ^++\Gamma +2 \mathrm{i} \Omega \right)^2}. \end{array} \end{eqnarray} \tag{ 20 }$$

Since the remaining expressions for the time evolution of the occupations are too lengthy we write only their steady-state values

$$\begin{eqnarray*} {\langle a_1^\dagger a_1\rangle}&=&\frac{2 \Gamma ^+ \left(\Gamma ^++\Gamma \right)+8 \Gamma \left(\Gamma -\Gamma ^+\right) g^2}{\left(\Gamma -3 \Gamma ^+\right) \left(\Gamma ^++\Gamma \right)}, \\ {\langle a_2^\dagger a_2\rangle}&=&\frac{2 \Gamma ^+ \left(\Gamma ^++\Gamma \right)+4 \left(\Gamma -\Gamma ^+\right) \left(3 \Gamma ^++\Gamma \right) g^2}{\left(\Gamma -3 \Gamma ^+\right) \left(\Gamma ^++\Gamma \right)}. \end{eqnarray*}$$

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 non-unitary part of the evolution.

As in the previous section we decompose the fermionic Liouvillian $\skew3\hat{\mathcal {F}}$ in terms of creation and annihilation maps (5)

$$\begin{eqnarray} &&\begin{array}{ll} \displaystyle \skew3\hat{\mathcal{F}}=\sum_{m,n}\skew3\hat{\mathcal{F}\,}^{(m,n)}\skew3\hat{\mathcal{P}}^{m+n},\\ \displaystyle \skew3\hat{\mathcal{F}\,}^{(m,n)}=\sum_{\underline{j},\underline{k}}L^{(m,n)}_{j_1\ldots j_n,k_1\ldots k_m}\skew3\hat{f}_{j_1}^+\ldots\skew3\hat{f}_{j_m}^+\skew3\hat{f}_{k_1}\ldots\skew3\hat{f}_{k_n}. \end{array} \end{eqnarray} \tag{ 21 }$$

The fermionic creation maps $\skew3\hat{f}_{j}^+$ defined in (5) left-annihilate the identity operator, namely

$$\begin{equation} \left ( \mathbbm{1}\right | \skew3\hat{f}^+_j=0, \end{equation} \tag{ 22 }$$

which again implies, by using the trace preservation property $\left ( \mathbbm{1}\right | \skew3\hat{\mathcal {F}}=0$, that in the sum (21) we always have m > 0. We define the anti-symmetric p-point correlation tensors as

$$\begin{equation} Z^{(p,\alpha)}_{j_1,j_2,\ldots,j_p}(t)=\left ( 1\right | \skew3\hat{f}_{j_1}\skew3\hat{f}_{j_2}\ldots\skew3\hat{f}_{j_p}\skew3\hat{\mathcal{P}\,}^\alpha \left | \rho(t)\right \rangle \end{equation} \tag{ 23 }$$

with α = 0,1 determining the presence of the super-operator $\skew3\hat{\mathcal {P}}$. An additional correlation tensor is necessary because of a parity super-operator in the fermionic Liouvillian (21). The time evolution of correlation tensors Z^(p,α) is induced by

$$\begin{equation} \frac{{\rm d}}{{\rm d} t} Z^{(p,\alpha)}_{j_1,j_2,\ldots,j_p}(t)=\left ( 1\right | \skew3\hat{f}_{j_1}\skew3\hat{f}_{j_2}\ldots\skew3\hat{f}_{j_p}\skew3\hat{\mathcal{P}\,}^\alpha\skew3\hat{\mathcal{F}}\left | \rho(t)\right \rangle. \end{equation} \tag{ 24 }$$

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-operators $\skew3\hat{\mathcal {F}\,}^{(m,n)}$ and $\skew3\hat{\mathcal {F}\,}^{(m,n)}\skew3\hat{\mathcal {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-operator $\skew3\hat{\mathcal {P}}$ in Z^p,1, which satisfies the following relations:

$$\begin{equation} \skew3\hat{\mathcal{P}}\skew3\hat{f}_{j}^+=-\skew3\hat{f}_{N+j}\skew3\hat{\mathcal{P}},\quad \skew3\hat{\mathcal{P}}\skew3\hat{f}_{N+j}^+=-\skew3\hat{f}_{j}\skew3\hat{\mathcal{P}},\quad \skew3\hat{\mathcal{P}\,}^2=\mathbbm{1}. \end{equation} \tag{ 26 }$$

Hence, the super-operators $\skew3\hat{\mathcal {F}\,}^{(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-vanishing $\skew3\hat{\mathcal {F}\,}^{(m,n)}$ in (21) ; α = 0;
• m + n : even and m − n ⩽ 0 for all non-vanishing $\skew3\hat{\mathcal {F}\,}^{(m,n)}$ in (21) ; α = 1;
• ( m − n = K, where K is an odd integer for all non-vanishing $\skew3\hat{\mathcal {F}\,}^{(m,n)}\skew3\hat{\mathcal {P}}$ in (21) with odd m + n) or ( m − n = 0 for all non-vanishing $\skew3\hat{\mathcal {F}\,}^{(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–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

$$\begin{equation} L_{\underline{j},\underline{k}}=c^\dagger_{j_1}c^\dagger_{j_2}\ldots c^\dagger_{j_\mu} c_{k_1}c_{k_2}\ldots c_{k_\nu}, \end{equation} \tag{ 27 }$$

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

$$\begin{equation} G^{(\mu,\nu;\mu',\nu')}_{\underline{j},\underline{k};\underline{j'},\underline{k'}} +\Phi(\mu,\nu,\mu'\nu')\bar{G}^{(\nu,\mu;\nu',\mu')}_{\underline{k},\underline{j};\underline{k'},\underline{j'}} =0, \end{equation} \tag{ 28 }$$

where we define the phase factor $\Phi (\mu ,\nu ,\mu '\nu ')=(-1)^{\lfloor \frac {\mu +\nu }{2}\rfloor +\lfloor \frac {\mu '+\nu '}{2}\rfloor +1 + \nu \nu ' + \mu ' \nu + \mu ' \nu ' + \mu \mu ' + \mu \nu + \mu \nu '}$. In the simplest nontrivial example with the coupling operator basis L_j,k∈{c^†_jc_k,c^†_jc^†_k,c_jc_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=\sum _{j,k}A_{jk}c_j^\dagger 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.

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

$$\begin{equation} H_{\rm Hub}=\sum_{j,k=1}^N t_{jk}c_{j,\sigma}^\dagger c_{k\sigma}+ \sum_{j=1}^NU_j c^\dagger_{j,\uparrow}c^\dagger_{j,\downarrow}c_{j,\uparrow}c_{j,\downarrow}. \end{equation} \tag{ 29 }$$

The Fermi–Hubbard model introduced in [30–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

$$\begin{equation} H_{\rm int}=U c^\dagger_{1}c^\dagger_{2}c_{1}c_{2}. \end{equation} \tag{ 30 }$$

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-operators

$$\begin{equation} \skew3\hat{\rm ad}_{H_{\rm int}}=U(\skew3\hat{f}_1^+\skew3\hat{f}_4\skew3\hat{f}_1\skew3\hat{f}_2-\skew3\hat{f}_2^+\skew3\hat{f}_3\skew3\hat{f}_1\skew3\hat{f}_2\skew3\hat{f}^+_3\skew3\hat{f}_3\skew3\hat{f}_4\skew3\hat{f}_2-\skew3\hat{f}^+_4\skew3\hat{f}_3\skew3\hat{f}_4\skew3\hat{f}_1) +\mbox{rest}, \end{equation} \tag{ 31 }$$

where the rest contains the hierarchy-preserving terms. Next, we wish to add dissipation of the form (13) with the coupling operator basis $L_{\underline {j},\underline {k}}\in \{c_{1},c_2, c^\dagger _{2}c_2 c_1,c^\dagger _{1}c_1 c_2, c^\dagger _{1}c_2^\dagger c_2,c^\dagger _{2}c_1^\dagger c_1,c_1^\dagger ,c_2^\dagger \}$ 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 two-point correlations vanish aside from the occupation numbers ${\langle c_{j,\sigma }^\dagger c_{j,\sigma }\rangle }=\frac {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.

Zoom In Zoom Out Reset image size

We anew employ the same strategy as for bosons and fermions and rewrite the complete mixed Liouvillian with $\skew3\hat{\mathcal {B}}^{(m',n')}$ and $\skew3\hat{\mathcal {F}\,}^{(m,n)}$ as defined in (6) and (21), respectively,

$$\begin{equation} \skew3\hat{\mathcal{L}}=\sum_{m,n,m'n'}\skew3\hat{\mathcal{B}}^{(m',n')}\skew3\hat{\mathcal{F}\,}^{(m,n)}\skew3\hat{\mathcal{P}}^{m+n}. \end{equation} \tag{ 32 }$$

The fermionic and the bosonic Liouvillian commute, $[\skew3\hat{\mathcal {B}},\skew3\hat{\mathcal {F}}]=0$. We define the mixed (p + q)-point correlation tensor as

$$\begin{equation} Z^{(p,q,\alpha)}_{j_1,\ldots,j_p,k_1,\ldots,k_q}(t)=\left ( 1\right | \hat{b}_{j_1}\ldots\hat{b}_{j_p}\skew3\hat{f}_{k_1}\ldots\skew3\hat{f}_{k_q}\skew3\hat{\mathcal{P}\,}^\alpha\left | \rho(t)\right \rangle \end{equation} \tag{ 33 }$$

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

$$\begin{equation} \frac{{\rm d}}{{\rm d} t} Z^{(p,q,\alpha)}_{j_1,\ldots j_p,k_1,\ldots k_q}(t)=\left ( 1\right | \hat{b}_{j_1}\ldots\hat{b}_{j_p}\skew3\hat{f}_{k_1}\ldots\skew3\hat{f}_{k_q}\skew3\hat{\mathcal{P}\,}^\alpha\skew3\hat{\mathcal{L}}\left | \rho(t)\right \rangle. \end{equation} \tag{ 34 }$$

The action of the mixed Liouvillian (32) on the correlation tensors can be explained by the following diagrams:

Again due to an additional phase super-operator in the diagram (36), we have the opposite sign of m and n in comparison to the diagram (35). On the other hand, the bosonic part (i.e. m' and n') does not change the sign. Taking into account all transitions depicted in the diagrams (35) and (36) we deduce that the equations for the correlators Z^(p,q,α) are closed under the following conditions (mixed closeness conditions):

• m + n : even and n − m ⩽ m' − n' for all non-vanishing $\skew3\hat{\mathcal {F}\,}^{(m,n)}\skew3\hat{\mathcal {B}}^{(m',n')}$ in (32) ; α = 0;
• m + n : even and m − n ⩽ m' − n' for all non-vanishing $\skew3\hat{\mathcal {F}\,}^{(m,n)}\skew3\hat{\mathcal {B}}^{(m',n')}$ in (32) ; α = 1;
• ( |m₁ − m₂ − n₁ + n₂| ⩽ m₁' + m₂' − n₁' − n₂' for all non-vanishing paris $\skew3\hat{\mathcal {F}\,}^{(m_1,n_1)}\skew3\hat{\mathcal {B}}^{(m_1',n_1')}\skew3\hat{\mathcal {P}}$ and $\skew3\hat{\mathcal {F}\,}^{(m_2,n_2)}\skew3\hat{\mathcal {B}}^{(m_2',n_2')}\skew3\hat{\mathcal {P}}$ in (32) with odd m₁ + n₁ and m₂ + n₂) or ( |m − n| ⩽ m' − n' for all non-vanishing $\skew3\hat{\mathcal {F}\,}^{(m,n)}\skew3\hat{\mathcal {B}}^{(m',n')}$ in (32) with even m + n) ; α = 0,1.

If the only non-vanishing $\skew3\hat{\mathcal {F}\,}^{(m,n)}$ in (32) has m + n = 0, above conditions reduce to simple bosonic closeness conditions and similarly if the only non-vanishing $\skew3\hat{\mathcal {B}}^{(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

$$\begin{equation} {\bf G}=\left(\begin{array}{@{}cccc@{}} G_{11} & 0 & G_{13} & 0\\ 0 & G_{11} & 0 & \bar{G}_{13}\\ \bar{G}_{13} & 0 & G_{44}+2\delta & 0\\ 0 & G_{13} & 0 & G_{44} \end{array}\right). \end{equation} \tag{ 37 }$$

The positivity of G implies |G₁₃| < G₁₁G₄₄ and G₁₁,G₄₄ > 0, whereas the stability of the system implies δ > 0.

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

$$\begin{equation} H_{\rm Rabi}=\omega a^\dagger a+\Delta c^\dagger c + g (c+c^\dagger)(a+a^\dagger). \end{equation} \tag{ 38 }$$

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-operators

$$\begin{eqnarray}&&\fl \skew3\hat{\rm ad}_{{H}_{\rm Rabi}}=\omega(\hat{b}_1^+\hat{b}_1-\hat{b}_2^+\hat{b}_2)+\Delta(\skew3\hat{f}_1^+\skew3\hat{f}_1-\skew3\hat{f}_2^+\skew3\hat{f}_2)+\frac{g}{\sqrt{2}}\Big(\hat{b}^+_1(\skew3\hat{f}_1-\skew3\hat{f}_2+\skew3\hat{f}_2^++\skew3\hat{f}_1^+) \nonumber\\ &&-\hat{b}^+_2(\skew3\hat{f}_1-\skew3\hat{f}_2+\skew3\hat{f}_1^+-\skew3\hat{f}_2^+)+2(\hat{b}_1+\hat{b}_2)(\skew3\hat{f}^+_2-\skew3\hat{f}_1^+)\Big)\skew3\hat{\mathcal{P}}. \end{eqnarray} \tag{ 39 }$$

We observe that the closeness conditions are not satisfied due to the terms of the form $\skew3\hat{f}^+_j\hat{b}_k\skew3\hat{\mathcal {P}}$, 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

$$\begin{equation} {\bf G}=\left(\begin{array}{@{}cccc@{}} G_{11} & G_{1,2} & 0 & -2\,{\rm i} \, g \\ G_{21} & G_{22} & 0 & -2\,{\rm i} \, g \\ 0 & 0 & G_{33} & G_{3,4} \\ 2\,{\rm i} \,g & 2\,{\rm i}\,g & G_{43} & G_{44} \end{array}\right). \end{equation} \tag{ 40 }$$

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}=\Gamma >2 \sqrt {2}|g|$, G₃₃ = Γ + 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 ( $\Gamma , \Gamma +2\delta , \Gamma - 2\sqrt {2}\,g, \Gamma + 2\sqrt {2}\,g$). The full Liouvillian $\skew3\hat{\mathcal {L}}_{\mathrm { Rabi}}=-{\mathrm { i}}\,\skew3\hat{\mathrm { ad}}_{H_{\mathrm { Rabi}}}+\skew3\hat{\mathcal {D}}_{\mathrm { Rabi}}$ expressed in terms of bosonic (4) and fermionic (5) creation and annihilation maps is

$$\begin{eqnarray} &&\fl \skew3\hat{\mathcal{L}}_{\rm Rabi}=-{\rm i}\Bigg( \omega(\hat{b}_1^+\hat{b}_1-\hat{b}_2^+\hat{b}_2)+\Delta(\skew3\hat{f}_1^+\skew3\hat{f}_1-\skew3\hat{f}_2^+\skew3\hat{f}_2)+\sqrt{2} \,g\, (\hat{b}_1^+ \skew3\hat{f}_1 - \hat{b}_1^+ \skew3\hat{f}_2 - \hat{b}_2^+ \skew3\hat{f}_1 + \hat{b}_2^+ \skew3\hat{f}_2 \nonumber\\ &&+ \hat{b}_1^+\skew3\hat{f}_1^+ - \hat{b}_1^+ \skew3\hat{f}_2^++ \hat{b}_2^+ \skew3\hat{f}^+_1 - \hat{b}_2^+ \skew3\hat{f}_2^+)\skew3\hat{\mathcal{P}}\Bigg) +\Gamma(\hat{b}^+_1\hat{b}^+_2-\skew3\hat{f}_1^+\skew3\hat{f}_1-\skew3\hat{f}_2^+\skew3\hat{f}_2)-\delta(\hat{b}_1^+\hat{b}_1+\hat{b}_2^+\hat{b}_2). \nonumber\\ \end{eqnarray} \tag{ 41 }$$

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:

$$\begin{eqnarray} &&\begin{array}{lll} \displaystyle {\langle c\rangle}={\langle a\rangle}={\langle c^\dagger c-cc^\dagger\rangle}=0, \\ \displaystyle {\langle ac\rangle}=\frac{2\, {\rm i} \,g \left(-1+\mathrm{e}^{-t (\Gamma +\delta +{\rm i} (\Delta +\omega ))}\right)}{\Gamma +\delta +{\rm i} (\Delta +\omega )},\\ \displaystyle {\langle a^\dagger c\rangle} = \frac{- 2\,{\rm i}\, g \left(-1+\mathrm{e}^{-t (\Gamma +\delta +{\rm i} (\Delta -\omega ))}\right)}{\Gamma +\delta +{\rm i} (\Delta-{\rm i}\omega) }. \end{array} \end{eqnarray} \tag{ 42 }$$

Since the remaining expressions for the time evolution of correlations are too lengthy we provide only their steady-state values.

$$\begin{eqnarray} &&\begin{array}{ll} \displaystyle {\langle a^\dagger a\rangle}=\frac{\Gamma +\frac{4\,g^2 (\Gamma +\delta )}{(\Gamma +\delta )^2+(\Delta -\omega )^2}+\frac{4\,g^2 (\Gamma +\delta )}{(\Gamma +\delta )^2+(\Delta +\omega )^2}}{2 \delta },\\ \displaystyle {\langle a^2\rangle}=-\frac{4\,g^2 (\Gamma +\delta +{\rm i} \omega )}{(\delta +{\rm i} \omega ) (\Gamma +\delta -{\rm i} (\Delta -\omega )) (\Gamma +\delta +{\rm i} (\Delta +\omega ))}. \end{array} \end{eqnarray} \tag{ 43 }$$

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:

$$\begin{eqnarray} L_1&=&\sqrt{\Gamma_{\rm flip}}(c-c^\dagger),\quad L_2=\sqrt{\Gamma_{\rm amp}} a,\nonumber\\ L_3&=&\sqrt{\Gamma_{\rm c_1}}\left(\frac{\mathrm{i}}{\sqrt{2}}(c+c^\dagger)+a^\dagger\right),\\ L_4&=&\sqrt{\Gamma_{\rm c_2}}\left(\frac{\mathrm{i}}{\sqrt{2}}(c+c^\dagger)-a^\dagger\right).\nonumber \end{eqnarray} \tag{ 44 }$$

The first Lindblad operator describes a bit-phase flip channel with the rate $\Gamma _{\mathrm { flip}}=\frac {\Gamma }{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 $\Gamma _{\mathrm { c_1}}=\frac {1}{2}\Gamma +\sqrt {2} g$ and $\Gamma _{\mathrm { c_2}}=\frac {1}{2}\Gamma -\sqrt {2} g$. 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². 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 $\Gamma _{\mathrm { flip}}>\frac {\Gamma }{2}$ and Γ_amp > Γ + 2δ, respectively. Hence, only the rates corresponding to collective channels have to be fine-tuned, i.e. $\Gamma _{\mathrm { c_1}}-\Gamma _{\mathrm { c_2}}=2\sqrt {2}g$ and Γ_c1, Γ_c2 > 0.

A general bosonic Hamiltonian can be written in the form

$$\begin{equation} H=\sum_{\mu,\nu}\sum_{\underline{j}\underline{k}}H_{\underline{j},\underline{k}}^{(\mu,\nu)}a^\dagger_{j_1}a^\dagger_{j_2}\ldots a^\dagger_{j_m}a_{k_1}a_{k_2}\ldots d_{k_n}, \end{equation} \tag{ A.1 }$$

where the indices μ and ν determine the lengths of the vectors $\underline {j}$ and $\underline {k}$, respectively. It is clear that due to the trace preservation and the property $\left \langle \mathbbm{1}\right | \hat{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 $\skew3\hat{\mathrm { ad}}_{H}\,\rho =[H,\rho ]$ that contain only one bosonic creation map $\hat{b}_{j_p}^+$ or $\hat{b}_{N+j_p}^+$

$$\begin{eqnarray} &&\fl \skew3\hat{\rm ad}_{H}=\sum_{\mu,\nu}\sum_{p=1}^{\mu}\sum_{\underline{j},\underline{k}}H^{(\mu,\nu)}_{\underline{j},{\underline{k}}}\hat{b}_{j_p}^+\prod_{i=1,i\neq p}^\mu\hat{b}_{N+j_i}\prod_{i=1}^{\nu}\hat{b}_{k_i}-\sum_{\mu,\nu}\sum_{p=1}^{\nu}\sum_{\underline{j},\underline{k}}H^{(\mu,\nu)}_{\underline{j},{\underline{k}}}\hat{b}_{N+k_p}^+\nonumber\\ &&\times\prod_{i=1,i\neq p}^{\nu}\hat{b}_{k_i}\prod_{i=1}^{\mu}\hat{b}_{N+j_i}+{\rm rest}. \end{eqnarray} \tag{ A.2 }$$

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\underline {j},P'\underline {k}}=H_{\underline {j},\underline {k}}$, where P and P' represent arbitrary permutations of elements of vectors $\underline {j}$ and $\underline {k}$, respectively, the explicit sums in equation (A.2) cannot vanish. Therefore, the map $\skew3\hat{\mathrm { ad}}_{H}$ always contains terms with only one creation super-operator and ν + μ − 1 annihilation super-operators, which leads to the following closeness condition: $M_{\mathrm { max}}=\max (\nu +\mu -2)\leqslant 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.

In this subsection we discuss the closeness conditions for general fermionic Hamiltonians

$$\begin{equation} H=\sum_{\mu,\nu}\sum_{\underline{j}\underline{k}}H_{\underline{j},\underline{k}}^{(\mu,\nu)}c^\dagger_{j_1}c^\dagger_{j_2}\ldots c^\dagger_{j_\mu}c_{k_1}c_{k_2}\ldots c_{k_\nu}. \end{equation} \tag{ A.3 }$$

Explicitly writing only terms with exactly one fermionic creation map $\skew3\hat{f}_{j_p}^+$ or $\skew3\hat{f}_{N+j_p}^+$ the adjoint map of the fermionic Hamiltonian (A.3) reads

$$\begin{eqnarray} &&\fl \skew3\hat{\rm ad}_H=\sum_{\mu,\nu}\sum_{p=1}^{\mu}\frac{(-1)^{\lfloor (\nu+\mu)/2\rfloor+\mu+p+1}}{2^{(\nu+\mu)/2}}\sum_{\underline{j},\underline{k}}H^{(\mu,\nu)}_{\underline{j},{\underline{k}}}\skew3\hat{f}_{j_p}^+\prod_{i=1;i\neq p}^\mu\skew3\hat{f}_{N+j_i}\prod_{i=1}^{\nu}\skew3\hat{f}_{k_i}\skew3\hat{\mathcal{P}}^{\mu+\nu} \nonumber\\ &&-\sum_{\mu,\nu}\sum_{p=1}^{\mu}\frac{(-1)^{\mu-p }}{2^{(\nu+\mu)/2}}\sum_{\underline{j},\underline{k}}H^{(\mu,\nu)}_{\underline{j},{\underline{k}}}\prod_{i=\nu}^{1}\skew3\hat{f}_{k_i}\prod_{i=\mu;i\neq p}^1\skew3\hat{f}_{N+j_i} \, \skew3\hat{f}_{j_p}^+\skew3\hat{\mathcal{P}}^{\mu+\nu} \nonumber\\ &&+\sum_{\mu,\nu}\sum_{p=1}^{\nu}\frac{(-1)^{\lfloor (\mu+\nu)/2\rfloor+\mu+\nu-p }}{2^{(\nu+\mu)/2}}\sum_{\underline{j},\underline{k}}H^{(\mu,\nu)}_{\underline{j},{\underline{k}}}\prod_{i=1}^{\mu}\skew3\hat{f}_{N+j_i}\prod_{i=1;i\neq p}^\nu\skew3\hat{f}_{k_i} \, \skew3\hat{f}_{N+k_p}^+\skew3\hat{\mathcal{P}}^{\mu+\nu} \nonumber\\ &&-\sum_{\mu,\nu}\sum_{p=1}^{\nu}\frac{(-1)^{\nu+\mu-p+1}}{2^{(\nu+\mu)/2}}\sum_{\underline{j},\underline{k}}H^{(\mu,\nu)}_{\underline{j},{\underline{k}}}\skew3\hat{f}_{N+k_p}^+\prod_{i=\nu;i\neq p}^1\skew3\hat{f}_{k_i}\prod_{i=\mu}^{1}\skew3\hat{f}_{N+j_i}\skew3\hat{\mathcal{P}}^{\mu+\nu}+{\rm rest}. \end{eqnarray} \tag{ A.4 }$$

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

$$\begin{eqnarray}&&\fl \skew3\hat{\rm ad}_{H}=\sum_{\mu,\nu}\sum_{p=1}^{\mu}\frac{(-1)^{\lfloor (\nu+\mu)/2\rfloor+\mu+p+1}}{2^{(\nu+\mu)/2-1}}\sum_{\underline{j},\underline{k}}H^{(\mu,\nu)}_{\underline{j},{\underline{k}}}\skew3\hat{f}_{j_p}^+\prod_{i=1;i\neq p}^\mu\skew3\hat{f}_{N+j_i}\prod_{i=1}^{\nu}\skew3\hat{f}_{k_i}\skew3\hat{\mathcal{P}}^{\mu+\nu} \nonumber\\ &&+\sum_{\mu,\nu}\sum_{p=1}^{\nu}\frac{(-1)^{\lfloor (\mu+\nu)/2\rfloor+\mu+\nu-p}}{2^{(\nu+\mu)/2-1}}\sum_{\underline{j},\underline{k}}H^{(\mu,\nu)}_{\underline{j},{\underline{k}}}\prod_{i=1}^{\mu}\skew3\hat{f}_{N+j_i}\prod_{i=1;i\neq p}^\nu\skew3\hat{f}_{k_i} \, \skew3\hat{f}_{N+k_p}^+\skew3\hat{\mathcal{P}}^{\mu+\nu} +{\rm rest}.\nonumber\\ \end{eqnarray} \tag{ A.5 }$$

Using the symmetry $H_{P\underline {j},P'\underline {k}}=\varepsilon _{P\underline {j}}\varepsilon _{P'\underline {k}}H_{\underline {j},\underline {k}}$, with P and P' being arbitrary permutations of elements of vectors $\underline {j}$ and $\underline {k}$, respectively, and $\varepsilon _{\underline {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 maps

$$\begin{eqnarray} &&\fl \skew3\hat{\rm ad}_{H}=\sum_{\mu,\nu}\sum_{p=1}^{\mu}\frac{(-1)^{\lfloor (\nu+\mu)/2\rfloor+\mu+p+1}}{2^{(\nu+\mu)/2-1}}\sum_{\underline{j},\underline{k}}H^{(\mu,\nu)}_{\underline{j},{\underline{k}}}(1-(-1)^{\mu+\nu})\prod_{i=1}^\mu\skew3\hat{f}_{N+j_i}\prod_{i=1}^{\nu}\skew3\hat{f}_{k_i}\skew3\hat{\mathcal{P}}^{\mu+\nu},+{\rm rest},\nonumber\\ \end{eqnarray} \tag{ 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

$$\begin{eqnarray*} &&\fl \skew3\hat{\rm ad}_{H}=\sum_{\mu,\nu}\sum_{p=1}^{\mu}\frac{(-1)^{\lfloor (\nu+\mu)/2\rfloor+\mu+p+1}}{2^{(\nu+\mu)/2}}(1+(-1)^{\mu+\nu})\sum_{\underline{j},\underline{k}}H^{(\mu,\nu)}_{\underline{j},{\underline{k}}}\skew3\hat{f}_{N+j_p}\prod_{i=1;i\neq p}^\mu\skew3\hat{f}_{j_i}^+\prod_{i=1}^{\nu}\skew3\hat{f}_{N+k_i}^+\skew3\hat{\mathcal{P}}^{\mu+\nu} \nonumber\\ &&+\sum_{\mu,\nu}\sum_{p=1}^{\nu}\frac{(-1)^{\lfloor (\mu+\nu)/2\rfloor+\mu+\nu-p}}{2^{(\nu+\mu)/2}}(1\!+\!(-1)^{\mu+\nu})\sum_{\underline{j},\underline{k}}H^{(\mu,\nu)}_{\underline{j},{\underline{k}}}\prod_{i=1}^{\mu}\skew3\hat{f}_{j_i}^+\!\prod_{i=1;i\neq p}^\nu\skew3\hat{f}_{N+k_i} \skew3\hat{f}_{k_p}\skew3\hat{\mathcal{P}}^{\mu+\nu} \nonumber\\ &&+{\rm rest}, \end{eqnarray*}$$

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.

In this subsection we determine the closeness condition for the mixed Hamiltonian of the form

$$\begin{equation} H=\sum_{\mu,\nu,\mu',\nu'}\sum_{\underline{j},\underline{k},\underline{t},\underline{v}}H^{(\mu,\nu,\mu',\nu')}_{\underline{j},\underline{k},\underline{t},\underline{v}}c^\dagger_{j_1}\ldots c^\dagger_{j_\mu}c_{k_1}\ldots c_{k_\nu}a^\dagger_{t_1}\ldots a^\dagger_{t_\mu}a_{v_1}\ldots a_{v_\nu} \end{equation} \tag{ A.7 }$$

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.

In general the bosonic coupling operators can be written as

$$\begin{equation} L_\alpha=a^\dagger_{j_1}a^\dagger_{j_2}\ldots a^\dagger_{j_m} a_{k_1}a_{k_2}\ldots a_{k_n}. \end{equation} \tag{ B.2 }$$

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 numbers 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 maps

$$\begin{eqnarray} &&\fl \skew3\hat{\mathcal{D}}=\frac{1}{2}\sum_{\underline{j},\underline{k},\underline{t},\underline{v}}G^{(\mu,\nu;\mu',\nu')}_{\underline{j},\underline{k};\underline{t},\underline{v}}\left(\sum_{p=1}^\mu \hat{b}_{j_p}^\dagger \prod_{i=1}^\nu\hat{b}_{k_i}\prod_{i=1}^{\mu'}\hat{b}_{t_i}\prod_{i=1,i\neq p}^\mu\hat{b}_{N+j_i}\prod_{i=1}^{\nu'}\hat{b}_{N+v_i}\right. \nonumber\\ &&- \sum_{p=1}^{\nu'} \hat{b}_{v_p}^\dagger \prod_{i=1}^\nu\hat{b}_{k_i}\prod_{i=1}^{\mu'}\hat{b}_{t_i}\prod_{i=1}^\mu\hat{b}_{N+j_i}\prod_{i=1, i\neq p}^{\nu'}\hat{b}_{N+v_i} \nonumber\\ &&+ \sum_{p=1}^{\mu'} \hat{b}_{N+t_p}^\dagger \prod_{i=1}^\nu\hat{b}_{k_i}\prod_{i=1,i\neq p}^{\mu'}\hat{b}_{t_i}\prod_{i=1}^\mu\hat{b}_{N+j_i}\prod_{i=1}^{\nu'}\hat{b}_{N+v_i} \nonumber\\ &&\left.- \sum_{p=1}^\nu \hat{b}_{N+k_p}^\dagger \prod_{i=1 ,i\neq p}^\nu\hat{b}_{k_i}\prod_{i=1}^{\mu'}\hat{b}_{t_i}\prod_{i=1}^\mu\hat{b}_{N+j_i}\prod_{i=1}^{\nu'}\hat{b}_{N+v_i} \right)+{\rm rest}. \end{eqnarray} \tag{ B.3 }$$

Since the tensor G is symmetric with respect to permutations of elements of each multi-index $\underline {j},\underline {k},\underline {t}$ and $\underline {v}$, we obtain the conditions

$$\begin{eqnarray} &&\begin{array}{ll} \fl\displaystyle \sum_{\underline{j},\underline{k},\underline{t},\underline{v}}\left(G^{(\mu,\nu;\mu',\nu')}_{\underline{j},\underline{k};\underline{t},\underline{v}}-G^{(\nu',\nu;\mu',\mu)}_{\underline{v},\underline{k};\underline{t},\underline{j}}\right)\sum_{p=1}^{\mu} \hat{b}_{j_p}^\dagger \prod_{i=1}^{\nu}\hat{b}_{k_i}\prod_{i=1}^{\mu'}\hat{b}_{t_i}\prod_{i=1 , i\neq p}^{\mu}\hat{b}_{N+j_i}\prod_{i=1}^{\nu'}\hat{b}_{N+v_i} =0, \\ \fl\displaystyle \sum_{\underline{j},\underline{k},\underline{t},\underline{v}}\left(G^{(\mu,\nu;\mu',\nu')}_{\underline{j},\underline{k};\underline{t},\underline{v}}-G^{(\mu,\mu';\nu,\nu')}_{\underline{j},\underline{t};\underline{k},\underline{v}}\right)\sum_{p=1}^{\nu} \hat{b}_{N+k_p}^\dagger \prod_{i=1 , i\neq p}^{\nu}\hat{b}_{k_i}\prod_{i=1}^{\mu'}\hat{b}_{t_i}\prod_{i=1}^{\mu}\hat{b}_{N+j_i}\prod_{i=1}^{\nu'}\hat{b}_{N+v_i} =0. \end{array} \end{eqnarray} \tag{ B.4 }$$

Now we introduce the vectors $\underline {x}=(\underline {j},\underline {v})$ and $\underline {y}=(t_2,t_3\ldots ,t_{m'},k_1,k_2,\ldots ,k_{n})$ and extract from equation (B.4) the following symmetry conditions for the rate matrix  G:

$$\begin{eqnarray} 0&=&\sum_{\mu+\nu'=n_x}\sum_{\nu+\mu'=n_y}\sum_{P \underline{x},P'\underline{y}}\mu\left( G^{(\mu,\nu;\mu',\nu')}_{\underline{j},\underline{k};\underline{t},\underline{v}}-G^{(\mu,\mu';\nu,\nu')}_{\underline{j},\underline{t};\underline{k},\underline{v}} \right)\nonumber\\ &=&\sum_{\mu+\nu'=n_x}\sum_{\nu+\mu'=n_y}\sum_{P \underline{x},P'\underline{y}}\mu\left( G^{(\mu,\nu;\mu',\nu')}_{\underline{j},\underline{k};\underline{t},\underline{v}}-G^{(\nu',\mu';\nu,\mu)}_{\underline{v},\underline{t};\underline{k},\underline{j}} \right)\\ &=&\sum_{\mu+\nu'=n_x}\sum_{\nu+\mu'=n_y}\sum_{P \underline{x},P'\underline{y}}\mu\left( G^{(\mu,\nu;\mu',\nu')}_{\underline{j},\underline{k};\underline{t},\underline{v}}-\bar{G}^{(\nu,\mu;\nu',\mu')}_{\underline{k},\underline{j};\underline{v},\underline{t}} \right),\nonumber \end{eqnarray} \tag{ B.5 }$$

which have to be satisfied for all variations of $\underline {x}$, y and t₁. The simplest solution of equations (B.5) is if all terms in the sum vanish separately, namely if $G^{(\mu ,\nu ;\mu ',\nu ')}_{\underline {j},\underline {k};\underline {t},\underline {v}}-\bar {G}^{(\nu ,\mu ;\nu ',\mu ')}_{\underline {k},\underline {j};\underline {v},\underline {t}}$. The question remains whether this is in general the only solution which satisfies also the positivity constraint. Equations (B.5) are a necessary condition for the closeness of the hierarchy. However, they are sufficient only in 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 whether they can always be satisfied. Nevertheless, for a general quadratic noise determined by coupling operators from the set {a_ja_k, a^†_ja^†_k, a^†_ja_k}, equations (B.5) are sufficient for the closeness of the hierarchy and simplify to

$$\begin{eqnarray} &&\fl 0= 4(G^{(2,0;2,0)}_{j_1,j_2;t_1,t_2}-\bar{G}^{(0,2;0,2)}_{j_1,j_2;t_1,t_2})+ G^{(1,1;1,1)}_{j_1,j_2;t_1,t_2}+G^{(1,1;1,1)}_{t_2,j_2;t_1,j_1}-\bar{G}^{(1,1;1,1)}_{j_2,j_1;t_2,t_1}-\bar{G}^{(1,1;1,1)}_{j_2,t_2;j_1,t_1},\nonumber\\ &&\fl 0= G^{(1,1;2,0)}_{j_1,k_1;t_1,t_2}+G^{(1,1;2,0)}_{j_1,t_2;t_1,k_1}-\bar{G}^{(1,1;0,2)}_{k_1,j_1;t_1,t_2}+G^{(1,1;2,0)}_{t_2,j_1;k_1,t_1}+ G^{(0,2;1,1)}_{j_1,k_1;t_1,t_2}-\bar{G}^{(2,0;1,1)}_{j_1,k_1;t_2,t_1},\\ &&\fl 0= G^{(0,2;2,0)}_{k_1,k_2;t_1,t_2}-\bar{G}^{(2,0;0,2)}_{k_1,k_2;t_1,t_2}.\nonumber \end{eqnarray} \tag{ B.6 }$$

These equations have to be satisfied for all values of j₁,j₂,t₁,t₂,k₁,k₂ and determine the most general set of quadratic Markovian bosonic noise which satisfies the bosonic closeness condition.

In this section we study general fermionic dissipators of the form (B.1) with the coupling operators

$$\begin{equation} L_\alpha=c^\dagger_{j_1}c^\dagger_{j_2}\ldots c^\dagger_{j_m} c_{k_1}c_{k_2}\ldots c_{k_n}. \end{equation} \tag{ B.7 }$$

We shall study only the simplest equations leading to necessary and sufficient symmetries of G for the first fermionic closeness condition in the case of quadratic noise. As noted in the main text, we need to ensure that all terms with nontrivial phase super-operators vanish. After a tedious calculation and bookkeeping of the sign we obtain the following expression for the fermionic dissipator:

$$\begin{eqnarray} &&\fl \skew3\hat{\mathcal{D}}=\sum_{\underline{j},\underline{k},\underline{t},\underline{v}}\frac{(-1)^{\lfloor \frac{\mu+\nu}{2} \rfloor+\mu}}{2^{\frac{\mu+\nu+\mu'+\nu'}{2}}}G^{(\mu,\nu;\mu',\nu')}_{\underline{j},\underline{k};\underline{t},\underline{v}} \left(\sum_{p=1}^\mu(-1)^{p+1+\nu'+(m-1)(\mu'+\nu)}\skew3\hat{f}_{j_p}^+ \prod_{i=1}^\nu\skew3\hat{f}_{k_i}\prod_{i=1}^{\mu'}\skew3\hat{f}_{t_i}\prod_{i=1;i\neq p}^\mu\skew3\hat{f}_{\nu+j_i} \right. \nonumber\\ &&\times\prod_{i=1}^{\nu'}\skew3\hat{f}_{\nu+v_i}+\sum_{p=1}^{\nu'}(-1)^{p+\mu+\nu+\mu'+\nu'+\mu(\mu'+\nu)}\skew3\hat{f}_{v_p}^+ \prod_{i=1}^\nu\skew3\hat{f}_{k_i}\prod_{i=1}^{\mu'}\skew3\hat{f}_{t_i}\prod_{i=1}^\mu\skew3\hat{f}_{\nu+j_i}\prod_{i=1; i\neq p}^{\nu'}\skew3\hat{f}_{\nu+v_i} \nonumber\\ &&+\sum_{p=1}^{n}(-1)^{p+1+\mu+\nu'+\mu(\mu'+\nu-1)}\skew3\hat{f}_{\nu+k_p}^+ \prod_{i=1;i\neq p}^\nu\skew3\hat{f}_{k_i}\prod_{i=1}^{\mu'}\skew3\hat{f}_{t_i}\prod_{i=1}^\mu\skew3\hat{f}_{\nu+j_i}\prod_{i=1}^{\nu'}\skew3\hat{f}_{\nu+v_i} \nonumber\\ &&+\left .\sum_{p=1}^{\mu'}(-1)^{p+\mu+\nu+\nu'+\mu(m'-1+\nu)}\skew3\hat{f}_{\nu+t_p}^+ \prod_{i=1}^\nu\skew3\hat{f}_{k_i}\prod_{i=1;i\neq p}^{\mu'}\skew3\hat{f}_{t_i}\prod_{i=1}^\mu\skew3\hat{f}_{\nu+j_i}\prod_{i=1}^{\nu'}\skew3\hat{f}_{\nu+v_i} \right)\nonumber\\ &&\times\skew3\hat{\mathcal{P}}^{\mu'+\nu'+\mu+\nu}+{\rm rest}. \end{eqnarray} \tag{ B.8 }$$

In order for the first and third fermionic closeness conditions to be satisfied we have to demand that the explicit sums in (B.7) vanish. This gives us the following conditions:

$$\begin{eqnarray}&&\fl \sum_{\mu+\nu'=n_x}\sum_{\nu+\mu'=n_y}\sum_{P\underline{x},P\underline{y}}\mu(-1)^{1+\mu+(\mu-1)(\mu'+\nu)+\nu'}\epsilon_{\underline{x}}\epsilon_{\underline{y}}\Bigg((-1)^{\lfloor \frac{\mu+\nu}{2}\rfloor} G^{(\mu,\nu;\mu',\nu')}_{\underline{j},\underline{k};\underline{t},\underline{v}}\nonumber\\ &&+(-1)^{\lfloor \frac{\mu'+\nu'}{2}\rfloor+1 + \nu \nu' + \mu' \nu + \mu' \nu' + \mu \mu' + \mu \nu + \mu \nu'}\bar{G}^{(\nu,\mu;\nu',\mu')}_{\underline{k},\underline{j};\underline{v},\underline{t}} \Bigg)=0, \end{eqnarray} \tag{ B.9 }$$

where $\underline {y}=(j_2,j_3,\ldots ,j_{m},v_1,v_2,\ldots , v_{n'})$, $\underline {y}=(\underline {k},\underline {t})$, and $\epsilon _{\underline {x}}~(\epsilon _{\underline {y}})$ is a completely antisymmetric tensor. Equations (B.9) have to be satisfied for any variation of $\underline {x}$ and $\underline {y}$ and any j₁ and are in general not sufficient to ensure the closeness of the hierarchy. However they are sufficient in the case of general quadratic Markovian fermionic noise determined by the coupling operator basis L_j,k∈{c_jc_k,c^†_jc^†_k,c^†_jc_k}. In this case equations (B.9) simplify to

$$\begin{eqnarray} &&\fl 0= 4(G^{(2,0;2,0)}_{j_1,j_2;t_1,t_2}-\bar{G}^{(0,2;0,2)}_{j_1,j_2;t_1,t_2})+ G^{(1,1;1,1)}_{j_1,j_2;t_1,t_2}-G^{(1,1;1,1)}_{j_1,t_1;j_2,t_2}-\bar{G}^{(1,1;1,1)}_{j_2,j_1;t_2,t_1}+\bar{G}^{(1,1;1,1)}_{t_1,j_1;t_2,j_2}, \nonumber\\ &&\fl 0= -2(G^{(1,1;2,0)}_{j_1,k_1;t_1,t_2}+G^{(1,1;2,0)}_{j_1,t_1;k_1,t_2}+G^{(1,1;2,0)}_{j_1,t_2;k_1,t_1})+2(-\bar{G}^{(1,1;0,2)}_{k_1,j_1;t_1,t_2}+G^{(1,1;2,0)}_{t_1,j_1;k_1,t_2}+G^{(1,1;2,0)}_{t_2,j_1;t_1,k_1}),\nonumber\\ &&+2G^{(1,1;0,2)}_{j_1,k_1;t_2,t_1}+2\bar{G}^{(1,1;2,0)}_{k_1,j_1;t_2,t_1}, \nonumber\\ &&\fl 0= G^{(2,0;0,2)}_{j_1,j_2;t_1,t_2}+G^{(2,0;0,2)}_{j_1,t_1;t_2,j_2}+G^{(2,0;0,2)}_{j_1,t_2;j_2,t_1}-\bar{G}^{(0,2;2,0)}_{j_1,j_2;t_1,t_2}-\bar{G}^{(0,2;2,0)}_{j_1,t_1;t_2,j_2}-\bar{G}^{(0,2;2,0)}_{j_1,t_2;j_2,t_1}, \end{eqnarray} \tag{ B.10 }$$

and have to be satisfied for all different values of j₁,j₂,k₁,t₁,t₂.

The 'higher-order' requirements can be obtained in a similar manner by writing the vanishing conditions for unwanted terms in the dissipators.

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:

$$\begin{eqnarray} {\bf G}=\left(\begin{array}{@{}cccc@{}} G_{11} & 0 & G_{13} & 0\\ 0 & G_{22} & 0 & G_{24}\\ \bar{G}_{13} & 0 & G_{33} & 0\\ 0 & \bar{G}_{24} & 0 & G_{44} \end{array}\right), \end{eqnarray} \tag{ B.11 }$$

which for |G₂₄|² < |G₂₂G₄₄| and |G₁₃|² < |G₁₁G₃₃| 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₁₁ = G₂₂ and $\bar {G}_{13}=G_{24}$. For simplicity we write G₃₃ = G₄₄ + 2δ. In this case the dissipator may be written as

$$\begin{equation} \fl \skew3\hat{\mathcal{D}}_{\rm mixed}=G_{11}(\skew3\hat{f}_{1}^+\skew3\hat{f}_1-\skew3\hat{f}_2^+\skew3\hat{f}_2)-\delta(\hat{b}_1^+\hat{b}_1+\hat{b}^+_2\hat{b}_2)+G_{44}\hat{b}_1^+\hat{b}_2^++\sqrt{2}(G_{13}\hat{b}^+_1\skew3\hat{f}^+_2+\bar{G}_{13}\hat{b}^+_2\skew3\hat{f}^+_1)\skew3\hat{\mathcal{P}} \end{equation} \tag{ B.12 }$$

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.